Dear Berwin,
Simply stacking the problems and treating the resulting observations as
independent will give you the correct coefficients, but incorrect
coefficient variances and artificially zero covariances.
The approach that I suggested originally -- testing a linear hypothesis
using the coefficient estimates and covariances from the multivariate linear
model -- seems simple enough. For example, to test that all three
coefficients are the same across the two equations,
b <- as.vector(coef(x.mlm))
V <- vcov(x.mlm)
L <- c(1, 0, 0,-1, 0, 0,
0, 1, 0, 0,-1, 0,
0, 0, 1, 0, 0,-1)
L <- matrix(L, nrow=3, byrow=TRUE)
t(L %*% b) %*% (L %*% V %*% t(L)) %*% (L %*% b)
The test statistic is chi-square with 3 df under the null hypothesis. (Note:
not checked carefully.)
(BTW, it's a bit unclear to me how much of this exchange was on r-help, but
I'm copying to r-help since at least one of Ulrich's messages referring to
alternative approaches appeared there. I hope that's OK.)
Regards,
John
--------------------------------
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox
--------------------------------
> -----Original Message-----
> From: Berwin A Turlach
> [mailto:[EMAIL PROTECTED] On Behalf Of Berwin
> A Turlach
> Sent: Tuesday, July 18, 2006 9:28 PM
> To: Andrew Robinson
> Cc: Ulrich Keller; John Fox
> Subject: Re: [R] Test for equality of coefficients in
> multivariate multipleregression
>
> G'day all,
>
> >>>>> "AR" == Andrew Robinson <[EMAIL PROTECTED]> writes:
>
> AR> I suggest that you try to rewrite the model system into a
> AR> single mixed-effects model, [...] Why a mixed-effect
> model, wouldn't a fixed effect be o.k. too?
>
> Something like:
>
> > 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)
> > coef(x.mlm)
> y1 y2
> (Intercept) -0.08885266 -0.05749196
> x1 0.33749086 0.60395258
> x2 0.72017894 1.11932077
>
>
> > 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)
> > fm1 <- lm(y~wh + x11 + x21 + x12 + x22, DF2)
> > fm1
>
> Call:
> lm(formula = y ~ wh + x11 + x21 + x12 + x22, data = DF2)
>
> Coefficients:
> (Intercept) wh x11 x21
> x12 x22
> -0.08885 0.03136 0.33749 0.72018
> 0.60395 1.11932
>
> > fm2 <- lm(y~wh + x1 + x21 + x22, DF2)
> > anova(fm2,fm1)
> Analysis of Variance Table
>
> Model 1: y ~ wh + x1 + x21 + x22
> Model 2: y ~ wh + x11 + x21 + x12 + x22
> Res.Df RSS Df Sum of Sq F Pr(>F)
> 1 195 246.919
> 2 194 246.031 1 0.888 0.6998 0.4039
>
>
> Cheers,
>
> Berwin
>
> ========================== Full address ============================
> Berwin A Turlach Tel.: +61 (8) 6488 3338
> (secr)
> School of Mathematics and Statistics +61 (8) 6488 3383
> (self)
> The University of Western Australia FAX : +61 (8) 6488 1028
> 35 Stirling Highway
> Crawley WA 6009 e-mail: [EMAIL PROTECTED]
> Australia http://www.maths.uwa.edu.au/~berwin
>
______________________________________________
[email protected] 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.
