[R] Linear model: contrasts

Ian Fiske ianfiske at gmail.com
Fri Feb 6 18:38:16 CET 2009


The problem comes from mixing up general linear model (LM) theory to compute
B with the classical anova estimators.  The two methods use different
approaches to solving the normal equations.  LM theory uses any generalized
inverse of X'X to solve the normal equations.  Yours comes from ginv() which
uses SVD.  On the other hand, classical anova theory uses imposed contraints
to give the system a unique solution.  I'm not sure which generalized
inverse would correspond to the one that imposes the classical anova
constraints, but it must not be the one resulting from SVD in your example. 
Without the imposed anova constraints (called cell means), you won't have
the interactions summing to zero as you want.

hope that helps,
Ian



Stefan Kraemer-2 wrote:
> 
> Hey,
> I am modelling a linear regression Y=X*B+E. To compute the effect of
> “group” the B-values of the regressors/columns that code the interaction
> effects (col. 5-8 and col. 11-14, see below) have to be weighted with
> non-zero elements within the contrast "Group 1" minus "Group 2" (see
> below). My first understanding was that the interaction effects add up to
> zero in each group. Why do they get non-zero contrast weights? Do you
> recommend any literature discussing this topic? 
> 
> Many Thanks,
> Stefan
> 
> The columns of matrix X dummy-codes
> 1.	Data of Group1 (G1)
> 2.	Data of Group2 (G2)
> 3.	Data condition A1 (factor A)
> 4.	Data condition A2 (factor A)
> 5.	Data in cell A1G1
> 6.	Data in cell A2G1 
> 7.	Data in cell A1G2
> 8.	Data in cell A2G2
> 9.	Data condition B1 (factor B)
> 10.	Data condition B2 (factor B)
> 11.	Data in cell B1G1
> 12.	Data in cell B2G1
> 13.	Data in cell B1G2
> 14.	Data in cell B2G2
> 15.	Mean
> 
> Example:
> #create X
> dimnames = list( paste("sub",c(1:8)),
> c("Group1","Group2","A1","A2","A1G1","A2G1", "A1G2","A2G2","B1","B2",
> "B1G1","B2G1", "B1G2","B2G2","mean" ))
> X=matrix(nrow=8,ncol=15,dimnames=dimnames)
> X[,1]= c(rep(1,4), rep(0,4))
> X[,2]=c(rep(0,4),rep(1,4))
> X[,3]=c(1,1,0,0,1,1,0,0);
> X[,4]=c(0,0,1,1,0,0,1,1);
> X[,5]=c(1,1,0,0,0,0,0,0);
> X[,6]=c(0,0,1,1,0,0,0,0);
> X[,7]=c(0,0,0,0,1,1,0,0);
> X[,8]=c(0,0,0,0,0,0,1,1);
> X[,9]=c(1, 0, 1, 0, 1, 0,1, 0);
> X[,10]=c(0, 1, 0, 1, 0, 1, 0, 1);
> X[,11]=c(1, 0, 1, 0,0, 0, 0, 0);
> X[,12]=c( 0, 1, 0, 1,0, 0, 0, 0);
> X[,13]=c( 0, 0, 0, 0,  1, 0, 1, 0)
> X[,14]=c(0, 0, 0, 0,    0, 1, 0, 1)
> X[,15]=c(1,1,1,1,1,1,1,1)
> #y data
> Y=c(1.2,5.1,1.7,7,6.1,8,1.1,2.5)
> 
> library(MASS)
> ginv()
> B=(ginv(t(X)%*%X) %*% t(X)) %*% Y
> 
> #contrast
> c1=rep(0,15)
> c1[1]=1
> c1[c(5,6,11,12)]=1/2
> c2=rep(0,15)
> c2[2]=-1
> c2[c(7,8,13,14)]=-1/2
> 
> #Group effect
> mean(Y[1:4])-mean(Y[5:8])
> c1%*%B+c2%*%B
> 
> -- 
> 
> 
> -- 
> Jetzt 1 Monat kostenlos! GMX FreeDSL - Telefonanschluss + DSL 
> für nur 17,95 Euro/mtl.!* http://dsl.gmx.de/?ac=OM.AD.PD003K11308T4569a
> 
> ______________________________________________
> R-help at r-project.org 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.
> 
> 

-- 
View this message in context: http://www.nabble.com/Linear-model%3A-contrasts-tp21871201p21876891.html
Sent from the R help mailing list archive at Nabble.com.




More information about the R-help mailing list