Statistics 850 Spring 2005Example of estimating coefficients and fitted values in a 2x2 anova with one observationper cell using different main effect contrasts.>> resp <- rnorm(4,0,1)> resp[1] 0.2064763 0.1232091 -1.1272979 1.5425342>> ## First parameterization> data <- data.frame(resp=resp,+ A=c(1,1,-1,-1),+ B=c(1,-1,1,-1))> dataresp A B1 0.2064763 1 12 0.1232091 1 -13 -1.1272979 -1 14 1.5425342 -1 -1>> lm.out <- lm(resp~A*B,data=data)> ## note that this model is equivalent to resp ~ A + B + A:B>> fit1 <- lm.out$fit>> X1 <- model.matrix(lm.out)>> X1(Intercept) A B A:B1 1 1 1 12 1 1 -1 -13 1 -1 1 -14 1 -1 -1 1attr(,"assign")[1] 0 1 2 3>> t(X1)%*%X1(Intercept) A B A:B(Intercept) 4 0 0 0A 0 4 0 0B 0 0 4 0A:B 0 0 0 4>> solve(t(X1)%*%X1)(Intercept) A B A:B(Intercept) 0.25 0.00 0.00 0.00A 0.00 0.25 0.00 0.00B 0.00 0.00 0.25 0.00A:B 0.00 0.00 0.00 0.25>Page 1Statistics 850 Spring 2005> solve(t(X1)%*%X1) %*% t(X1)1 2 3 4(Intercept) 0.25 0.25 0.25 0.25A 0.25 0.25 -0.25 -0.25B 0.25 -0.25 0.25 -0.25A:B 0.25 -0.25 -0.25 0.25>> coef1 <- solve(t(X1)%*%X1) %*% t(X1)%*%resp> coef1[,1](Intercept) 0.18623043A -0.02138770B -0.64664121A:B 0.68827483>> coef(lm.out)(Intercept) A B A:B0.18623043 -0.02138770 -0.64664121 0.68827483>> cbind(fit1, X1%*%coef1)fit11 0.2064763 0.20647632 0.1232091 0.12320913 -1.1272979 -1.12729794 1.5425342 1.5425342>Second parameterization> data <- data.frame(resp=resp,+ A=1/2*c(1,1,-1,-1),+ B=1/2*c(1,-1,1,-1))> dataresp A B1 0.2064763 0.5 0.52 0.1232091 0.5 -0.53 -1.1272979 -0.5 0.54 1.5425342 -0.5 -0.5>> lm.out <- lm(resp~A*B,data=data)>> X2 <- model.matrix(lm.out)>> X2(Intercept) A B A:B1 1 0.5 0.5 0.252 1 0.5 -0.5 -0.253 1 -0.5 0.5 -0.254 1 -0.5 -0.5 0.25attr(,"assign")[1] 0 1 2 3>Page 2Statistics 850 Spring 2005> t(X2)%*%X2(Intercept) A B A:B(Intercept) 4 0 0 0.00A 0 1 0 0.00B 0 0 1 0.00A:B 0 0 0 0.25>> solve(t(X2)%*%X2)%*%t(X2)1 2 3 4(Intercept) 0.25 0.25 0.25 0.25A 0.50 0.50 -0.50 -0.50B 0.50 -0.50 0.50 -0.50A:B 1.00 -1.00 -1.00 1.00>> coef2 <- solve(t(X2)%*%X2)%*%t(X2)%*%resp> coef2[,1](Intercept) 0.18623043A -0.04277541B -1.29328242A:B 2.75309932>> fit2 <- X2%*%coef2>> coef1/coef2[,1](Intercept) 1.00A 0.50B 0.50A:B 0.25>> X1/X2(Intercept) A B A:B1 1 2 2 42 1 2 2 43 1 2 2 44 1 2 2 4attr(,"assign")[1] 0 1 2 3>> as.matrix(fit1/fit2)[,1]1 12 13 14 1>Page
View Full Document