[R] beginner's questions about lme, fixed and random effects
Prof Brian Ripley
ripley at stats.ox.ac.uk
Mon Dec 3 13:17:12 CET 2001
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list