From: Ulrich Keller <uhkeller_at_web.de>

Date: Wed 19 Jul 2006 - 20:09:36 EST

DF2$x1 <- with(DF, c(x1, x1))

DF2$wh <- rep(c(0,1), each=100)

> fixef(mm1)

(Intercept) wh x11 x21 x12 x22 0.07800993 0.15234579 0.52936947 0.13853332 0.37285132 0.46048418 > coef(x.mlm)

> mm2 <- update(mm1, y~wh + x1 + x12 + x22) > anova(mm1, mm2)

R-help@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Wed Jul 19 20:22:51 2006

Date: Wed 19 Jul 2006 - 20:09:36 EST

Hello and thank you for your answers, Andrew and Berwin. If I'm not mistaken, the mixed-model version of Berwin's approach would be:

#My stuff:

DF<-data.frame(x1=rep(c(0,1),each=50),x2=rep(c(0,1),50))
tmp<-rnorm(100)

DF$y1<-tmp+DF$x1*.5+DF$x2*.3+rnorm(100,0,.5) DF$y2<-tmp+DF$x1*.5+DF$x2*.7+rnorm(100,0,.5) x.mlm<-lm(cbind(y1,y2)~x1+x2,data=DF)

#1st part of Andrew's suggestion:

DF2 <- with(DF, data.frame(y=c(y1,y2))) DF2$x11 <- with(DF, c(x1, rep(0,100))) DF2$x21 <- with(DF, c(x2, rep(0,100))) DF2$x12 <- with(DF, c(rep(0,100), x1)) DF2$x22 <- with(DF, c(rep(0,100), x2))

DF2$x1 <- with(DF, c(x1, x1))

DF2$wh <- rep(c(0,1), each=100)

#Mixed version of models:

> DF2$unit <- rep(c(1:100), 2) > library(nlme) > mm1 <- lme(y~wh + x11 + x21 + x12 + x22, random= ~1 | unit, DF2,method="ML")

> fixef(mm1)

(Intercept) wh x11 x21 x12 x22 0.07800993 0.15234579 0.52936947 0.13853332 0.37285132 0.46048418 > coef(x.mlm)

y1 y2 (Intercept) 0.07800993 0.2303557 x1 0.52936947 0.3728513 x2 0.13853332 0.4604842

> mm2 <- update(mm1, y~wh + x1 + x12 + x22) > anova(mm1, mm2)

Model df AIC BIC logLik Test L.Ratio p-value mm1 1 8 523.6284 550.0149 -253.8142 mm2 2 7 522.0173 545.1055 -254.0086 1 vs 2 0.388908 0.5329

This seems to be correct. What anova() tells me is that the effect of x1 is the same for y1 and y2. What I don't understand then is why the coefficients for x12 and x22 differ so much between mm1 and mm2:

> fixef(mm2)

(Intercept) wh x1 x12 x22
0.1472766 0.1384474 0.5293695 -0.1565182 0.3497476
> fixef(mm1)

(Intercept) wh x11 x21 x12 x22
0.07800993 0.15234579 0.52936947 0.13853332 0.37285132 0.46048418

Sorry for being a bit slow here, I'm (obviously) not a statistician. Thanks again,

Uli

Andrew Robinson wrote:

*> G'day Berwin,
**>
**> my major reason for preferring the mixed-effect approach is that, as
**> you can see below, the residual df for your models are 195 and 194,
**> respectively. The 100 units are each contributing two degrees of
**> freedom to your inference, and presumably to your estimate of the
**> variance. I feel a little sqeueamish about that, because it's not
**> clear to me that they can be assumed to be independent. I worry about
**> the effects on the size of the test. With a mixed-effects model each
**> unit could be a cluster of two observations, and I would guess the
**> size would be closer to nominal, if not nominal.
**>
**> Cheers
**>
**> Andrew
*

R-help@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Wed Jul 19 20:22:51 2006

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.1.8, at Wed 19 Jul 2006 - 22:16:57 EST.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*