MCMCglmm
Dear Jarrod Hadfield I followed your advice but got this message: m5_06.mcmc<- MCMCglmm(score ~ 1, random=~us(0+phase):marker+candidate+batch, data=mg2006_sub) Warning in MCMCglmm(score ~ 1, random = ~us(0 + phase):marker + candidate + : some combinations in us(0 + phase):marker do not exist and 118 missing records have been generated Error in MCMCglmm(score ~ 1, random = ~us(0 + phase):marker + candidate + : ill-conditioned G/R structure: use proper priors if you haven't or rescale data if you have could you please help me? Dr. Iasonas Lamprianou Assistant Professor (Educational Research and Evaluation) Department of Education Sciences European University-Cyprus P.O. Box 22006 1516 Nicosia Cyprus Tel.: +357-22-713178 Fax: +357-22-590539 Honorary Research Fellow Department of Education The University of Manchester Oxford Road, Manchester M13 9PL, UK Tel. 0044 161 275 3485 iasonas.lamprianou at manchester.ac.uk
--- On Thu, 8/4/10, Jarrod Hadfield <j.hadfield at ed.ac.uk> wrote:
From: Jarrod Hadfield <j.hadfield at ed.ac.uk> Subject: Re: [R-sig-ME] MCMCglmm To: "Iasonas Lamprianou" <lamprianou at yahoo.com> Cc: r-sig-mixed-models at r-project.org Date: Thursday, 8 April, 2010, 15:10 Dear Jason, These should work: m4.mcmc <- MCMCglmm(score ~ 1, random=~marker+marker:day+candidate +batch, data=mg2006) m5 .mcmc<- MCMCglmm(score ~ 1, random=~us(1+day):marker+candidate +batch, data=mg2006) You may need to code day as a factor for m4, and as numeric for m5? depending on the model you actually want to fit. Cheers, Jarrod On 8 Apr 2010, at 14:58, Iasonas Lamprianou wrote:
Dear all, I recently experimented with MCMCglmm and I loved
(really loved) the?
fact that it will give me confidence intervals for the
variance of?
the random effects. It seems that MCMC is a reasonable
method to do?
so, in contrast to REML which seems to have problems
on this front.?
However, MCMCglmm is painfully slower than lmer which
is more?
familiar to me. The good news is that the point
estimates of lmer?
are near the centre of the confidence intervals by
MCMCglmm.
I reduced my sample size a bit and managed to fit
those two models?
with lmer (the second would not fit because it needed
1.5GB or RAM).?
Both seem to have a reasonable fit (at least at first
look).
m4 <- lmer(score ~
1+(1|marker/day)+(1|candidate)+(1|batch), mg2006)
m5 <- lmer(score ~
1+(1+day|marker)+(1|candidate)+(1|batch), mg2006)
I would like to run these two models above with
MCMCglmm. Does?
anyone know how to do it? Thank you for the help Jason Dr. Iasonas Lamprianou Assistant Professor (Educational Research and
Evaluation)
Department of Education Sciences European University-Cyprus P.O. Box 22006 1516 Nicosia Cyprus Tel.: +357-22-713178 Fax: +357-22-590539 Honorary Research Fellow Department of Education The University of Manchester Oxford Road, Manchester M13 9PL, UK Tel. 0044? 161 275 3485 iasonas.lamprianou at manchester.ac.uk --- On Thu, 8/4/10, r-sig-mixed-models-request at r-project.org
<r-sig-mixed-models-request at r-project.org
wrote:
From: r-sig-mixed-models-request at r-project.org
<r-sig-mixed-models-request at r-project.org
Subject: R-sig-mixed-models Digest, Vol 40, Issue
15
To: r-sig-mixed-models at r-project.org Date: Thursday, 8 April, 2010, 12:13 Send R-sig-mixed-models mailing list submissions to ? ???r-sig-mixed-models at r-project.org To subscribe or unsubscribe via the World Wide
Web, visit
? ???https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models or, via email, send a message with subject or body
'help'
to ? ???r-sig-mixed-models-request at r-project.org You can reach the person managing the list at ? ???r-sig-mixed-models-owner at r-project.org When replying, please edit your Subject line so it
is more
specific than "Re: Contents of R-sig-mixed-models
digest..."
Today's Topics: ? ? 1. Re: Multi-level models Odds ratio
(E
T) ? ? 2. Re: Multi-level models Odds ratio (ONKELINX, Thierry) ? ? 3. Re: Multi-level models Odds ratio
(E
T) ? ? 4. Re: Multi-level models Odds ratio (Andy Fugard (Work)) ? ? 5. Re: Multi-level models Odds ratio
(E
T)
----------------------------------------------------------------------
Message: 1 Date: Thu, 8 Apr 2010 11:27:49 +0100 From: E T <2nuzzbot at gmail.com> To: Daniel Ezra Johnson <danielezrajohnson at gmail.com> Cc: "r-sig-mixed-models at r-project.org" ? ???<r-sig-mixed-models at r-project.org> Subject: Re: [R-sig-ME] Multi-level models Odds
ratio
Message-ID: ? ???<l2w706f8d1f1004080327re8708f46mc7d334f41ae19a10 at mail.gmail.com> Content-Type: text/plain odds.ratios = exp(coefs(model)) Thanks, however unfortunately when I try the above
command
I receive the following error: Error: could not find function "coefs" Regards Et On Wed, Apr 7, 2010 at 5:47 PM, Daniel Ezra
Johnson <
danielezrajohnson at gmail.com> wrote:
something like odds.ratios =
exp(coefs(model))
On Apr 7, 2010, at 12:28 PM, E T <2nuzzbot at gmail.com>
wrote:
???Hi all,
Apologies for the simplicity of my
question....
however any advice is
greatly appreciated. Thanks Is there a specific command available to
obtain
the odds ratios produced
from a multilevel logistic model? I have estimated a multi-level logistic
model
using the lme4 package. I
can obtain results using the 'summary'
command,
however I would like to obtain
the computed odds ratios. (Similar to the output that can be
produced for
logistic GLM using the
logistic.display command from the epicalc
package).
? ? [[alternative HTML version
deleted]]
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models ? ???[[alternative HTML version deleted]] ------------------------------ Message: 2 Date: Thu, 8 Apr 2010 12:32:01 +0200 From: "ONKELINX, Thierry" <Thierry.ONKELINX at inbo.be> To: "E T" <2nuzzbot at gmail.com>, "Daniel Ezra Johnson" ? ???<danielezrajohnson at gmail.com> Cc: r-sig-mixed-models at r-project.org Subject: Re: [R-sig-ME] Multi-level models Odds ratio Message-ID: ? ???<2E9C414912813E4EB981326983E0A104071B69A6 at inboexch.inbo.be> Content-Type: text/plain; charset="us-ascii" It should be exp(coef(model)) Without the "s" HTH, Thierry ------------------------------------------------------------------------ ---- ir. Thierry Onkelinx Instituut voor natuur- en bosonderzoek team Biometrie & Kwaliteitszorg Gaverstraat 4 9500 Geraardsbergen Belgium Research Institute for Nature and Forest team Biometrics & Quality Assurance Gaverstraat 4 9500 Geraardsbergen Belgium tel. + 32 54/436 185 Thierry.Onkelinx at inbo.be www.inbo.be To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of. ~ Sir Ronald Aylmer Fisher The plural of anecdote is not data. ~ Roger Brinner The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data. ~ John Tukey -----Oorspronkelijk bericht----- Van: r-sig-mixed-models-bounces at r-project.org [mailto:r-sig-mixed-models-bounces at r-project.org] Namens E T Verzonden: donderdag 8 april 2010 12:28 Aan: Daniel Ezra Johnson CC: r-sig-mixed-models at r-project.org Onderwerp: Re: [R-sig-ME] Multi-level models Odds ratio odds.ratios = exp(coefs(model)) Thanks, however unfortunately when I try the above command I receive the following error: Error: could not find function "coefs" Regards Et On Wed, Apr 7, 2010 at 5:47 PM, Daniel Ezra Johnson < danielezrajohnson at gmail.com> wrote: something like odds.ratios = exp(coefs(model)) On Apr 7, 2010, at 12:28 PM, E T <2nuzzbot at gmail.com> wrote: ???Hi all, Apologies for the simplicity of my question.... however any advice is greatly appreciated. Thanks Is there a specific command available to obtain the odds ratios produced from a multilevel logistic model? I have estimated a multi-level logistic model using the lme4 package. I can obtain results using the 'summary' command, however I would like to obtain the computed odds ratios. (Similar to the output that can be produced for logistic GLM using the logistic.display command from the epicalc package). ? ? [[alternative HTML version deleted]] _______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models ? ???[[alternative HTML version deleted]] _______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models Druk dit bericht a.u.b. niet onnodig af. Please do not print this message unnecessarily. Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer en binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd is door een geldig ondertekend document. The views expressed in? this message and any annex are purely those of the writer and may not be regarded as stating an official position of INBO, as long as the message is not confirmed by a duly signed document. ------------------------------ Message: 3 Date: Thu, 8 Apr 2010 11:35:03 +0100 From: E T <2nuzzbot at gmail.com> To: Daniel Ezra Johnson <danielezrajohnson at gmail.com> Cc: "r-sig-mixed-models at r-project.org" ? ???<r-sig-mixed-models at r-project.org> Subject: Re: [R-sig-ME] Multi-level models Odds ratio Message-ID: ? ???<i2y706f8d1f1004080335q84f61b78u8c7b656b67a08a8e at mail.gmail.com> Content-Type: text/plain If I use the command coef(model) this extracts the coefficients in the model, however if I try exp(coef(model)) I receive an error: Error in exp(coef(model)) : Non-numeric argument to mathematical function I could manually get the exp of each factor in my model..... but as I have a large model (and also have numerous other models to produce), I was wondering if there was an automated method Regards Et On Thu, Apr 8, 2010 at 11:27 AM, E T <2nuzzbot at gmail.com> wrote: odds.ratios = exp(coefs(model)) Thanks, however unfortunately when I try the above command I receive the following error: Error: could not find function "coefs" Regards Et On Wed, Apr 7, 2010 at 5:47 PM, Daniel Ezra Johnson < danielezrajohnson at gmail.com> wrote: something like odds.ratios = exp(coefs(model)) On Apr 7, 2010, at 12:28 PM, E T <2nuzzbot at gmail.com> wrote: ???Hi all, Apologies for the simplicity of my question.... however any advice is greatly appreciated. Thanks Is there a specific command available to obtain the odds ratios produced from a multilevel logistic model? I have estimated a multi-level logistic model using the lme4 package. I can obtain results using the 'summary' command, however I would like to obtain the computed odds ratios. (Similar to the output that can be produced for logistic GLM using the logistic.display command from the epicalc package). ? ? [[alternative HTML version deleted]] _______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models ? ???[[alternative HTML version deleted]] ------------------------------ Message: 4 Date: Thu, 08 Apr 2010 12:48:30 +0200 From: "Andy Fugard (Work)" <andy.fugard at sbg.ac.at> To: E T <2nuzzbot at gmail.com> Cc: "r-sig-mixed-models at r-project.org" ? ???<r-sig-mixed-models at r-project.org> Subject: Re: [R-sig-ME] Multi-level models Odds ratio Message-ID: <4BBDB47E.8030305 at sbg.ac.at> Content-Type: text/plain; charset=ISO-8859-1 Here's another example, borrowed from the help for "lmer": gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd), ? ? ? ? ? ? ???family = binomial, data = cbpp) As you say, coef works: coef(gm1) $herd ? ? (Intercept) period2???period3???period4 1???-0.8085096 -0.9923347 -1.128675 -1.580374 2???-1.6974292 -0.9923347 -1.128675 -1.580374 3???-0.9922697 -0.9923347 -1.128675 -1.580374 4???-1.3592525 -0.9923347 -1.128675 -1.580374 5???-1.5885461 -0.9923347 -1.128675 -1.580374 6???-1.7987950 -0.9923347 -1.128675 -1.580374 7???-0.5091313 -0.9923347 -1.128675 -1.580374 8???-0.7991613 -0.9923347 -1.128675 -1.580374 9???-1.6361848 -0.9923347 -1.128675 -1.580374 10? -1.9394614 -0.9923347 -1.128675 -1.580374 11? -1.4831632 -0.9923347 -1.128675 -1.580374 12? -1.4633469 -0.9923347 -1.128675 -1.580374 13? -2.0884474 -0.9923347 -1.128675 -1.580374 14? -0.4278151 -0.9923347 -1.128675 -1.580374 15? -1.9290041 -0.9923347 -1.128675 -1.580374 But note the "$herd" bit.? Since this model has a varying intercept by herd, you get a column in the resulting data frame called "herd". So you could try, for this example: exp(coef(gm1)$herd) ? ? (Intercept)???period2???period3???period4 1? ? 0.4455216 0.3707102 0.3234614 0.2058981 2? ? 0.1831538 0.3707102 0.3234614 0.2058981 3? ? 0.3707343 0.3707102 0.3234614 0.2058981 4? ? 0.2568527 0.3707102 0.3234614 0.2058981 5? ? 0.2042223 0.3707102 0.3234614 0.2058981 6? ? 0.1654982 0.3707102 0.3234614 0.2058981 7? ? 0.6010174 0.3707102 0.3234614 0.2058981 8? ? 0.4497060 0.3707102 0.3234614 0.2058981 9? ? 0.1947215 0.3707102 0.3234614 0.2058981 10???0.1437814 0.3707102 0.3234614 0.2058981 11???0.2269188 0.3707102 0.3234614 0.2058981 12???0.2314603 0.3707102 0.3234614 0.2058981 13???0.1238793 0.3707102 0.3234614 0.2058981 14???0.6519320 0.3707102 0.3234614 0.2058981 15???0.1452928 0.3707102 0.3234614 0.2058981 Since the slopes don't vary by herd, you might also want just the fixed effects: exp(fixef(gm1)) (Intercept)? ???period2 ? ? period3? ???period4 0.2469585???0.3707102???0.3234614???0.2058981 HTH, Andy E T wrote: If I use the command coef(model) this extracts the coefficients in the model, however if I try exp(coef(model)) I receive an error: Error in exp(coef(model)) : Non-numeric argument to mathematical function I could manually get the exp of each factor in my model..... but as I have a large model (and also have numerous other models to produce), I was wondering if there was an automated method Regards Et On Thu, Apr 8, 2010 at 11:27 AM, E T <2nuzzbot at gmail.com> wrote: odds.ratios = exp(coefs(model)) Thanks, however unfortunately when I try the above command I receive the following error: Error: could not find function "coefs" Regards Et On Wed, Apr 7, 2010 at 5:47 PM, Daniel Ezra Johnson < danielezrajohnson at gmail.com> wrote: something like odds.ratios = exp(coefs(model)) On Apr 7, 2010, at 12:28 PM, E T <2nuzzbot at gmail.com> wrote: ???Hi all, Apologies for the simplicity of my question.... however any advice is greatly appreciated. Thanks Is there a specific command available to obtain the odds ratios produced from a multilevel logistic model? I have estimated a multi-level logistic model using the lme4 package. I can obtain results using the 'summary' command, however I would like to obtain the computed odds ratios. (Similar to the output that can be produced for logistic GLM using the logistic.display command from the epicalc package). ? ? [[alternative HTML version deleted]] _______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models ? ???[[alternative HTML version deleted]] _______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models -- Andy Fugard, Postdoctoral researcher, ESF LogICCC project "Modeling human inference within the framework of probability logic" Department of Psychology, University of Salzburg, Austria http://www.andyfugard.info ------------------------------ Message: 5 Date: Thu, 8 Apr 2010 12:13:23 +0100 From: E T <2nuzzbot at gmail.com> To: "Andy Fugard (Work)" <andy.fugard at sbg.ac.at> Cc: "r-sig-mixed-models at r-project.org" ? ???<r-sig-mixed-models at r-project.org> Subject: Re: [R-sig-ME] Multi-level models Odds ratio Message-ID: ? ???<p2p706f8d1f1004080413k3014fe10q444382c927c2f90e at mail.gmail.com> Content-Type: text/plain exp(coef(model)$group) exp(fixef(model)) Thanks.... yes this worked successfully :o) Et On Thu, Apr 8, 2010 at 11:48 AM, Andy Fugard (Work) <andy.fugard at sbg.ac.at>wrote: Here's another example, borrowed from the help for "lmer": gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd), ? ? ? ? ? ? ???family = binomial, data = cbpp) As you say, coef works: coef(gm1) $herd ? ? (Intercept) period2???period3???period4 1???-0.8085096 -0.9923347 -1.128675 -1.580374 2???-1.6974292 -0.9923347 -1.128675 -1.580374 3???-0.9922697 -0.9923347 -1.128675 -1.580374 4???-1.3592525 -0.9923347 -1.128675 -1.580374 5???-1.5885461 -0.9923347 -1.128675 -1.580374 6???-1.7987950 -0.9923347 -1.128675 -1.580374 7???-0.5091313 -0.9923347 -1.128675 -1.580374 8???-0.7991613 -0.9923347 -1.128675 -1.580374 9???-1.6361848 -0.9923347 -1.128675 -1.580374 10? -1.9394614 -0.9923347 -1.128675 -1.580374 11? -1.4831632 -0.9923347 -1.128675 -1.580374 12? -1.4633469 -0.9923347 -1.128675 -1.580374 13? -2.0884474 -0.9923347 -1.128675 -1.580374 14? -0.4278151 -0.9923347 -1.128675 -1.580374 15? -1.9290041 -0.9923347 -1.128675 -1.580374 But note the "$herd" bit.? Since this model has a varying intercept by herd, you get a column in the resulting data frame called "herd". So you could try, for this example: exp(coef(gm1)$herd) ? ? (Intercept)???period2???period3???period4 1? ? 0.4455216 0.3707102 0.3234614 0.2058981 2? ? 0.1831538 0.3707102 0.3234614 0.2058981 3? ? 0.3707343 0.3707102 0.3234614 0.2058981 4? ? 0.2568527 0.3707102 0.3234614 0.2058981 5? ? 0.2042223 0.3707102 0.3234614 0.2058981 6? ? 0.1654982 0.3707102 0.3234614 0.2058981 7? ? 0.6010174 0.3707102 0.3234614 0.2058981 8? ? 0.4497060 0.3707102 0.3234614 0.2058981 9? ? 0.1947215 0.3707102 0.3234614 0.2058981 10???0.1437814 0.3707102 0.3234614 0.2058981 11???0.2269188 0.3707102 0.3234614 0.2058981 12???0.2314603 0.3707102 0.3234614 0.2058981 13???0.1238793 0.3707102 0.3234614 0.2058981 14???0.6519320 0.3707102 0.3234614 0.2058981 15???0.1452928 0.3707102 0.3234614 0.2058981 Since the slopes don't vary by herd, you might also want just the fixed effects: exp(fixef(gm1)) (Intercept)? ???period2 ? ? period3? ???period4 0.2469585???0.3707102???0.3234614???0.2058981 HTH, Andy E T wrote: If I use the command coef(model) this extracts the coefficients in the model, however if I try exp(coef(model)) I receive an error: Error in exp(coef(model)) : Non-numeric argument to mathematical function I could manually get the exp of each factor in my model..... but as I have a large model (and also have numerous other models to produce), I was wondering if there was an automated method Regards Et On Thu, Apr 8, 2010 at 11:27 AM, E T <2nuzzbot at gmail.com> wrote: odds.ratios = exp(coefs(model)) Thanks, however unfortunately when I try the above command I receive the following error: Error: could not find function "coefs" Regards Et On Wed, Apr 7, 2010 at 5:47 PM, Daniel Ezra Johnson < danielezrajohnson at gmail.com> wrote: something like odds.ratios = exp(coefs(model)) On Apr 7, 2010, at 12:28 PM, E T <2nuzzbot at gmail.com> wrote: ???Hi all, Apologies for the simplicity of my question.... however any advice is greatly appreciated. Thanks Is there a specific command available to obtain the odds ratios produced from a multilevel logistic model? I have estimated a multi-level logistic model using the lme4 package. I can obtain results using the 'summary' command, however I would like to obtain the computed odds ratios. (Similar to the output that can be produced for logistic GLM using the logistic.display command from the epicalc package). ? ? [[alternative HTML version deleted]] _______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models ? ? ? ? [[alternative HTML version deleted]] _______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models -- Andy Fugard, Postdoctoral researcher, ESF LogICCC project "Modeling human inference within the framework of probability logic" Department of Psychology, University of Salzburg, Austria http://www.andyfugard.info ? ???[[alternative HTML version deleted]] ------------------------------ _______________________________________________ R-sig-mixed-models mailing list R-sig-mixed-models at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models End of R-sig-mixed-models Digest, Vol 40, Issue 15 ************************************************** _______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models -- The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336.