# [R] how to estimate treatment-interaction contrasts

szhan at uoguelph.ca szhan at uoguelph.ca
Mon Jul 16 18:22:53 CEST 2007

```Hello, R experts,
Could I use the command suggested by Chunk to get the estimates for
treat-interaction contrasts (non-orthognal) which indicates below:
summary(lm(score ~ A*B))
or is there any other ways to get the estimates?
Again, Chunk, thank you for your help!
Joshua
Quoting Chuck Cleland <ccleland at optonline.net>:

> szhan at uoguelph.ca wrote:
>> Hello, Chuck,
>> Thank you very much for your help! But the contrasts I want to do
>> simutaneously is
>>   contrasts(B)
>>     [,1] [,2] [,3] [,4]
>> b1   -4   -3   -2   -1
>> b2    1   -3   -2   -1
>> b3    1    2   -2   -1
>> b4    1    2    3   -1
>> b5    1    2    3    4
>>
>> Could you please show me how to calculate estimates for ALL
>> intearaction constrasts using THESE contrasts? Say C2: c(-3, -3, 2, 2,
>> 2) as an example. I used the ortholognal constrasts as you suggest,
>> estimate for interaction contrast C2 is still -24.1.
>> Joshua
>
> Joshua:
>   I use a variation on my contrasts to show how you might get the
> estimates.  I then confirm the estimates by calculating the differences
> of differences in means "by hand".
>
> 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/5,  1/2,  0,  0,
>                          3/5, -1/2,  0,  0,
>                         -2/5,  0,  2/3,  0,
>                         -2/5,  0, -1/3,  1/2,
>                         -2/5,  0, -1/3, -1/2), 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
>
> # Estimates with tests that match those above
>
> summary(lm(score ~ A*B))
>
> Call:
> lm(formula = score ~ A * B)
>
> Residuals:
>      Min       1Q   Median       3Q      Max
> -1.16667 -0.29167  0.03333  0.29167  0.73333
>
> Coefficients:
>             Estimate Std. Error t value Pr(>|t|)
> (Intercept)   6.4933     0.1187  54.705  < 2e-16 ***
> Aa2           0.6533     0.1679   3.892 0.000905 ***
> B1            0.2889     0.2423   1.192 0.247086
> B2            0.4000     0.3754   1.066 0.299271
> B3            0.1333     0.3251   0.410 0.686038
> B4           -1.5333     0.3754  -4.085 0.000577 ***
> Aa2:B1       -1.3389     0.3426  -3.907 0.000873 ***
> Aa2:B2        0.3667     0.5308   0.691 0.497665
> Aa2:B3       -1.3833     0.4597  -3.009 0.006932 **
> Aa2:B4        0.3667     0.5308   0.691 0.497665
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> Residual standard error: 0.4597 on 20 degrees of freedom
> Multiple R-Squared: 0.8038,     Adjusted R-squared: 0.7156
> F-statistic: 9.107 on 9 and 20 DF,  p-value: 2.324e-05
>
> # Confirm that estimates match differences of differences in means
>
> diff(rowMeans(tapply(score, list(A,B), mean)[,1:2]) -
>      rowMeans(tapply(score, list(A,B), mean)[,3:5]))
>        a2
> -1.338889
>
> diff(tapply(score, list(A,B), mean)[,1] -
>      tapply(score, list(A,B), mean)[,2])
>        a2
> 0.3666667
>
> diff(tapply(score, list(A,B), mean)[,3] -
>      rowMeans(tapply(score, list(A,B), mean)[,4:5]))
>        a2
> -1.383333
>
> diff(tapply(score, list(A,B), mean)[,4] -
>      tapply(score, list(A,B), mean)[,5])
>        a2
> 0.3666667
>
>   So, applying the same strategy to your contrasts would give the following:
>
> contrasts(B) <- matrix(c(-4/5, -3/5, -2/5, -1/5,
>                           1/5, -3/5, -2/5, -1/5,
>                           1/5,  2/5, -2/5, -1/5,
>                           1/5,  2/5,  3/5, -1/5,
>                           1/5,  2/5,  3/5,  4/5), 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 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
>
> summary(lm(score ~ A*B))
>
> Call:
> lm(formula = score ~ A * B)
>
> Residuals:
>      Min       1Q   Median       3Q      Max
> -1.16667 -0.29167  0.03333  0.29167  0.73333
>
> Coefficients:
>               Estimate Std. Error   t value Pr(>|t|)
> (Intercept)  6.493e+00  1.187e-01    54.705  < 2e-16 ***
> Aa2          6.533e-01  1.679e-01     3.892 0.000905 ***
> B1          -4.000e-01  3.754e-01    -1.066 0.299271
> B2          -3.815e-16  3.754e-01 -1.02e-15 1.000000
> B3          -9.000e-01  3.754e-01    -2.398 0.026373 *
> B4           1.533e+00  3.754e-01     4.085 0.000577 ***
> Aa2:B1      -3.667e-01  5.308e-01    -0.691 0.497665
> Aa2:B2       6.000e-01  5.308e-01     1.130 0.271719
> Aa2:B3       1.567e+00  5.308e-01     2.951 0.007893 **
> Aa2:B4      -3.667e-01  5.308e-01    -0.691 0.497665
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> Residual standard error: 0.4597 on 20 degrees of freedom
> Multiple R-Squared: 0.8038,     Adjusted R-squared: 0.7156
> F-statistic: 9.107 on 9 and 20 DF,  p-value: 2.324e-05
>
>   Because your contrasts are correlated, the order in which they are
> entered makes a difference.  Thus, the two summaries only agree for the
> last interaction contrast (Aa2:B4).
>   Perhaps more intuitively, you also can see that the tests are
> sensitive to order for your contrasts as follows (no output shown):
>
> # Second Contrast First
>
> anova(lm(score ~
>          A * I(B %in% c("b1","b2")) +
>                  A * I(B %in% "b1") +
>     A * I(B %in% c("b1","b2","b3")) +
>  A * I(B %in% c("b1","b2","b3","b4"))), test="F")
>
> # Second Contrast Last
>
> anova(lm(score ~ A * I(B %in% "b1") +
>     A * I(B %in% c("b1","b2","b3")) +
> A * I(B %in% c("b1","b2","b3","b4")) +
>           A * I(B %in% c("b1","b2"))), test="F")
>
>   So if you really do want to evaluate your correlated contrasts
> simultaneously, what makes the most sense to me is:
>
> summary(lm(score ~ A*B))
>
>   I would be interested in other people's thoughts on this, particularly
> how to use something like estimable() in the gmodels package to achieve
> these contrasts.
>
> hope this helps,
>
> Chuck
>
>> Quoting Chuck Cleland <ccleland at optonline.net>:
>>
>>> 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
>>>>
>>>> interaction contrast and corresponding F value?
>>>> Joshua
>>>   If 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
>>>> 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
>>>
>> ______________________________________________
>> R-help at stat.math.ethz.ch mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> 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
>

```