2
$\begingroup$

I have a data as follows.

data <- structure(list(Loc = structure(c(17L, 27L, 33L, 6L, 39L, 10L, 
                                         9L, 11L, 36L, 4L, 2L, 5L, 34L, 45L, 18L, 15L, 29L, 31L, 14L, 
                                         40L, 44L, 19L, 1L, 22L, 3L, 32L, 16L, 21L, 20L, 13L, 12L, 25L, 
                                         35L, 8L, 24L, 41L, 30L, 23L, 7L, 37L, 42L, 26L, 28L, 38L, 43L, 
                                         17L, 27L, 33L, 6L, 39L, 10L, 9L, 11L, 36L, 4L, 2L, 5L, 34L, 45L, 
                                         18L, 15L, 29L, 31L, 14L, 40L, 44L, 19L, 1L, 22L, 3L, 32L, 16L, 
                                         21L, 20L, 13L, 12L, 25L, 35L, 8L, 24L, 41L, 30L, 23L, 7L, 37L, 
                                         42L, 26L, 28L, 38L, 43L), .Label = c("Hpan100786", "Hpan100794", 
                                                                              "Hpan107006", "Hpan109517", "Hpan109543", "Hpan109546", "Hpan109549p", 
                                                                              "Hpan109551", "Hpan109552", "Hpan109553", "Hpan109556", "Hpan112828", 
                                                                              "Hpan113401", "Hpan114525", "Hpan116055", "Hpan14479", "Hpan25167", 
                                                                              "Hpan2573", "Hpan34342", "Hpan36372", "Hpan39020", "Hpan39039", 
                                                                              "Hpan4496", "Hpan6103", "Hpan61434", "Hpan65773", "Hpan76756", 
                                                                              "Hpan93602", "Hpan95269", "Hpan95278", "Lpan15539", "Lpan16561", 
                                                                              "Lpan16823", "Lpan21755", "Lpan24060", "Lpan24994", "Lpan25763", 
                                                                              "Lpan3", "Lpan501186", "Lpan501220", "Lpan501263", "Lpan501286", 
                                                                              "Lpan501468", "Lpan501493", "Lpan501531"), class = "factor"), 
                       Class = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
                                           1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
                                           2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
                                           3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
                                           1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
                                           2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
                                           3L, 3L, 3L, 3L, 3L), .Label = c("Good", "Medium", "Poor"), class = "factor"), 
                       Rep = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
                                         1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
                                         1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
                                         1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
                                         2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
                                         2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
                                         2L, 2L, 2L, 2L, 2L), .Label = c("R1", "R2"), class = "factor"), 
                       y = c(1202, 686, 921, 705, 1060, 642, 996, 1154, 1195, 692, 
                             1045, 1008, 1107, 1195, 1267, 703, 1371, 1122, 992, 1721, 
                             989, 503, 1023, 740, 1169, 716, 625, 1403, 783, 1182, 1411, 
                             851, 1191, 1444, 1643, 1185, 1635, 1525, 1162, 1472, 626, 
                             575, 721, 849, 867, 1178, 703, 845, 689, 1011, 659, 1022, 
                             1203, 1202, 716, 1004, 1004, 1082, 1176, 1248, 719, 1328, 
                             1185, 966, 1685, 1050, 510, 998, 721, 1173, 718, 614, 1392, 
                             767, 1035, 1430, 902, 1225, 1451, 1615, 1157, 1659, 1506, 
                             1151, 1466, 618, 583, 736, 885, 876)), .Names = c("Loc", 
                                                                               "Class", "Rep", "y"), class = "data.frame", row.names = c(NA, 
                                                                                                                                         -90L))

I am trying to do ANOVA with planned contrasts in R using aov. The 45 Loc belong to 3 different Class. I want to compare the Loc in one class with those in another Class. I am trying to follow A (sort of) Complete Guide to Contrasts in R steps for "Running Fewer than J-1 Contrasts for J Groups".

#Creation of contrast matrix
# reorder the levels of `Loc` according to `Class`
data$Loc <- factor(data$Loc, levels = unique(data[order(data$Class), "Loc"]))

c1 <- c(rep(1/15, 15), rep(-1/15, 15), rep(0, 15)) # Good vs Med
c2 <- c(rep(1/15, 15), rep(0, 15), rep(-1/15, 15)) # Good vs Poor
c3 <- c(rep(0, 15), rep(1/15, 15), rep(-1/15, 15)) # Med vs Poor

# Specify the contrasts
mat <- cbind(c1, c2, c3)

options(contrasts = c("contr.helmert", "contr.poly"))
m1 <- aov(y ~ Loc, data = data)
summary(m1)
            Df  Sum Sq Mean Sq F value Pr(>F)    
Loc         44 8547697  194266   292.3 <2e-16 ***
Residuals   45   29910     665                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

m2 <- aov(y ~ Loc, data = data, contrasts=list(Loc = mat))
summary.aov(m2, split=list(Loc=list("Good Vs Med"=1,
                                     "Good Vs Poor" = 2,
                                     "Med Vs Poor"=3)))

                    Df  Sum Sq Mean Sq F value Pr(>F)  
Loc                  2  483311  241656   2.597 0.0802 .
  Loc: Good Vs Med   1    1363    1363   0.015 0.9039  
  Loc: Good Vs Poor  1  481948  481948   5.180 0.0253 *
  Loc: Med Vs Poor   1                                 
Residuals           87 8094296   93038                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The third contrast is missing.

Also the contrast c3 effect is missing

names(m2$coefficients)
[1] "(Intercept)" "Locc1"       "Locc2"       "Locc3"      
names(m2$effects)
 [1] "(Intercept)" "Locc1"       "Locc2"       ""            ""            ""           
 [7] ""            ""            ""            ""            ""            ""           
[13] ""            ""            ""            ""            ""            ""           
[19] ""            ""            ""            ""            ""            ""           
[25] ""            ""            ""            ""            ""            ""           
[31] ""            ""            ""            ""            ""            ""           
[37] ""            ""            ""            ""            ""            ""           
[43] ""            ""            ""           

How to get the appropirate comparisons here? What am I doing wrong in specifying the contrast matrix?

In SPSS I can do the comparisons using the following syntax

UNIANOVA y BY Loc
/METHOD=SSTYPE(3)
/INTERCEPT=INCLUDE
/CRITERIA=ALPHA(0.05)
/LMATRIX = "Good_vs_Med" y -15, 15, -15, 15, 15, 15, 0, 0, 15, 15, 15, 0, -15, -15, -15, -15, 15, 15, -15, -15, -15, -15, 0, 0, 0, 0, 15, 0, -15, 0, -15, -15, 15, 15, 0, 15, 0, 0, 15, -15, 0, 0, 0, -15, 15;
/LMATRIX = "Good_vs_Poor" y 0, 15, 0, 15, 15, 15, -15, -15, 15, 15, 15, -15, 0, 0, 0, 0, 15, 15, 0, 0, 0, 0, -15, -15, -15, -15, 15, -15, 0, -15, 0, 0, 15, 15, -15, 15, -15, -15, 15, 0, -15, -15, -15, 0, 15;
/LMATRIX = "Med_vs_Poor" y 15, 0, 15, 0, 0, 0, -15, -15, 0, 0, 0, -15, 15, 15, 15, 15, 0, 0, 15, 15, 15, 15, -15, -15, -15, -15, 0, -15, 15, -15, 15, 15, 0, 0, -15, 0, -15, -15, 0, 15, -15, -15, -15, 15, 0;
/PRINT = TEST(LMATRIX)
/DESIGN=y
$\endgroup$
9
  • $\begingroup$ At a guess R insists on only using 2 df for comparing 3 groups and SPSS does not. You might get better answers on an R site. $\endgroup$
    – mdewey
    Jun 22 '17 at 17:14
  • 1
    $\begingroup$ You have only 2 degrees of freedom. 3 groups, 1 grand mean, 2 degrees of freedom. Are you testing the hypothesis $\mu_1=\mu_2=\mu_3$ (sane) or $\mu_1=\mu_2=\mu_3=0$ (not clear). $\endgroup$
    – AdamO
    Jun 22 '17 at 17:30
  • 1
    $\begingroup$ @AdamO There are 45 levels in Loc, so 44 df (as is clear in m1). Something is going wrong when aov() estimates m2. Interesting problem! $\endgroup$ Jun 22 '17 at 17:35
  • 1
    $\begingroup$ A related post on SO (unfortunately without a useful answer): https://stackoverflow.com/questions/31818174/custom-contrasts-in-r-contrast-coefficient-matrix-or-contrast-matrix-coding-s $\endgroup$ Jun 22 '17 at 17:36
  • 1
    $\begingroup$ @Crops then the 2 level ANOVA is correct for hypothesis one, and Tukey post hoc for the next hypotheses. 2-4 are redundant. $\endgroup$
    – AdamO
    Jun 22 '17 at 18:27
4
$\begingroup$

Interesting problem! The issue is that you're basically re-coding Loc into a variable with three levels (corresponding to your three levels of Class in this case), which leaves it with only 2 degrees of freedom. R can't estimate three contrasts for a factor with three levels in a model with an intercept.

You can see this because you get 44 df in m1 when R uses its automatic contrasts, but only 2 when you specify your own contrast matrix in m2. The problem is that the contrasts you're specifying create a singular matrix:

> contrasts(data$Loc) <- mat
Error in `contrasts<-`(`*tmp*`, value = c(0.0666666666666667, 0.0666666666666667,  : 
  singular contrast matrix

If you modify your contrast matrix so it's not a complete coding scheme for three levels, R will go back to filling in the missing df with complementary orthogonal contrasts to the ones you specify:

c4 <- c(rep(0, 5), rep(1/15, 15), rep(-1/15, 15), rep(0, 10)) # made up c4
mat <- cbind(c1, c2, c3, c4)
m2 <- aov(y ~ Loc, data = data, contrasts=list(Loc = mat))
summary.aov(m2, split=list(Loc=list("Good Vs Med"=1,
                                    "Good Vs Poor" = 2,
                                    "Med vs Poor"=3)))

Results:

                    Df  Sum Sq Mean Sq F value Pr(>F)  
Loc                  3  683581  227860   2.482 0.0663 .
  Loc: Good Vs Med   1    1363    1363   0.015 0.9033  
  Loc: Good Vs Poor  1  481948  481948   5.250 0.0244 *
  Loc: Med vs Poor   1  200270  200270   2.182 0.1433  
Residuals           86 7894026   91791                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Note that when you specify the contrasts within the model call (aov), it limits the df for the predictor to the number of contrasts you ran - 1, and any additional SS goes into the residual SS. Above, I used a pretty gross hack to get 3 df (so we could see all of your comparisons) by specifying 4 contrasts.

If you specify your contrasts as an attribute of the variable Loc before running the model, R will attempt to auto-fill the rest of the contrast matrix out to the full 44 df, but it still can't in this case because the matrix is singular.

A better approach

It sounds like you're actually trying to get comparisons among the levels of Class, not Loc, but you want all three possible comparisons so you're trying to get Loc to do it for you since Class can't. Instead of running contrasts, you should be looking at pairwise comparisons for Class. You can do that with pairwise.t.test from the stats package, but I prefer the lsmeans package for nicer formatting and more flexibility.

m3 <- aov(y ~ Class, data = data)
library(lsmeans)

You can then use lsmeans functions to compare your different classes on y. The default output gives you the estimated marginal means (which is the same as the raw means, in this case, since you don't have any covariates) with 95%CI around them:

> lsmeans(m3, ~ Class)
 Class     lsmean       SE df  lower.CL upper.CL
 Good    987.2333 55.68898 87  876.5454 1097.921
 Medium  996.7667 55.68898 87  886.0788 1107.455
 Poor   1147.2333 55.68898 87 1036.5454 1257.921

You can see all pairwise comparisons by wrapping that with the pairs function:

> pairs(lsmeans(m3, ~ Class))
 contrast         estimate       SE df t.ratio p.value
 Good - Medium   -9.533333 78.75612 87  -0.121  0.9920
 Good - Poor   -160.000000 78.75612 87  -2.032  0.1107
 Medium - Poor -150.466667 78.75612 87  -1.911  0.1418

P value adjustment: tukey method for comparing a family of 3 estimates 

It defaults to a sensible post hoc adjustment (Tukey's HSD), but you can change that with the adjust argument if you want. If you turn off the adjustment for multiple comparisons, you'll see that one of your comparisons is significant with uncorrected p-values:

> pairs(lsmeans(m3, ~ Class), adjust = "none")
 contrast         estimate       SE df t.ratio p.value
 Good - Medium   -9.533333 78.75612 87  -0.121  0.9039
 Good - Poor   -160.000000 78.75612 87  -2.032  0.0452
 Medium - Poor -150.466667 78.75612 87  -1.911  0.0594

This is a much more sensible way to compare each of the three classes in your data.

$\endgroup$
3
  • $\begingroup$ How do I set c4? Is it arbitrary? If I am using c4 <- c(rep(-1/15, 15), rep(1/30, 30)), again it is a singular matrix. $\endgroup$
    – Crops
    Jun 23 '17 at 7:14
  • $\begingroup$ If I am using each contrast individually, I am getting the sum of squares as in SPSS. mat <- cbind(c3) m2 <- aov(y ~ Loc, data = data, contrasts=list(Loc = mat)) summary.aov(m2, split=list(Loc=list("Good Vs Med"=1, "Good Vs Poor" = 2, "Med vs Poor"=3))) $\endgroup$
    – Crops
    Jun 23 '17 at 7:16
  • $\begingroup$ @Crops as long as you include c1, c2, and c3 as you have them, you'll end up with a singular matrix. Adding an arbitrary c4 will get R to print results for you, but the F-ratios will depend on the residual and so may change depending on the c4 you create --- not ideal. I don't recommend you test the comparisons with contrasts. Instead, use pairwise comparisons. $\endgroup$ Jun 23 '17 at 16:50
2
$\begingroup$

My advice is to forget contrasts in aov and use the lsmeans or multcomp packages instead. They allow you to run as many contrasts as you would like without being limited by degrees of freedom. They also allow you to run a test on contrasts that span more than one line.

The follow page has examples. http://rcompanion.org/rcompanion/h_01.html. (Caveat: I am the author of this page.)

$\endgroup$
0

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