D'oh (on the rounding)!
But on mcmcsamp(), I'm still confused. I changed the example slightly
and got the same problem:
> y <- (1:20)*pi
> x <- (1:20)^2
> group <- rep (1:2, each=10)
> M1 <- lmer (y ~ 1 + (1 | group))
> M2 <- lmer (y ~ 1 + x + (1 + x | group))
> mcmcsamp (M1, saveb=TRUE)
> mcmcsamp (M2, saveb=TRUE)
Error: inconsistent degrees of freedom and dimension
Error in t(.Call("mer_MCMCsamp", object, saveb, n, trans, PACKAGE =
"Matrix")) :
unable to find the argument 'x' in selecting a method for
function 't'
It really should be able to work (actually, the earlier example should
work too), but maybe it gets hung up when there are only two groups.
Indeed, when I change 20 and 2 above to 30 and 3, it works fine.
So I guess, as a practical matter, this is fine. The domain where
lmer() excels is large datasets with many groups; conversely, Bugs works
best with small datasets with few groups. However, I do like to use
lmer() as a starting point, so I hope that at some point it will fully
work in the above example also.
Also, since I have you on the line, so to speak, I noticed that coef()
gives estimated group-level coefficients, and ranef() gives these
coefficients centered at zero:
> coef(M2)
$group
(Intercept) x
1 6.885045 0.2701600
2 23.586015 0.1010110
> ranef(M2)
An object of class "lmer.ranef"
[[1]]
(Intercept) x
1 -8.350485 0.08457451
2 8.350485 -0.08457451
I was just wondering: is one or the other of these parameterizations
preferred by Doug Bates et al.? I wanted to know because we discuss
lmer() in our book, and I'd like our examples to remain relevant after
the book appears and lmer() continues to be developed.
Peter Dalgaard wrote:
>Andrew Gelman <[EMAIL PROTECTED]> writes:
>
>
>
>>I've been trying to keep track with lmer, and now I have a couple of
>>questions with the latest version of Matrix (0.995-4). I fit 2 very
>>similar models, and the results are severely rounded in one case and
>>rounded not at all in the other.
>>
>> > y <- 1:10
>> > group <- rep (c(1,2), c(5,5))
>> > M1 <- lmer (y ~ 1 + (1 | group))
>> > coef(M1)
>>$group
>> (Intercept)
>>1 3.1
>>2 7.9
>>
>> > x <- rep (c(1,2), c(3,7))
>> > M2 <- lmer (y ~ 1 + x + (1 + x | group))
>> > coef(M2)
>>$group
>> (Intercept) x
>>1 -0.755102 2.755102
>>2 0.616483 3.640738
>>
>>I can't figure out why everything is rounded for the first model but not
>>for the second. Also, mcmcsamp() works for M1 but not for M2:
>>
>>
>
>
>Well,
>
>
>
>>dput(coef(M1)[[1]])
>>
>>
>structure(list("(Intercept)" = c(3.10000000006436, 7.89999999993564
>)), .Names = "(Intercept)", row.names = c("1", "2"), class =
>"data.frame")
>
>
>>c(3.10000000006436, 7.89999999993564)
>>
>>
>[1] 3.1 7.9
>
>I.e., if you pass fake data, sometimes you get *results* that can be
>rounded to a few significant digits. R tries to get rid of trailing
>zeros in its print routines.
>
>
>
>
>> > mcmcsamp(M1)
>> (Intercept) log(sigma^2) log(grop.(In))
>>
>>
> > [1,] 9.099073 0.5711817 3.246981
>
>
>>attr(,"mcpar")
>>[1] 1 1 1
>>attr(,"class")
>>[1] "mcmc"
>>
>> > mcmcsamp(M2)
>>Error: inconsistent degrees of freedom and dimension
>>Error in t(.Call("mer_MCMCsamp", object, saveb, n, trans, PACKAGE =
>>"Matrix")) :
>> unable to find the argument 'x' in selecting a method for
>>function 't'
>>
>>
>
>Looks like a buglet, but
>
>
>
>>x
>>
>>
> [1] 1 1 1 2 2 2 2 2 2 2
>
>
>>group
>>
>>
> [1] 1 1 1 1 1 2 2 2 2 2
>
>Effects of x can (seemingly?) only be detected within group 1. I.e.
>the random variation of the effect of x is based on a sample of size
>1, so I'm actually more surprised that you get a fit at all...
>
>
>
--
Andrew Gelman
Professor, Department of Statistics
Professor, Department of Political Science
[EMAIL PROTECTED]
www.stat.columbia.edu/~gelman
Tues, Wed, Thurs:
Social Work Bldg (Amsterdam Ave at 122 St), Room 1016
212-851-2142
Mon, Fri:
International Affairs Bldg (Amsterdam Ave at 118 St), Room 711
212-854-7075
Mailing address:
1255 Amsterdam Ave, Room 1016
Columbia University
New York, NY 10027-5904
212-851-2142
(fax) 212-851-2164
______________________________________________
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
