[R-sig-ME] Fwd: FW: Suspicious output from lme4-mcmcsamp

[Douglas Bates]

I meant to cc: this reply to the list.


---------- Forwarded message ----------
From: Douglas Bates <bates at stat.wisc.edu>
Date: Wed, Oct 8, 2008 at 12:27 PM
Subject: Re: [R-sig-ME] FW: Suspicious output from lme4-mcmcsamp
To: "De Woody J.A." <j.dewoody at soton.ac.uk>


On Wed, Oct 8, 2008 at 9:41 AM, De Woody J.A. <j.dewoody at soton.ac.uk> wrote:
> Hello, lme4 community,
>
> Sorry for the re-post, but I hadn't realized there was a list specific to lme4 when I posted to the R-help earlier.
>
> I have been using the lmer and mcmcsamp functions in R with some difficulty.  I do not believe this is my code or data, however, because my attempts to use the sample code and 'sleepstudy' data provided with the lme4 packaged (and used on several R-Wiki pages) do not return the same results as those indicated in the help pages.  For instance:
>
>> sessionInfo()
> R version 2.7.2 (2008-08-25)
> i386-pc-mingw32
>
> locale:
> LC_COLLATE=English_United Kingdom.1252;LC_CTYPE=English_United Kingdom.1252;LC_MONETARY=English_United Kingdom.1252;LC_NUMERIC=C;LC_TIME=English_United Kingdom.1252
>
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base
>
> other attached packages:
> [1] lme4_0.999375-26   Matrix_0.999375-11 lattice_0.17-13
>
> loaded via a namespace (and not attached):
> [1] grid_2.7.2
>
>> fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
>> sm1 <- mcmcsamp(fm1, 5000)
> Error in .local(object, n, verbose, ...) :
>  Code for non-trivial theta_T not yet written
>
> ##
> I cannot find exactly what this theta_T error means, although I do find it mentioned in what I believe to be source code.  Regardless, I cannot understand why the mcmcsamp returns the error for this data set.

As Dieter Menne indicated, this is one of the error messages that is
meant more as a reminder to myself that I should write the code for
this case.  The rather crytic terminology relates to the particular
way that the covariance matrix of the random effects is parameterized,
as sigma^2 * T %*% S %*% S %*% t(T) where S is a diagonal (S)cale
matrix with non-negative diagonal elements and T is a unit,
lower-(T)riangular matrix.  If the random effects are uncorrelated
then T is the identity matrix and there is no need to sample from the
posterior distribution of the parameters determining T.  I haven't yet
written the code to sample from that posterior distribution when the
random effects are correlated.

> Even when I change the model and the mcmcsamp appears to run, the output is not as expected:

Expected by whom?  I don't know of examples where summary of an
merMCMC object has had a form other than what is shown below.

Such objects are in an S4 class and, in the fullness of time, as Bill
Venables is fond of saying, there will be extractor methods for them
that will provide more meaningful output.  At present there is not a
lot of useful output.  Mostly I look at the plots produced by

xyplot(sm2)

to see how the chains evolved.  At present, even that has a few
problems for models with multiple, uncorrelated random effects per
level of the grouping factor.

>
>> fm2 <- lmer(Reaction ~ Days + (1|Subject), sleepstudy)
>> sm2 <- mcmcsamp(fm2, 5000)
>> summary(sm2)
>  Length   Class    Mode
>      1 merMCMC      S4
>> str(sm2)
> Formal class 'merMCMC' [package "lme4"] with 9 slots
>  ..@ Gp      : int [1:2] 0 18
>  ..@ ST      : num [1, 1:5000] 1.198 0.932 0.835 0.826 0.933 ...
>  ..@ call    : language lmer(formula = Reaction ~ Days + (1 | Subject), data = sleepstudy)
>  ..@ deviance: num [1:5000] 1794 1794 1796 1798 1798 ...
>  ..@ dims    : Named int [1:17] 1 180 2 18 1 1 1 2 5 1 ...
>  .. ..- attr(*, "names")= chr [1:17] "nf" "n" "p" "q" ...
>  ..@ fixef   : num [1:2, 1:5000] 251.4  10.5 253.3  11.0 259.5 ...
>  .. ..- attr(*, "dimnames")=List of 2
>  .. .. ..$ : chr [1:2] "(Intercept)" "Days"
>  .. .. ..$ : NULL
>  ..@ nc      : int 1
>  ..@ ranef   : num[1:18, 0 ]
>  ..@ sigma   : num [1, 1:5000] 31.0 29.7 30.4 28.4 38.1 ...
>
> ##
> As I understand it, the call >summary(sm2) should return information of the results of the mcmcsamp distribution.    In addition, I am expecting the str(sm2) to show the 'fixef' slot to have something resembling "log(sigma^2)" and "log(Subject.(In))".  Am I wrong?  Are all of the outputs in the correct form?

I don't understand why the fixef slot would have log(sigma^2), etc.
The fixef slot is the values of the fixed-effects parameters in the
sample.

The sigma slot is the sample from the distribution of the common scale
parameter, which is written as sigma.  The ST slot is the sample from
the distribution of the relative standard deviation of the random
effects.  To get the standard deviation of the random effects you
multiply ST by sigma.

I am not confident of the results from mcmcsamp at present and
regrettably I won't have time to look at it in more detail for a
couple of months.  If it were better to avoid confusion for the time
being I could disable it.  If someone else wants to experiment with it
send me email off-list and I will outline what should be happending in
that code.

>
> Has anyone else had this problem?  Could this be related to the possible 'mistake in the mcmcsamp function at present' mentioned in the recent postings regarding the $ST and $sigma slots (Re: mcmcsamp(lme4): What is contained in $ST and $sigma?)?
>
> Any thoughts, suggestions, or directions would, of course, be most appreciated.
>
> Many thanks!
> Jenn
>
> ****************************
> Jennifer DeWoody
> University of Southampton
> School of Biological Sciences
> Building 62, Room 6007, Boldrewood Campus Southampton  SO16 7PX United Kingdom
> Voice: +44 (0)23 8059 4286
> Email: j.dewoody at soton.ac.uk
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>