Suspicious output from lme4-mcmcsamp
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