Statistics 850 Spring 2005Example of “treatment contrasts” used by R in estimating ANOVA coefficientsThe first example shows a simple numerical design matrix in R (no factors) for the groups “1”, “a”, “b”,“ab”.> resp <- rnorm(4,0,1)> data <- data.frame(resp=resp,+ A=c(0,1,0,1),+ B=c(0,0,1,1))> dataresp A B1 0.3112628 0 02 -0.7630832 1 03 -1.7276357 0 14 -0.5787976 1 1> lm1.out <- lm(resp~A*B,data=data)Note that the columns of the design matrix have zeros where our effect contrasts had “-1”s. The results(fitted values and tests) are equivalent although this design matrix does not have orthogonal columns andthe contrast coefficients do not sum to zero as can be seen from XTX.> X <- model.matrix(lm1.out)> X(Intercept) A B A:B1 1 0 0 02 1 1 0 03 1 0 1 04 1 1 1 1> X%*%t(X)1 2 3 41 1 1 1 12 1 2 1 23 1 1 2 24 1 2 2 4> coef <- solve(t(X)%*%X)%*%t(X)%*%data$resp> coef[,1](Intercept) 0.3112628A -1.0743461B -2.0388985A:B 2.2231841> lm1.out$coef(Intercept) A B A:B0.3112628 -1.0743461 -2.0388985 2.2231841Page 1Statistics 850 Spring 2005An equivalent example with factors rather than numeric design> data <- data.frame(resp=resp,+ A=c("f","m","f","m"),+ B=c("f","f","m","m"))> dataresp A B1 0.3112628 f f2 -0.7630832 m f3 -1.7276357 f m4 -0.5787976 m mR drops the first group. Since f comes before m in the alphabet the contrast for ”f” is dropped.> contrasts(data$A)mf 0m 1> contrasts(data$B)mf 0m 1> lm2.out <- lm(resp~A*B,data=data)> X <- as.num(model.matrix(lm2.out))> X[,1] [,2] [,3] [,4][1,] 1 0 0 0[2,] 1 1 0 0[3,] 1 0 1 0[4,] 1 1 1 1Note that the coefficients calculated using the model matrix and those from the R fit match and match theprevious example.> coef <- solve(t(X)%*%X)%*%t(X)%*%data$resp> coef[,1][1,] 0.3112628[2,] -1.0743461[3,] -2.0388985[4,] 2.2231841> lm2.out$coef(Intercept) Am Bm Am:Bm0.3112628 -1.0743461 -2.0388985 2.2231841> lm1.out$coef(Intercept) A B A:B0.3112628 -1.0743461 -2.0388985 2.2231841Page 2Statistics 850 Spring 2005Example of a 3x2 factorial. A has levels ”low”, ”med” and ”high” and B has levels ”low” and ”high”.> resp <- rnorm(6,0,1)> data <- data.frame(resp=resp,+ A=c("low","med","high",+ "low","med","high"),+ B=c("low","low","low",+ "high","high","high"))I’m specifying the factor levels here so we get them in order.> data$A <- factor(data$A,levels=c("low","med","high"))> data$B <- factor(data$B,levels=c("low","high"))> dataresp A B1 -0.2601737 low low2 0.1730960 med low3 -0.7590816 high low4 0.9682281 low high5 -0.4736121 med high6 -0.3177872 high high>> contrasts(data$A)med highlow 0 0med 1 0high 0 1> contrasts(data$B)highlow 0high 1>> lm.out <- lm(resp~A*B,data=data)>> X <- model.matrix(lm.out)> X(Intercept) Amed Ahigh Bhigh Amed:Bhigh Ahigh:Bhigh1 1 0 0 0 0 02 1 1 0 0 0 03 1 0 1 0 0 04 1 0 0 1 0 05 1 1 0 1 1 06 1 0 1 1 0 1attr(,"contrasts")attr(,"contrasts")$A[1] "contr.treatment"attr(,"contrasts")$B[1] "contr.treatment"Page 3Statistics 850 Spring 2005> X%*%t(X)1 2 3 4 5 61 1 1 1 1 1 12 1 2 1 1 2 13 1 1 2 1 1 24 1 1 1 2 2 25 1 2 1 2 4 26 1 1 2 2 2 4> coef <- solve(t(X)%*%X)%*%t(X)%*%data$resp> coef[,1](Intercept) -0.2601737Amed 0.4332697Ahigh -0.4989079Bhigh 1.2284018Amed:Bhigh -1.8751099Ahigh:Bhigh -0.7871074> lm.out$coef(Intercept) Amed Ahigh Bhigh Amed:Bhigh Ahigh:Bhigh-0.2601737 0.4332697 -0.4989079 1.2284018 -1.8751099 -0.7871074> summary(lm.out)Residuals:ALL 6 residuals are 0: no residual degrees of freedom!Coefficients:Estimate Std. Error t value Pr(>|t|)(Intercept) -0.2602 NA NA NAAmed 0.4333 NA NA NAAhigh -0.4989 NA NA NABhigh 1.2284 NA NA NAAmed:Bhigh -1.8751 NA NA NAAhigh:Bhigh -0.7871 NA NA NAIf A is a blocking variable and if we assume no iteration we can do testing.> lm.out <- lm(resp~A+B,data=data)> summary(lm.out)Coefficients:Estimate Std. Error t value Pr(>|t|)(Intercept) 0.1835 0.5436 0.338 0.768Amed -0.5043 0.6658 -0.757 0.528Ahigh -0.8925 0.6658 -1.340 0.312Bhigh 0.3410 0.5436 0.627 0.595> anova(lm.out)Analysis of Variance TableResponse: respDf Sum Sq Mean Sq F value Pr(>F)A 2 0.80098 0.40049 0.9035 0.5254B 1 0.17442 0.17442 0.3935 0.5945Residuals 2 0.88655 0.44328Page 4Statistics 850 Spring 2005R can create many kinds of c ontrasts. Note that the polynomial and Helmert contrasts are orthogonal.> contr.treatment(2)21 02 1> contr.treatment(5)2 3 4 51 0 0 0 02 1 0 0 03 0 1 0 04 0 0 1 05 0 0 0 1> con <- contr.poly(5)> con.L .Q .C ^4[1,] -0.6324555 0.5345225 -3.162278e-01 0.1195229[2,] -0.3162278 -0.2672612 6.324555e-01 -0.4780914[3,] 0.0000000 -0.5345225 -4.095972e-16 0.7171372[4,] 0.3162278 -0.2672612 -6.324555e-01 -0.4780914[5,] 0.6324555 0.5345225 3.162278e-01 0.1195229> round(t(tt)%*%tt,8).L .Q .C ^4.L 1 0 0 0.Q 0 1 0 0.C 0 0 1 0^4 0 0 0 1>> contr.sum(5)[,1] [,2] [,3] [,4]1 1 0 0 02 0 1 0 03 0 0 1 04 0 0 0 15 -1 -1 -1 -1>> con <- contr.helmert(5)[,1] [,2] [,3] [,4]1 -1 -1 -1 -12 1 -1 -1 -13 0 2 -1 -14 0 0 3 -15 0 0 0 4Page 5Statistics 850 Spring 2005Here is the same example with polynomial contrasts for B.> contrasts(data$A) <- contr.poly(3)> contrasts(data$B) <- contr.treatment(2)> dataresp A B1 -0.2601737 low low2 0.1730960 med low3 -0.7590816 high low4 0.9682281 low high5 -0.4736121 med high6 -0.3177872 high high> lm.out <- lm(resp~A+B,data=data)> round(model.matrix(lm.out),4)(Intercept) A.L A.Q B21 1 -0.7071 0.4082 02 1 0.0000 -0.8165 03 1 0.7071 0.4082 04 1 -0.7071 0.4082 15 1 0.0000 -0.8165 16 1 0.7071 0.4082 1attr(,"contrasts")attr(,"contrasts")$A.L .Qlow -7.071068e-01 0.4082483med -7.850462e-17 -0.8164966high 7.071068e-01 0.4082483attr(,"contrasts")$B2low 0high 1> summary(lm.out)Coefficients:Estimate Std. Error t value Pr(>|t|)(Intercept) -0.2821 0.3844 -0.734 0.539A.L -0.6311 0.4708 -1.340 0.312A.Q 0.0474 0.4708 0.101 0.929B2 0.3410 0.5436 0.627 0.595Residual standard error: 0.6658 on 2 degrees of freedomMultiple R-Squared: 0.5239, Adjusted R-squared: -0.1904F-statistic: 0.7335 on 3 and 2 DF, p-value: 0.6208> anova(lm.out)Analysis of Variance TableResponse: respDf Sum Sq Mean Sq F value Pr(>F)A 2 0.80098 0.40049 0.9035 0.5254B 1 0.17442 0.17442 0.3935 0.5945Residuals 2 0.88655 0.44328Page 6Statistics 850 Spring 2005Here is an example with m ore than one observation per “cell”.> resp <- rnorm(12,0,1)> data <- data.frame(resp=resp,+ A=c( "low","low"+ ,"med","med",+ "high","high",+ "low","low"+ ,"med","med",+ "high","high"),+ B=c(+ "low","low","low",+ "low","low","low",+ "high","high","high",+ "high","high","high"))>> data$A <- factor(data$A,levels=c("low","med","high"))> data$B