11
$\begingroup$

I am doing a one way ANOVA (per species) with custom contrasts.

     [,1] [,2] [,3] [,4]
0.5    -1    0    0    0
5       1   -1    0    0
12.5    0    1   -1    0
25      0    0    1   -1
50      0    0    0    1

where I compare intensity 0.5 against 5, 5 against 12.5 and so on. These is the data I'm working on

enter image description here

with the following results

Generalized least squares fit by REML
  Model: dark ~ intensity 
  Data: skofijski.diurnal[skofijski.diurnal$species == "niphargus", ] 
       AIC      BIC    logLik
  63.41333 67.66163 -25.70667

Coefficients:
            Value Std.Error  t-value p-value
(Intercept) 16.95 0.2140872 79.17334  0.0000
intensity1   2.20 0.4281744  5.13809  0.0001
intensity2   1.40 0.5244044  2.66970  0.0175
intensity3   2.10 0.5244044  4.00454  0.0011
intensity4   1.80 0.4281744  4.20389  0.0008

 Correlation: 
           (Intr) intns1 intns2 intns3
intensity1 0.000                      
intensity2 0.000  0.612               
intensity3 0.000  0.408  0.667        
intensity4 0.000  0.250  0.408  0.612 

Standardized residuals:
       Min         Q1        Med         Q3        Max 
-2.3500484 -0.7833495  0.2611165  0.7833495  1.3055824 

Residual standard error: 0.9574271 
Degrees of freedom: 20 total; 15 residual

16.95 is the global mean for "niphargus". In intensity1, I'm comparing means for intensity 0.5 against 5.

If I understood this right, the coefficient for intensity1 of 2.2 should be half the difference between means of intensity levels 0.5 and 5. However, my hand calculations don't match those of the summary. Can anyone chip in what am I doing wrong?

ce1 <- skofijski.diurnal$intensity
levels(ce1) <- c("0.5", "5", "0", "0", "0")
ce1 <- as.factor(as.character(ce1))
tapply(skofijski.diurnal$dark, ce1, mean)
       0    0.5      5 
  14.500 11.875 13.000 
diff(tapply(skofijski.diurnal$dark, ce1, mean))/2
      0.5       5 
  -1.3125  0.5625
Improve this question
$\endgroup$
5
  • $\begingroup$ Could you provide the lm() function from R that you used to estimate. How exactly did you use the contrasts function? $\endgroup$
    – Philippe
    Commented Sep 11, 2012 at 0:50
  • $\begingroup$ btw geom_points(position=position_dodge(width=0.75)) will fix the way the points in your plot don't line up with the boxes. $\endgroup$
    – flies
    Commented Sep 10, 2018 at 13:49
  • $\begingroup$ @flies since my question, there's been an introduction of geom_jitter, which is a shortcut for all the geom_point() parameters that jitter. $\endgroup$
    – Roman Luštrik
    Commented Sep 11, 2018 at 16:06
  • $\begingroup$ I didn't notice the jitter there. does geom_jitter(position_dodge) work? I've been using geom_points(position_jitterdodge) to add dots to boxplots with dodging. $\endgroup$
    – flies
    Commented Sep 12, 2018 at 21:21
  • $\begingroup$ @flies see the docs for geom_jitter here. In my experience since my above answer, I find it unnecessary to use boxplots. Ever. If I have many points, I use violin plots which show density of points in much finer details than boxplots. Boxplots were invented back when plotting many points or their densities was not convenient. Perhaps it's time we start thinking of dropping this (handicapped) visualization. $\endgroup$
    – Roman Luštrik
    Commented Sep 14, 2018 at 7:19

2 Answers 2

Reset to default
13
$\begingroup$

The matrix you specified for the contrasts is correct in principle. To convert it into an appropriate contrast matrix, you need to calculate the generalized inverse of your original matrix.

If M is your matrix:

M

#     [,1] [,2] [,3] [,4]
#0.5    -1    0    0    0
#5       1   -1    0    0
#12.5    0    1   -1    0
#25      0    0    1   -1
#50      0    0    0    1

Now, calculate the generalized inverse using ginv and transpose the result using t:

library(MASS)
t(ginv(M))

#     [,1] [,2] [,3] [,4]
#[1,] -0.8 -0.6 -0.4 -0.2
#[2,]  0.2 -0.6 -0.4 -0.2
#[3,]  0.2  0.4 -0.4 -0.2
#[4,]  0.2  0.4  0.6 -0.2
#[5,]  0.2  0.4  0.6  0.8

The result is identical to the one of @Greg Snow. Use this matrix for your analysis.

This is a much easier way than doing it manually.

There is an even easier way to generate a matrix of sliding differences (a.k.a. repeated contrasts). This can be done with the function contr.sdif and the number of factor levels as a parameter. If you have five factor levels, like in your example:

library(MASS)
contr.sdif(5)

#   2-1  3-2  4-3  5-4
#1 -0.8 -0.6 -0.4 -0.2
#2  0.2 -0.6 -0.4 -0.2
#3  0.2  0.4 -0.4 -0.2
#4  0.2  0.4  0.6 -0.2
#5  0.2  0.4  0.6  0.8
Improve this answer
$\endgroup$
6
$\begingroup$

If the matrix at the top is how you are encoding the dummy variables (what you are passing to the C or contrast function in R) then they the first one is comparing the 1st level to the others (actually 0.8 times the 1st subtracted from 0.2 times the sum of the others).

The second term compares the 1st 2 levels to the last 3. The 3rd compares the 1st 3 levels to the last2 and the 4th compares the 1st 4 levels to the last one.

If you want to do the comparisons that you describe (compare each pair) then the dummy variable encoding that you want is:

      [,1] [,2] [,3] [,4]
[1,] -0.8 -0.6 -0.4 -0.2
[2,]  0.2 -0.6 -0.4 -0.2
[3,]  0.2  0.4 -0.4 -0.2
[4,]  0.2  0.4  0.6 -0.2
[5,]  0.2  0.4  0.6  0.8
Improve this answer
$\endgroup$
2
  • $\begingroup$ Is the use of this generalized inverese matrix also necessary when using aov() instead of lm()? I'm asking, because I have read several tutorials, in which contrast matrices for aov() are constructed just as the one given by Roman. E.g. see Chapter 5 in cran.r-project.org/doc/contrib/Vikneswaran-ED_companion.pdf $\endgroup$
    – crsh
    Commented Oct 18, 2012 at 16:14
  • 2
    $\begingroup$ The aov function calls the lm function to do the main computations, so things like contrast matrices will have the same effect in both. $\endgroup$
    – Greg Snow
    Commented Oct 18, 2012 at 16:26

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.