[R-sig-ME] lmer model and mcmcsamp

[Jarrod Hadfield]

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.