[R] beginner's questions about lme, fixed and random effects
Peter Dalgaard BSA
p.dalgaard at biostat.ku.dk
Mon Dec 3 17:00:00 CET 2001
baron at cattell.psych.upenn.edu (Jonathan Baron) writes:
> I'm trying to understand better the differences between fixed and
> random effects by running very simple examples in the nlme
> package. My first attempt was to try doing a t-test in lme.
> This is very similar to the Rail example that comes with nlme,
> but it has two groups instead of five.
>
> So I try
>
> a1 <- 1:10
> a2 <- 7:16
> t.test(a2,a1)
>
> getting t(18)=4.43, p=.0003224. Then I try to do it with lme:
>
> a12 <- c(a1,a2)
> grp <- factor(rep(1:2,c(10,10)))
>
> Now, at this point, I think I should be able to do something like
> this:
> lme(a12~grp)
> or
> lme(a12~1|grp)
> but I keep getting an error message, "Invalid formula for
> groups." So I tried making a groupedData object:
>
> data1 <- as.data.frame(cbind(a12,grp))
> gd1 <- groupedData(a12~1|grp,data=as.data.frame(cbind(a12,grp)))
>
> Now I can do
> lme(gd1)
> or
> lme(gd1,random=1|grp)
> or many other things, but nothing seems to yield anything like
> the t test, and I'm not even sure what the fixed effect test
> (with a p of .011 with summary(lme(gd1))) is testing. (It
> doesn't seem to be about whether the grand mean of a12 is greater
> than zero.) I've been studying the relevant documentaion,
> including Pinheiro and Bates's book, but I'm still stumped. I'm
> sure I'm being very dense about something very simple, like,
> "This doesn't make any sense." But why not?
>
> All this is leading up to a real application to a much more
> complicated problem, but I think I need to understand the simple
> stuff first.
I think I have to disagree a little with previous correspondents. It
would be useful to have lme fit a model with no random effects, but it
currently will not. You can fool it in two ways to produce the t-test:
indiv <- 1:20
summary(lme(a12~grp,random=~1|indiv))
one <- rep(1,20)
summary(lme(a12~grp,random=~-1|one))
Notice that lme in the first version doesn't notice the aliasing of
the individual and the residual variance:
(Intercept) Residual
StdDev: 2.834877 1.063079
where only the total variance (2.834877^2 + 1.063079^2) is really
identifiable. The fixed effect analysis is similar to lm():
Fixed effects: a12 ~ grp
Value Std.Error DF t-value p-value
(Intercept) 5.5 0.9574271 18 5.744563 <.0001
grp2 6.0 1.3540064 18 4.431294 3e-04
In the other case, you get an essentially arbitrary random component
for the overall level:
(Intercept) Residual
StdDev: 1.805342 3.027650
whereas the residual is correct. The random component causes the
intercept to have an increased variance:
Fixed effects: a12 ~ grp
Value Std.Error DF t-value p-value
(Intercept) 5.5 2.043508 18 2.691450 0.0149
grp2 6.0 1.354006 18 4.431294 0.0003
I think this is wrong. The random level is unidentifiable, so you
should probably get an infinite SE for the intercept.
--
O__ ---- Peter Dalgaard Blegdamsvej 3
c/ /'_ --- Dept. of Biostatistics 2200 Cph. N
(*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk) FAX: (+45) 35327907
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list