QK Box 8.12

Partridge and Farquhar (1981) studied the effect of number of mating partners on longevity of fruitflies. There were five treatments: one virgin female per day, one newly inseminated (pregnant) female per day, eight virgin females per day, eight newly inseminated (pregnant) females per day, and a control group with no females. There were 25 males, individually kept in separate vials, in each group. The thorax length of each individual fly was also recorded as a covariate. If thorax length explains some of the variation in longevity, then the evaluation of the effect of treatment on longevity adjusted for thorax length will be more precise. The raw data were extracted by reading from Figure 2 in the original paper (see also description and discussion in Hanley & Shapiro 1994). Our general H₀ was that there was no effect of partner treatment on longevity of male fruitflies, adjusting for thorax length.

This example was also used in the first edition; the data file is here.

Partridge, L. & Farquhar, M. (1981). Sexual activity and the lifespan of male fruitflies. Nature, 294, 580-81.

Preliminaries

First, load the required packages (car, effects)

Import partridge data file (partridge.csv)

partridge <- read.csv("../data/partridge.csv")
head(partridge,10)

partridge$treatment <- factor(partridge$treatment)

Check assumptions

boxplot(longev ~ treatment, data = partridge)

boxplot(thorax ~ treatment, data = partridge)

Boxplots don’t indicate any skewness or unequal variances

scatterplot(longev ~ thorax | treatment, data = partridge)

Lines look linear and close to parallel

Check that covariate not affected by treatments

cov.aov <- aov(thorax ~ treatment, data = partridge)
summary(cov.aov)

             Df Sum Sq Mean Sq F value Pr(>F)
treatment     4  0.030 0.00750    1.26   0.29
Residuals   120  0.714 0.00595               

Fit full model with interaction to evaluate homogeneous slopes.

All SS types produce same result for interaction

partridge.aov1 <- aov(longev ~ treatment + thorax + treatment*thorax, data = partridge)
plot(partridge.aov1)

Residual plot shows some pattern with smaller variances at each end - only one unusual value (66). longev values only cover small range so no transformation used

anova(partridge.aov1)

Analysis of Variance Table

Response: longev
                  Df Sum Sq Mean Sq F value  Pr(>F)    
treatment          4  11939    2985   26.20 1.9e-15 ***
thorax             1  13169   13169  115.59 < 2e-16 ***
treatment:thorax   4     43      11    0.09    0.98    
Residuals        115  13102     114                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Conclude no evidence for non-parallelism (P = 0.984)

Refit model leaving out interaction

Set contrasts to get Type III SS

partridge.aov2 <- aov(longev ~ treatment + thorax, contrasts = list(treatment = contr.sum), data = partridge)
anova(partridge.aov2)

Analysis of Variance Table

Response: longev
           Df Sum Sq Mean Sq F value  Pr(>F)    
treatment   4  11939    2985      27 5.9e-16 ***
thorax      1  13169   13169     119 < 2e-16 ***
Residuals 119  13145     110                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Default anova will be Type I SS, so order of factor and covariate in model matter; use Type III SS

Anova(partridge.aov2, type = 'III')

Anova Table (Type III tests)

Response: longev
            Sum Sq  Df F value  Pr(>F)    
(Intercept)   3070   1    27.8 6.1e-07 ***
treatment     9611   4    21.8 1.7e-13 ***
thorax       13169   1   119.2 < 2e-16 ***
Residuals    13145 119                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Compare to type II

Anova(partridge.aov2, type = 'II')

Anova Table (Type II tests)

Response: longev
          Sum Sq  Df F value  Pr(>F)    
treatment   9611   4    21.8 1.7e-13 ***
thorax     13169   1   119.2 < 2e-16 ***
Residuals  13145 119                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Compare MS Residual to a model without covariate

Only single predictor so SS types irrelevant

partridge.aov3 <- aov(longev ~ treatment, data = partridge)
anova(partridge.aov3)

Analysis of Variance Table

Response: longev
           Df Sum Sq Mean Sq F value  Pr(>F)    
treatment   4  11939    2985    13.6 3.5e-09 ***
Residuals 120  26314     219                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Get adjusted means

adjmeans <- effect("treatment", partridge.aov2, se = TRUE)
summary(adjmeans)

 treatment effect
treatment
 none preg1 preg8 virg1 virg8 
 61.5  64.2  65.4  54.5  41.6 

 Lower 95 Percent Confidence Limits
treatment
 none preg1 preg8 virg1 virg8 
 57.3  60.0  61.3  50.3  37.4 

 Upper 95 Percent Confidence Limits
treatment
 none preg1 preg8 virg1 virg8 
 65.7  68.3  69.6  58.7  45.8 

Now for contrasts

Need to centre covariate to get correct intercept (i.e. overall mean or reference category) for contrasts

partridge$cthorax <- scale(partridge$thorax, center = TRUE, scale = FALSE)

Get default contrasts by setting control group (no females) as reference

partridge$treatment <- relevel(partridge$treatment, ref = "none")

Refit model with centered covariate

partridge.aov4 <- aov(longev ~ treatment + cthorax, data = partridge)

get summary.lm output

summary.lm(partridge.aov4)

Call:
aov(formula = longev ~ treatment + cthorax, data = partridge)

Residuals:
    Min      1Q  Median      3Q     Max 
-26.189  -6.599  -0.989   6.408  30.244 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)       61.52       2.11   29.15   <2e-16 ***
treatmentpreg1     2.65       2.98    0.89     0.37    
treatmentpreg8     3.93       3.00    1.31     0.19    
treatmentvirg1    -7.02       2.97   -2.36     0.02 *  
treatmentvirg8   -19.95       3.01   -6.64    1e-09 ***
cthorax          135.82      12.44   10.92   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 10.5 on 119 degrees of freedom
Multiple R-squared:  0.656, Adjusted R-squared:  0.642 
F-statistic: 45.5 on 5 and 119 DF,  p-value: <2e-16

confint(partridge.aov4)

                2.5 % 97.5 %
(Intercept)     57.34  65.70
treatmentpreg1  -3.24   8.54
treatmentpreg8  -2.00   9.86
treatmentvirg1 -12.90  -1.13
treatmentvirg8 -25.90 -14.00
cthorax        111.19 160.45

Get specific planned contrasts

First, preg vs virg and then 1 virg vs 8 virg)

contrasts(partridge$treatment) <- cbind(c(0,0.5,0.5,-0.5,-0.5), c(0,0,0,1,-1))
contrasts(partridge$treatment)

      [,1] [,2]    [,3]   [,4]
none   0.0    0 -0.3855 -0.807
preg1  0.5    0  0.7344 -0.103
preg8  0.5    0 -0.5417  0.507
virg1 -0.5    1  0.0964  0.202
virg8 -0.5   -1  0.0964  0.202

partridge.aov5 <- aov(longev ~ treatment + thorax, data = partridge)
summary.lm(partridge.aov5)

Call:
aov(formula = longev ~ treatment + thorax, data = partridge)

Residuals:
    Min      1Q  Median      3Q     Max 
-26.189  -6.599  -0.989   6.408  30.244 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -54.06      10.26   -5.27  6.1e-07 ***
treatment1     16.77       2.10    7.98  1.0e-12 ***
treatment2      6.47       1.50    4.30  3.5e-05 ***
treatment3     -2.78       2.10   -1.32    0.189    
treatment4     -3.72       2.12   -1.76    0.081 .  
thorax        135.82      12.44   10.92  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 10.5 on 119 degrees of freedom
Multiple R-squared:  0.656, Adjusted R-squared:  0.642 
F-statistic: 45.5 on 5 and 119 DF,  p-value: <2e-16

confint(partridge.aov5)

             2.5 % 97.5 %
(Intercept) -74.37 -33.76
treatment1   12.61  20.94
treatment2    3.49   9.45
treatment3   -6.94   1.39
treatment4   -7.92   0.47
thorax      111.19 160.45

Do same contrasts but partition anova table and use F tests

pregvsvirg = c(0,0.5,0.5,-0.5,-0.5)
virg1vsvirg8 = c(0,0,0,1,-1)
matriz = cbind(pregvsvirg, virg1vsvirg8)
contrasts(partridge$treatment) = matriz
CList=list("pregvsvirg" = 1, "virg1vsvirg8" = 2)
partridge.aov5 <- aov(longev ~ treatment + thorax, data = partridge)
summary(partridge.aov5, split = list(treatment = CList))

                           Df Sum Sq Mean Sq F value  Pr(>F)    
treatment                   4  11939    2985    27.0 5.9e-16 ***
  treatment: pregvsvirg     1   6675    6675    60.4 3.0e-12 ***
  treatment: virg1vsvirg8   1   4068    4068    36.8 1.6e-08 ***
thorax                      1  13169   13169   119.2 < 2e-16 ***
Residuals                 119  13145     110                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Finally, compare to anova on longevity without covariate

partridge.aov6 <- aov(longev ~ treatment, data = partridge)
summary(partridge.aov6)

             Df Sum Sq Mean Sq F value  Pr(>F)    
treatment     4  11939    2985    13.6 3.5e-09 ***
Residuals   120  26314     219                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1