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.