[R] Suspicious output from lme4-mcmcsamp

[De Woody J.A.]

Hello, R community,

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.

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


> 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?

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