[R] beginner's questions about lme, fixed and random effects

[Peter Dalgaard BSA]

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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._