szhan at uoguelph.ca
2007-Jul-12 20:16 UTC
[R] how to estimate treatment-interaction contrasts
Hello, R experts,
Sorry for asking this question again again since I really want a help!
I have a two-factor experiment data and like to calculate estimates of
interation contrasts say factor A has levels of a1, a2, and B has
levels of b1, b2, b3, b4, and b5 with 3 replicates. I am not sure the
constrast estimate I got is right using the script below:
score<-c(7.2,6.5,6.9,6.4,6.9,6.1,6.9,5.3,7.2,5.7,5.1,5.9,7.6,6.9,6.8,
7.2,6.6,6.9,6.4,6.0,6.0,6.9,6.9,6.4,7.5,7.7,7.0,8.6,8.8,8.3)
A <- gl(2, 15, labels=c("a1", "a2"))
B <- rep(gl(5, 3, labels=c("b1", "b2", "b3",
"b4", "b5")), 2)
contrasts(B)<-cbind(c(-4,rep(1,4)),c(rep(-3,2),rep(2,3)),
+ c(rep(-2,3),rep(3,2)),c(rep(-1,4), rep(4,1)))
fit1 <- aov(score ~ A*B)
summary(fit1, split=list(B=1:4), expand.split = TRUE)
Df Sum Sq Mean Sq F value Pr(>F)
A 1 3.2013 3.2013 15.1483 0.0009054 ***
B 4 8.7780 2.1945 10.3841 0.0001019 ***
B: C1 1 0.0301 0.0301 0.1424 0.7099296
B: C2 1 2.0335 2.0335 9.6221 0.0056199 **
B: C3 1 1.2469 1.2469 5.9004 0.0246876 *
B: C4 1 5.4675 5.4675 25.8715 5.637e-05 ***
A:B 4 5.3420 1.3355 6.3194 0.0018616 **
A:B: C1 1 0.7207 0.7207 3.4105 0.0796342 .
A:B: C2 1 2.6068 2.6068 12.3350 0.0021927 **
A:B: C3 1 1.9136 1.9136 9.0549 0.0069317 **
A:B: C4 1 0.1008 0.1008 0.4771 0.4976647
Residuals 20 4.2267 0.2113
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05
'.' 0.1 ' ' 1
Now I like to get interaction contrast estimate for b1 and b2 vs b3, b4 and b5
cont <- c(1, -1)[A] * c(-3, -3, 2, 2, 2)[B]
estimat<-sum(cont*score) # value of the contrast estimate for A:B C2
> estimat
[1] -24.1
I am not sure the estimate for A:B C2 contrast (-24.1) is correct
because the F value given the output above(12.3350) does not equal to
those I calculate below (15.2684):
t.stat <- sum(cont*score)/se.contrast(fit1,
as.matrix(cont))> t.stat^2
Contrast 1
15.2684
Could you please help me calculate the correct the estimate of
interaction contrast and corresponding F value?
Thanks in advance!
Joshua
szhan at uoguelph.ca wrote:> Hello, R experts, > Sorry for asking this question again again since I really want a help! > > I have a two-factor experiment data and like to calculate estimates of > interation contrasts say factor A has levels of a1, a2, and B has > levels of b1, b2, b3, b4, and b5 with 3 replicates. I am not sure the > constrast estimate I got is right using the script below: > > score<-c(7.2,6.5,6.9,6.4,6.9,6.1,6.9,5.3,7.2,5.7,5.1,5.9,7.6,6.9,6.8, > 7.2,6.6,6.9,6.4,6.0,6.0,6.9,6.9,6.4,7.5,7.7,7.0,8.6,8.8,8.3) > > A <- gl(2, 15, labels=c("a1", "a2")) > B <- rep(gl(5, 3, labels=c("b1", "b2", "b3", "b4", "b5")), 2) > > contrasts(B)<-cbind(c(-4,rep(1,4)),c(rep(-3,2),rep(2,3)), > + c(rep(-2,3),rep(3,2)),c(rep(-1,4), rep(4,1))) > fit1 <- aov(score ~ A*B) > summary(fit1, split=list(B=1:4), expand.split = TRUE) > Df Sum Sq Mean Sq F value Pr(>F) > A 1 3.2013 3.2013 15.1483 0.0009054 *** > B 4 8.7780 2.1945 10.3841 0.0001019 *** > B: C1 1 0.0301 0.0301 0.1424 0.7099296 > B: C2 1 2.0335 2.0335 9.6221 0.0056199 ** > B: C3 1 1.2469 1.2469 5.9004 0.0246876 * > B: C4 1 5.4675 5.4675 25.8715 5.637e-05 *** > A:B 4 5.3420 1.3355 6.3194 0.0018616 ** > A:B: C1 1 0.7207 0.7207 3.4105 0.0796342 . > A:B: C2 1 2.6068 2.6068 12.3350 0.0021927 ** > A:B: C3 1 1.9136 1.9136 9.0549 0.0069317 ** > A:B: C4 1 0.1008 0.1008 0.4771 0.4976647 > Residuals 20 4.2267 0.2113 > --- > Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > > Now I like to get interaction contrast estimate for b1 and b2 vs b3, b4 and b5 > cont <- c(1, -1)[A] * c(-3, -3, 2, 2, 2)[B] > > estimat<-sum(cont*score) # value of the contrast estimate for A:B C2 > >> estimat > [1] -24.1 > > I am not sure the estimate for A:B C2 contrast (-24.1) is correct > because the F value given the output above(12.3350) does not equal to > those I calculate below (15.2684): > > t.stat <- sum(cont*score)/se.contrast(fit1, as.matrix(cont)) >> t.stat^2 > Contrast 1 > 15.2684 > > Could you please help me calculate the correct the estimate of > interaction contrast and corresponding F value? > Thanks in advance! > JoshuaIf the contrasts for B are orthogonal, then you get the result you expected: score <- c(7.2,6.5,6.9,6.4,6.9,6.1,6.9,5.3,7.2,5.7,5.1,5.9,7.6,6.9,6.8, 7.2,6.6,6.9,6.4,6.0,6.0,6.9,6.9,6.4,7.5,7.7,7.0,8.6,8.8,8.3) A <- gl(2, 15, labels=c("a1", "a2")) B <- rep(gl(5, 3, labels=c("b1", "b2", "b3", "b4", "b5")), 2) contrasts(B) <- matrix(c(3, -1, 0, 0, 3, 1, 0, 0, -2, 0, 2, 0, -2, 0, -1, 1, -2, 0, -1, -1), ncol=4, byrow=TRUE) fit1 <- aov(score ~ A*B) summary(fit1, split=list(B=1:4), expand.split = TRUE) Df Sum Sq Mean Sq F value Pr(>F) A 1 3.2013 3.2013 15.1483 0.0009054 *** B 4 8.7780 2.1945 10.3841 0.0001019 *** B: C1 1 1.0427 1.0427 4.9340 0.0380408 * B: C2 1 1.0208 1.0208 4.8304 0.0399049 * B: C3 1 1.2469 1.2469 5.9004 0.0246876 * B: C4 1 5.4675 5.4675 25.8715 5.637e-05 *** A:B 4 5.3420 1.3355 6.3194 0.0018616 ** A:B: C1 1 3.2267 3.2267 15.2684 0.0008734 *** A:B: C2 1 0.1008 0.1008 0.4771 0.4976647 A:B: C3 1 1.9136 1.9136 9.0549 0.0069317 ** A:B: C4 1 0.1008 0.1008 0.4771 0.4976647 Residuals 20 4.2267 0.2113 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Note that I put the contrast of interest for B in the first column of the contrast matrix. hope this helps, Chuck> ______________________________________________ > R-help at stat.math.ethz.ch 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.-- Chuck Cleland, Ph.D. NDRI, Inc. 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894