R help - Jul 2006 - Test for equality of coefficients in multivariate multipleregression

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
>

G'day John,
>>>>> "JF" == John Fox <jfox at mcmaster.ca>
writes:
    JF> Simply stacking the problems and treating the resulting
    JF> observations as independent will give you the correct
    JF> coefficients, but incorrect coefficient variances
Yes, after Andrew's (off-list) answer I realised this too.  If I am
not mistaken, all variances/covariances should be off by a factor of
1/2 or something like that.

    JF> and artificially zero covariances.
Well, I must admit that I misread Ulrich's code for most of the day.
I hadn't realised that the variable `tmp' introduces a correlation
between `y1' and `y2' in his code:
> 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)
for some reason, my brain kept parsing this as generate *one* random
intercept for y1 and *one* random intercept for y2, not that each
individual observation has a random intercept.  Under the model that
my brain kept parsing, one would have zero covariances. :)

Now I understand why Andrew suggested the use of mixed models and
would go down that way too.  But I believe your approach is valid too.

    JF> (BTW, it's a bit unclear to me how much of this exchange was
    JF> on r-help,
Easy, all those that have r-help either in the TO: or CC: field.
Those were Ulrich's original message and the answer by you and Andrew,
I kept all my mails so far off-list.

    JF> but I'm copying to r-help since at least one of Ulrich's
    JF> messages referring to alternative approaches appeared there.
Yes, I noticed that and answered off-list.  In that message, if I read
it correctly, he had confused Andrew and me.

    JF> I hope that's OK.)
Sure, why not? :)

Cheers,

        Berwin