[R] Test for equality of coefficients in multivariate multipleregression

John Fox jfox at mcmaster.ca
Wed Jul 19 13:56:13 CEST 2006


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
>



More information about the R-help mailing list