R help - Dec 2001 - beginner's questions about lme, fixed and random effects

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.

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

On Mon, 3 Dec 2001, Jonathan Baron wrote:
> 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?
Right, "This doesn't make any sense.".  To have a random effect
you need a
set of randomly selected groups.  You don't have one.

You can do a paired t test this way, as then the pairs are a randomly
selected group (or could be in principle).

y <- rnorm(20)
pairs <- factor(rep(1:10, 2), labels=LETTERS[1:10])
treat <- factor(rep(c("Y", "N"), rep(10, 2)))
t.test(y ~ treat, paired=T)
data:  y by treat
t = 0.9827, df = 9, p-value = 0.3514

sample estimates:
mean of the differences
              0.3775932


You do need to specify random to lme:
> lme(y ~ treat, random = ~1 | pairs)
Linear mixed-effects model fit by REML
  Data: NULL
  Log-restricted-likelihood: -25.97613
  Fixed: y ~ treat
(Intercept)      treatY
  0.2533409  -0.3775932

Random effects:
 Formula: ~1 | pairs
        (Intercept)  Residual
StdDev:   0.2794975 0.8592174

Number of Observations: 20
Number of Groups: 10
> summary(.Last.value)
Linear mixed-effects model fit by REML
 Data: NULL
       AIC      BIC    logLik
  59.95227 63.51375 -25.97613

Random effects:
 Formula: ~1 | pairs
        (Intercept)  Residual
StdDev:   0.2794975 0.8592174

Fixed effects: y ~ treat
                 Value Std.Error DF    t-value p-value
(Intercept)  0.2533409 0.2857225  9  0.8866678  0.3983
treatY      -0.3775932 0.3842537  9 -0.9826665  0.3514
....


-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272860 (secr)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595

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

Jonathan,

I don't think it makes a lof ot sense to try to use lme for a two-sample
t-test, since there is no random effect there. Maybe a comparison more in line
with what you were trying to do is to compare the results from lme with those
from a split-plot anova. Here is a silly example (where the value of the.data
is the sum of the treatment effects, the subject effects, and
"an.error"):


## Silly example
subject.effects <- rep(seq(1:5),rep(2,5))
subjects.ids <-
as.factor(rep(c("a","b","c","d","e"),rep(2,5)))
treatment.effects <- rep(c(0,4), 5)
treatment.ids <- as.factor(rep(c("trt1", "trt2"), 5))
an.error <- c(1,2,0,1,3,2,0,1,2,0)
the.data <- subject.effects + treatment.effects + an.error

summary(aov(the.data ~ Error(subjects.ids) + treatment.ids))
summary(lme(the.data ~ treatment.ids, random = ~1|subjects.ids))

## Of course, with this setup, those tests are equivalent to a paired-t:
t.test(the.data[treatment.ids == "trt1"], the.data[treatment.ids ==
"trt2"],
paired=TRUE)


Ramón

Jonathan Baron <baron at cattell.psych.upenn.edu> dijo:
> 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.
> 
> Jon Baron
>
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
> 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
>
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
> 


-- 



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

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