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