An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20110811/47df18ee/attachment.pl>
lmer model and mcmcsamp
4 messages · Emma Jones, Cyrus Shaoul, Dennis Murphy +1 more
An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20110811/09b0b039/attachment.pl>
Hi: See http://glmm.wikidot.com/faq and pay particular attention to the sections labeled * Why doesn't lme4 provide p values/denominator degrees of freedom? What other options do I have? * Implementations of MCMC and parametric bootstrap That will give you a fairly broad view of the state of affairs at the present time. HTH, Dennis
On Thu, Aug 11, 2011 at 9:36 AM, Emma Jones <stp08emj at sheffield.ac.uk> wrote:
I posted this message back in February and got no response....but am still puzzled. Can anyone help? Thanks Emma Hi, I am new to using the package lme4 and wondered if anyone can help me? I have fitted a model (note the data is unbalanced)- see below with the results:
/ ?a=lmer(RW~1+(1|Year),chron,na.action=na.omit)
/>/ ?a /Linear mixed model fit by REML Formula: RW ~ 1 + (1 | Year) ? ?Data: chron ? AIC ?BIC logLik deviance REMLdev ?8771 8790 ?-4383 ? ? 8761 ? ?8765 Random effects: ?Groups ? Name ? ? ? ?Variance Std.Dev. ?Year ? ? (Intercept) 0.48050 ?0.69318 ?Residual ? ? ? ? ? ? 0.52704 ?0.72598 Number of obs: 3719, groups: Year, 239 Fixed effects: ? ? ? ? ? ? Estimate Std. Error t value (Intercept) -0.02793 ? ?0.04803 -0.5815 This I am happy with and I know is correctly fitted- my problems don't lie here. What I would ideally like now is to produce 95% credible intervals on the variance components. I have found some code on a help sheet online and followed this... samp=mcmcsamp(a,n=10000) #gives a credible interval for mcmc HPDinterval(VarCorr(samp,type="varcov"),prob=0.95) ? ? ? ? ?lower ? ? upper [1,] 0.1958900 0.2678522 [2,] 0.5317379 0.5853902 attr(,"Probability") 0.95 As it can be seen, both of my point estimates lie outside the credible interval. I am confused by this, can anyone advise?? Many thanks in advance, Emma ? ? ? ?[[alternative HTML version deleted]]
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
2 days later
Hi,
I have had similar problems - I came to the conclusion that a) I had
to make some transformation such as log something/divide by sigma etc
b) the "locally-uniform" priors are i) not uniform or ii) uniform but
in the wrong locality or c) mcmcsamp is broken.
It would be good to get some clarity; mcmcsamp even seems to behave
oddly for the example data set:
fm1 <- lmer(Reaction ~ Days + (1|Subject), sleepstudy)
fm1.mcmc<-mcmcsamp(fm1, n=1000)
summary(fm1)
Linear mixed model fit by REML
Formula: Reaction ~ Days + (1 | Subject)
Data: sleepstudy
AIC BIC logLik deviance REMLdev
1794 1807 -893.2 1794 1786
Random effects:
Groups Name Variance Std.Dev.
Subject (Intercept) 1378.18 37.124
Residual 960.46 30.991
Number of obs: 180, groups: Subject, 18
Fixed effects:
Estimate Std. Error t value
(Intercept) 251.4051 9.7459 25.80
Days 10.4673 0.8042 13.02
HPDinterval(VarCorr(fm1.mcmc,type="varcov"),prob=0.95)
lower upper
[1,] 307.8815 868.9589
[2,] 862.6868 1363.5824
attr(,"Probability")
[1] 0.95
Cheers,
Jarrod
On 11 Aug 2011, at 17:36, Emma Jones wrote:
I posted this message back in February and got no response....but am still puzzled. Can anyone help? Thanks Emma Hi, I am new to using the package lme4 and wondered if anyone can help me? I have fitted a model (note the data is unbalanced)- see below with the results:
/ a=lmer(RW~1+(1|Year),chron,na.action=na.omit)
/>/ a
/Linear mixed model fit by REML
Formula: RW ~ 1 + (1 | Year)
Data: chron
AIC BIC logLik deviance REMLdev
8771 8790 -4383 8761 8765
Random effects:
Groups Name Variance Std.Dev.
Year (Intercept) 0.48050 0.69318
Residual 0.52704 0.72598
Number of obs: 3719, groups: Year, 239
Fixed effects:
Estimate Std. Error t value
(Intercept) -0.02793 0.04803 -0.5815
This I am happy with and I know is correctly fitted- my problems don't
lie here.
What I would ideally like now is to produce 95% credible intervals on
the variance components. I have found some code on a help sheet online
and followed this...
samp=mcmcsamp(a,n=10000)
#gives a credible interval for mcmc
HPDinterval(VarCorr(samp,type="varcov"),prob=0.95)
lower upper
[1,] 0.1958900 0.2678522
[2,] 0.5317379 0.5853902
attr(,"Probability")
0.95
As it can be seen, both of my point estimates lie outside the credible
interval.
I am confused by this, can anyone advise??
Many thanks in advance,
Emma
[[alternative HTML version deleted]]
_______________________________________________ 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.