John Fox

2006-Jul-19 11:56 UTC

### [R] 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 >

Berwin A Turlach

2006-Jul-19 12:43 UTC

### [R] Test for equality of coefficients in multivariate multipleregression

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