MCMCglmm
Hi, In order to update the covariance matrix it is much easier if every combination of phase and marker exist. MCMCglmm generates any missing combinations and treats the unknown responses as missing data. It is just a computational strategy and the warning message can be ignored (I may suppress it in future versions). The error message is unrelated. I presume phase is a factor (?) with n levels. At some iteration the nXn covariance matrix associated with us(0+phase):marker becomes singular (or close to). If n=2 the error message implies that a variance has hit zero, or a correlation has hit -1 and 1. If n>2 then this implies that one (or more) eigenvalues of the covariance matrix has hit zero. Numerical problems arise when these conditions occur so MCMCglmm terminates. In your analysis you have used the default flat priors, but if a proper prior is specified these conditions do not generally arise. You can choose from the standard inverse-Wishart prior or from the parameter expanded non- central F prior (see Gelman 2006 Bayesian Analysis 1 515-533, or the CourseNotes). The latter is particularly useful if the variances are close to zero because you can get dramatic improvements in mixing. Cheers, Jarrod
On 9 Apr 2010, at 13:03, Iasonas Lamprianou wrote:
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.
The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336.