[R] Test for equality of coefficients in multivariate multipleregression

[John Fox]

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:berwin at bossiaea.maths.uwa.edu.au] 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 <A.Robinson at ms.unimelb.edu.au> 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: berwin at maths.uwa.edu.au
> Australia                        http://www.maths.uwa.edu.au/~berwin
>