[R-sig-ME] lmer objects and mcmcsamp
Austin Frank
austin.frank at gmail.com
Thu Feb 26 23:23:43 CET 2009
Hugh--
Interesting question. I just tried it out and replicated the report on
R version 2.8.1 (2008-12-22) powerpc-apple-darwin8.11.1 with
lme4_0.999375-28 and Matrix_0.999375-21. I don't have the first idea
what's going on here, but I'm quite certain this question should also go
to R-sig-mixed-models, so I've copied that group in this response.
Thanks for the report,
/au
On Thu, Feb 26 2009, Hugh Rabagliati wrote:
> This is a bit of a general question, but it arose out of using the
> pvals.fnc function so I figured this forum might have a few ideas
> about it.
>
> The issue is whether there's a very odd bug in the mcmcsamp package of
> lme4? (or whether I still don't understand mcmc methods well enough).
>
> If I create a mer object using lmer, use it as an argument for
> mcmcsamp (sampling > 1 times), assign the output to a new mermcmc
> object and then examine my mer object again, I notice a rather
> peculiar thing. In particular, all of the variance/standard error
> terms change, as do the associated t values for the fixed effects.
> The estimated coefficients are unaffected. I figure this is a bug,
> because I can't see any reason why mcmcsamp would want to do this. I
> took a look through the code for mcmcsamp, but I don't speak C and
> nothing jumped out at me. Certainly, the code is only supposed to
> return a mermcmc object, so I have no idea why its messing with my
> mer.
>
> I've got this same result on a couple of different computers (all
> macs). It does, however, seem to be specific to either the version of
> lmer ( 0.999375-28) or of R (2.8.1) that I'm using. I tried it on an
> old PC version of R (2.5.1) using lme4 version 0.99875-9, and the
> same problems don't happen then. I've included the output from both
> the PC and mac versions below.
>
> Sage advice, comments, or confirmation that this is a known bug (I
> couldn't find mention of it elsewhere) would be welcome.
>
> Thanks very much,
>
> Hugh Rabagliati
>
>
> ########
> #
> #PC Code & output - R v. 2.5.1 & lme4 v. 0.99875-9.
> #
> (fm1 <- lmer(Reaction ~ Days + (1|Subject) + (0+Days|Subject),
> sleepstudy))
> #Linear mixed-effects model fit by REML
> #Formula: Reaction ~ Days + (1 | Subject) + (0 + Days | Subject)
> # Data: sleepstudy
> # AIC BIC logLik MLdeviance REMLdeviance
> # 1752 1764 -871.8 1752 1744
> #Random effects:
> # Groups Name Variance Std.Dev.
> # Subject (Intercept) 627.508 25.0501
> # Subject Days 35.858 5.9881
> # Residual 653.590 25.5654
> #number of obs: 180, groups: Subject, 18; Subject, 18
>
> #Fixed effects:
> # Estimate Std. Error t value
> #(Intercept) 251.405 6.885 36.51
> #Days 10.467 1.560 6.71
>
> #Correlation of Fixed Effects:
> # (Intr)
> #Days -0.184
>
> fm1 -> fm1.old
> samp0 <- mcmcsamp(fm1, n = 1000)
> fm1
> #Linear mixed-effects model fit by REML
> #Formula: Reaction ~ Days + (1 | Subject) + (0 + Days | Subject)
> # Data: sleepstudy
> # AIC BIC logLik MLdeviance REMLdeviance
> # 1752 1764 -871.8 1752 1744
> #Random effects:
> # Groups Name Variance Std.Dev.
> # Subject (Intercept) 627.508 25.0501
> # Subject Days 35.858 5.9881
> # Residual 653.590 25.5654
> #number of obs: 180, groups: Subject, 18; Subject, 18
>
> #Fixed effects:
> # Estimate Std. Error t value
> #(Intercept) 251.405 6.885 36.51
> #Days 10.467 1.560 6.71
>
> #Correlation of Fixed Effects:
> # (Intr)
> #Days -0.184
>
> #
> # As you can see, the estimates don't change here
> #
> ###########
>
>
> #########
> #
> # Mac R v. 2.8.1, ‘lme4’ version 0.999375-28
> #
> # I make two mers, which should be identical (I make two because there
> are some assignment weirdnesses going on here too, #which I haven't
> yet understood)
>
> (fm1.to_mcmc <- lmer(Reaction ~ Days + (1|Subject) + (0+Days|
> Subject), sleepstudy))
>
> #Linear mixed model fit by REML
> #Formula: Reaction ~ Days + (1 | Subject) + (0 + Days | Subject)
> # Data: sleepstudy
> # AIC BIC logLik deviance REMLdev
> # 1754 1770 -871.8 1752 1744
> #Random effects:
> # Groups Name Variance Std.Dev.
> # Subject (Intercept) 627.568 25.0513
> # Subject Days 35.858 5.9882
> # Residual 653.584 25.5653
> #Number of obs: 180, groups: Subject, 18
>
> #Fixed effects:
> # Estimate Std. Error t value
> #(Intercept) 251.405 6.885 36.51
> #Days 10.467 1.559 6.71
>
> #Correlation of Fixed Effects:
> # (Intr)
> #Days -0.184
>
>
> (fm1.not_to_mcmc <- lmer(Reaction ~ Days + (1|Subject) + (0+Days|
> Subject), sleepstudy))
>
> #Linear mixed model fit by REML
> #Formula: Reaction ~ Days + (1 | Subject) + (0 + Days | Subject)
> # Data: sleepstudy
> # AIC BIC logLik deviance REMLdev
> # 1754 1770 -871.8 1752 1744
> #Random effects:
> # Groups Name Variance Std.Dev.
> # Subject (Intercept) 627.568 25.0513
> # Subject Days 35.858 5.9882
> # Residual 653.584 25.5653
> #Number of obs: 180, groups: Subject, 18
>
> #Fixed effects:
> # Estimate Std. Error t value
> #(Intercept) 251.405 6.885 36.51
> #Days 10.467 1.559 6.71
>
> #Correlation of Fixed Effects:
> # (Intr)
> #Days -0.184
>
>
> samp0 <- mcmcsamp(fm1.to_mcmc, n = 1000)
> fm1.to_mcmc
>
> #Linear mixed model fit by REML
> #Formula: Reaction ~ Days + (1 | Subject) + (0 + Days | Subject)
> # Data: sleepstudy
> # AIC BIC logLik deviance REMLdev
> # 1763 1779 -876.7 1761 1753
> #Random effects:
> # Groups Name Variance Std.Dev.
> # Subject (Intercept) 868.398 29.469
> # Subject Days 49.619 7.044
> #Residual 904.398 30.073
> #Number of obs: 180, groups: Subject, 18
>
> #Fixed effects:
> # Estimate Std. Error t value
> #(Intercept) 251.405 5.260 47.79
> #Days 10.467 1.518 6.90
>
> # All the variances etc are different from before
>
> #Correlation of Fixed Effects:
> # (Intr)
> #Days -0.343
>
> fm1.not_to_mcmc
>
> #Linear mixed model fit by REML
> #Formula: Reaction ~ Days + (1 | Subject) + (0 + Days | Subject)
> # Data: sleepstudy
> # AIC BIC logLik deviance REMLdev
> # 1754 1770 -871.8 1752 1744
> #Random effects:
> # Groups Name Variance Std.Dev.
> # Subject (Intercept) 627.568 25.0513
> # Subject Days 35.858 5.9882
> # Residual 653.584 25.5653
> #Number of obs: 180, groups: Subject, 18
>
> #Fixed effects:
> # Estimate Std. Error t value
> #(Intercept) 251.405 6.885 36.51
> #Days 10.467 1.559 6.71
>
> # Variances are unaffected here
>
> #Correlation of Fixed Effects:
> # (Intr)
> #Days -0.184
--
Austin Frank
http://aufrank.net
GPG Public Key (D7398C2F): http://aufrank.net/personal.asc
More information about the R-sig-mixed-models
mailing list