5
$\begingroup$

I have conjured up a dataset here:

    df = data.frame(score = c( 6.03, 5.89, 5.75, 6.48, 6.50, 
           6.68, 6.27, 5.61, 6.82, 7.26, 7.83, 7.82, 6.63, 8.70, 
           8.85, 8.81, 8.60, 7.81, 8.95, 9.55,10.36,10.26, 
           9.50,11.00,11.71),
                    treatment = c(rep("treat1", 5), rep("treat2", 
            5), rep("treat3", 5), rep("treat4", 5), 
            rep("treat5", 5)))

Question 1: Is the below contrast matrix theoretically identical to asking for contrasts: "1,1,1,1,-4" which in English means: compare each treatment against "treatment5"?

    options(contrasts=c('contr.sum','contr.poly'))
    contrasts(df$day)

        [,1] [,2] [,3] [,4]
    fri    1    0    0    0
    mon    0    1    0    0
    thu    0    0    1    0
    tue    0    0    0    1
    wed   -1   -1   -1   -1

Question 2: Furthermore, if I wish to only compare "treatment1" and "treatment5", have I set the contrasts correctly in R? (Working shown below)

    contrastmatrix <- matrix(c(1, 0, 0, 0, -1, 0, 0, 0, 0,  0, 0, 
          0, 0, 0, 0, 0, 0, 0, 0, 0), nrow = 5, ncol = 4, 
         byrow = FALSE)
    contrasts(df$treatment) <- contrastmatrix 

    contrasts(df$treatment)

           [,1] [,2] [,3] [,4]
    treat1    1    0    0    0
    treat2    0    0    0    0
    treat3    0    0    0    0
    treat4    0    0    0    0
    treat5   -1    0    0    0

Question 3: Now I wish to test for a linear trend amongst the treatments, hence the contrast settings I wish to set it to is "-2,-1,0,1,2". How do I set this in R? (I.e. should my working be identical to what I have done in Question 2?)

$\endgroup$
3
  • 1
    $\begingroup$ Welcome to Cross Validated! Although questions about R per se are more on-topic at StackOverflow, I voted to keep your question open since setting up contrasts to perform specific comparisons is pretty statistical in nature. $\endgroup$ Jan 12, 2016 at 23:32
  • 1
    $\begingroup$ I agree with @MattKrause. I think this has sufficient statistical content to stay open. $\endgroup$ Jan 13, 2016 at 0:35
  • $\begingroup$ Thanks guys. I was contemplating on putting this on StackOverflow instead, but I think the root of my lack of understanding still lies in my inability to grasp the concept of contrasts hence why I posted my question here instead. For example, I don't really understand why or how contrasts such as "1,1,1,1,-4" can be broken down into the matrix shown above in R nor do I understand how the linear trend contrast (for 5 factors) "-2,-1,0,1,2" even makes sense in R. Any help would be greatly appreciated! $\endgroup$
    – j.cube
    Jan 13, 2016 at 1:00

2 Answers 2

4
$\begingroup$

Response to Question 1:

No. The contrast matrix you have illustrated asks for 4 contrasts: Fri vs Wed, Mon vs Wed, Thu vs Wed, and Tue vs Wed. Wed was chosen because R goes alphabetical. You could control this by ordering the df$day factor.

The contrast you have specified: "1,1,1,1,-4" requests the average of Fri-Mon-Thu-Tue versus Wed and could be specified as cbind(c(1,1,1,1,-4)).

Response to Question 2:

Yes. The contrast in column 1 is comparing Treatment 1 to Treatment 5.

Response to Question 3:

cbind(c(-2,-1,0,1,2),c(0,0,0,0,0),c(0,0,0,0,0),c(0,0,0,0,0)) would generate a matrix in which the first column has the linear trend contrast IF the coefficients are in order. Note that in your day of the week example, your order is Fri-Mon-Thu-Tue-Wed and so this would not be a time-linear-trend since the factor levels are not order. So beware of that and use ordered factors to keep from pulling your hair out.

$\endgroup$
1
  • $\begingroup$ Thank you very very much! I just realised that I accidentally wrote the "mon - fri" factors instead of "treat1 - treat5" for Q1. Also, I get that in a 'normal' one-way ANOVA, the entire matrix must be used for the multiple comparisons required, but in the case for Question 2, it feels redundant to have so many 0's. Question: does R have an easier way to input contrasts such as (1,1,1,1,-4), (1,0,0,0,-1) and (-2,-1,0,1,2) individually? Question: what does it statistically mean if I combine the 3 (just mentioned) contrasts into one matrix and conduct an ANOVA on that? $\endgroup$
    – j.cube
    Jan 13, 2016 at 23:12
6
$\begingroup$

Testing contrasts of factorial variables is notoriously difficult in R. Most of the things that go on beneath the surface in friendlier programs like SPSS must be spelled out specifically in R. This is both a bad and a good thing. Bad because it takes you a long time to learn how, but good because your understanding is vastly improved as a result.

I don't know where the Mon-Fri factor comes from. But, using the treat 1-5 levels from the original dataset supplied, I want to take issue with something @StatNoodle wrote in their (mostly excellent) answer. To do this we need to make the treatment variable a factor and specify the contr.sum set of contrasts.

df$treatment <- factor(df$treatment)
contrasts(df$treatment) <- contr.sum(5)
contrasts(df$treatment)

       [,1] [,2] [,3] [,4]
treat1    1    0    0    0
treat2    0    1    0    0
treat3    0    0    1    0
treat4    0    0    0    1
treat5   -1   -1   -1   -1

I want to clarify that this set of contrasts, does not compare treat1 to treat5, treat2 to treat5, treat3 to treat5 and treat4 to treat5. It is deviation coding (see here) and compares the mean for each of treatments 1 to 4 to the grand mean, that is the mean of all five groups.

This is easy to show if we just get the means for each level

(sumDF <- aggregate(score ~ treatment, data = df, mean))

  treatment  score
1    treat1  6.130
2    treat2  6.528
3    treat3  7.966
4    treat4  8.744
5    treat5 10.566

And calculate the grand mean by taking the average of the group means

(grandMean <- mean(sumDF$score))

[1] 7.9868

Now run an lm() to get the coefficients to compare to the descriptive means

contLM <- lm(score ~ treatment, df)
contLM

Coefficients:
(Intercept)   treatment1   treatment2   treatment3   treatment4  
     7.9868      -1.8568      -1.4588      -0.0208       0.7572

The intercept is our grand mean, 7.9868. The treatment1 coefficient is the average of scores for the treat1 group minus the grand mean 6.130-7.9868 = -1.8568. The treatment2 coefficient is the average of scores for the treat2 group minus the grand mean 6.528-7.9868 = -1.4588. The treatment3 coefficient is the average of scores for the treat3 group minus the grand mean 7.966-7.9868 = -0.0208. The treatment4 coefficient is the average of scores for the treat4 group minus the grand mean 8.744-7.9868 = 0.7572. There is no coefficient comparing the treat5 group to the grand mean (it is taken out of the matrix by the row of -1s).

If you want to compare the average of each group to a reference group, use the default treatment coding...

(contLM <- lm(score ~ treatment, df, contrasts = list(treatment = contr.treatment(5))))

Coefficients:
(Intercept)   treatment2   treatment3   treatment4   treatment5  
      6.130        0.398        1.836        2.614        4.436

...which if you do the arithmetic on the descriptive means above will show that intercept is the average of scores in treat1, the treatment2 coefficient the difference between the mean score in treat2 and the mean score in treat1, the treatment3 coefficient the difference between the mean score of treat3 and the mean score of treat1 etc.

Answers to Questions

If you want to perform custom (i.e. non-default) contrasts in R I would strongly recommend using a package designed specifically for such things. The lsmeans package is excellent.

library(lsmeans)

First get a reference grid using the lsmeans() function

(lsm <- lsmeans(contLM, ~treatment))

 treatment lsmean        SE df lower.CL  upper.CL
 treat1     6.130 0.3086882 20 5.486088  6.773912
 treat2     6.528 0.3086882 20 5.884088  7.171912
 treat3     7.966 0.3086882 20 7.322088  8.609912
 treat4     8.744 0.3086882 20 8.100088  9.387912
 treat5    10.566 0.3086882 20 9.922088 11.209912

The values in the lsmean column of the output are covariate adjusted means for each group. The descending order of these coeffcients - treat1 first to treat5 last` - is important because we use these relative positions to specify what contrasts we wish to test.

Then specify the contrast you want in a list. We want to test the first four treatments (positions 1 to 4) to treat5 (fifth position).

contr <- list("treats1-4 vs treat5" = c(-1, -1, -1, -1, 4))

We then pass this list into the lsmeans::contrast() function

contrast(lsm, contr)

 contrast            estimate       SE df t.ratio p.value
 treats1-4 vs treat5   12.896 1.380496 20   9.342  <.0001

This will yield the correct t ratio for the contrasts (square this value if you want the F ratio) and the correct p value. But notice that the estimate is not what we would expect if we perform the procedure manually on the raw group means, subtracting the mean for treat5 from the average mean scores across groups 1 to 4.

mean(c(6.130, 6.528, 7.966, 8.744)) - 10.566
[1] -3.224

In fact it the coefficient is out by a factor of exactly 4

4*-3.224
[1] -12.896

To get accurate estimates we need to transform the contrast coefficients so each of the two sides of the contrast sum to |1| and the two sides together sum to 0.

contr2 <- list("treats1-4 vs treat5" = c(-.25, -.25, -.25, -.25, 1))
contrast(lsm, contr2)

contrast               estimate        SE df t.ratio p.value
Question1: trt1-4 vs 5    3.224 0.3451239 20   9.342  <.0001

A great feature of lsmeans is that you can pass multiple contrasts into the list

contr3 <- list("Question1: trt1-4 vs 5" = c(-.25, -.25, -.25, -.25, 1),
               "Question2: trt1 vs trt5" = c(1, 0, 0, 0, -1),
               "Question3: linear trend" = c(-0.2, -0.1, 0, 0.1, 0.2))
contrast(lsm, contr3)

 contrast                estimate         SE df t.ratio p.value
 Question1: trt1-4 vs 5    3.2240 0.34512389 20   9.342  <.0001
 Question2: trt1 vs trt5  -4.4360 0.43655103 20 -10.161  <.0001
 Question3: linear trend   1.1088 0.09761578 20  11.359  <.0001

These custom contrasts should answer all of your questions, with accurate estimates. Note the transformed coefficients for the polynomial linear trend contrasts. These sorts of transformations require all kinds of weird voodoo which mostly go on under the hood in SPSS (and which I delved into briefly a few years ago, see here) but are necessary to get accurate estimates. In this case we divide the trend coefficients by their squared sum, then multiply by the incremental change across coefficients (they increase by 1 each time), thus

c(-2, -1, 0, 1, 2)/(sum(c(-2, -1, 0, 1, 2)^2))*1
[1] -0.2 -0.1  0.0  0.1  0.2

We can verify that this estimate of rate of change across groups, 1.109, is accurate because if we transform the treatment factor to a numerical variable

df$treatNum <- rep(1:5, each = 5)

And then run the regression

lm(score ~ treatNum, df)

Coefficients:
(Intercept)     treatNum  
      4.660        1.109  

The coefficient is exactly the same as the estimate returned by the polynomial trend contrast in the contr list, and also matches the rate of increase across factors we can see if we graph score against the numerical treatment variable

ggplot(df, aes(treatNum, score)) + 
       geom_line(stat = "summary", fun.y = "mean") + 
       geom_smooth(method = "lm")

UPDATE: the lsmeans package has been superseded by the emmeans package but they were created by the same person and do more or less the same thing.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.