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

[Prof Brian Ripley]

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