library(ggplot2) # for general plotting
library(car)     # for ANOVA (Type II used, better than Type I when there is an unbalanced design)

Information of data source

Enzyme inhibitors were used to treat mites, fecundity and enzyme activity following enzyme inhibition were assessed.

Inhibitors:

PBO - CYP (cytochrome P450) inhibitor

DEM - GST inhibitor

DEF - esterase inhibitor

Read in the data and view structure to identify any issues in data formatting

detox.inhibitor.data <- read.csv("~/Lab Stuff/Adapted mites/Tomato/Inhibitor assay/3 dpi feeding data/Inhibitor Adapted vs. non-adapted - Moneymaker/All trials/R data.csv", header = TRUE)

# Trial as a factor
detox.inhibitor.data$Trial <- factor(detox.inhibitor.data$Trial)
str(detox.inhibitor.data)
## 'data.frame':    144 obs. of  4 variables:
##  $ Trial      : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Mite.Strain: Factor w/ 2 levels "TU","TU-A": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Treatment  : Factor w/ 4 levels "DEF","DEM","PBO",..: 4 1 2 3 4 1 2 3 4 1 ...
##  $ Eggs.Mite  : num  6 2.13 3.12 4.3 8.74 ...
# Subset data for diffferent mite strains
detox.inhibitor.data.TU <- subset(detox.inhibitor.data, Mite.Strain =="TU")
detox.inhibitor.data.TA <- subset(detox.inhibitor.data, Mite.Strain =="TU-A")

Formulate hypothesis

H0: There will be no difference in fecundity upon treatment with inhibitors of any class in TU or TU-A mites.

HA: Treatment with detoxification inhibitors will result in degreased fecundity in TU-A mites.

Conduct data exploration

Outliers in the response variable (Eggs.Mite) within explanatory variables (Trial, Treatment, Mite.Strain).

ggplot(detox.inhibitor.data, aes(x = Trial, y = Eggs.Mite)) + geom_boxplot()+ theme_classic()

ggplot(detox.inhibitor.data, aes(x = Treatment, y = Eggs.Mite)) + geom_boxplot()+ theme_classic()

ggplot(detox.inhibitor.data, aes(x = Mite.Strain, y = Eggs.Mite)) + geom_boxplot()+ theme_classic()

Outier left in, probably represents real variability.

Collinearity of the explanatory variables

Does not apply, all explanatory variables are categorical/factorial.

Spatial/temporal or other hierarchical aspects of sampling design

No, I am treating Trial as a main effect to check for reproducibility (not a random effect/blocking factor).

Interactions (is the quality of the data good enough to include them?)

Interaction betweenTrial and Mite.Strain will be performed to test for reproducibility.

Interaction betweenTrial and Treatment will be performed to test for reproducibility.

Interaction betweenMite.Strain and Treatment will be performed to check if the mite strains responded differently to treatment.

Zero inflation in Y

No

Are categorical covariates balanced?

summary(detox.inhibitor.data)
##  Trial  Mite.Strain Treatment    Eggs.Mite     
##  1:48   TU  :72     DEF  :36   Min.   : 1.800  
##  2:48   TU-A:72     DEM  :36   1st Qu.: 4.683  
##  3:48               PBO  :36   Median : 7.528  
##                     Water:36   Mean   : 9.329  
##                                3rd Qu.:12.974  
##                                Max.   :24.200

Yes

Apply models

# fit linear model and display model fit information and ANOVA table
# full model including 3 way interaction term - to verify it is not significant, if it is, interpretation of hypothesis testing will be problematic
m.0 <- lm(Eggs.Mite ~ Mite.Strain + Treatment + Trial + Mite.Strain:Treatment + Treatment:Trial + Mite.Strain:Trial + Mite.Strain:Treatment:Trial, data = detox.inhibitor.data)
summary(m.0)
## 
## Call:
## lm(formula = Eggs.Mite ~ Mite.Strain + Treatment + Trial + Mite.Strain:Treatment + 
##     Treatment:Trial + Mite.Strain:Trial + Mite.Strain:Treatment:Trial, 
##     data = detox.inhibitor.data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.4176 -1.1786 -0.0242  1.1510  8.9606 
## 
## Coefficients:
##                                       Estimate Std. Error t value Pr(>|t|)
## (Intercept)                             3.7545     0.9730   3.859 0.000185
## Mite.StrainTU-A                         6.6683     1.3760   4.846  3.8e-06
## TreatmentDEM                            1.1572     1.3760   0.841 0.402008
## TreatmentPBO                            1.5353     1.3760   1.116 0.266734
## TreatmentWater                          1.5166     1.3760   1.102 0.272565
## Trial2                                  2.1104     1.3760   1.534 0.127723
## Trial3                                 -0.9497     1.3760  -0.690 0.491400
## Mite.StrainTU-A:TreatmentDEM           -0.6162     1.9459  -0.317 0.752055
## Mite.StrainTU-A:TreatmentPBO           -1.7607     1.9459  -0.905 0.367365
## Mite.StrainTU-A:TreatmentWater          0.7605     1.9459   0.391 0.696637
## TreatmentDEM:Trial2                    -2.3323     1.9459  -1.199 0.233060
## TreatmentPBO:Trial2                    -1.3607     1.9459  -0.699 0.485748
## TreatmentWater:Trial2                  -2.5451     1.9459  -1.308 0.193392
## TreatmentDEM:Trial3                     0.1771     1.9459   0.091 0.927653
## TreatmentPBO:Trial3                     0.5909     1.9459   0.304 0.761918
## TreatmentWater:Trial3                  -0.6214     1.9459  -0.319 0.750009
## Mite.StrainTU-A:Trial2                  4.2201     1.9459   2.169 0.032078
## Mite.StrainTU-A:Trial3                  1.9593     1.9459   1.007 0.316024
## Mite.StrainTU-A:TreatmentDEM:Trial2     2.5428     2.7519   0.924 0.357345
## Mite.StrainTU-A:TreatmentPBO:Trial2     3.3032     2.7519   1.200 0.232380
## Mite.StrainTU-A:TreatmentWater:Trial2   4.3728     2.7519   1.589 0.114696
## Mite.StrainTU-A:TreatmentDEM:Trial3    -0.4908     2.7519  -0.178 0.858747
## Mite.StrainTU-A:TreatmentPBO:Trial3    -0.5550     2.7519  -0.202 0.840521
## Mite.StrainTU-A:TreatmentWater:Trial3   2.3653     2.7519   0.860 0.391778
##                                          
## (Intercept)                           ***
## Mite.StrainTU-A                       ***
## TreatmentDEM                             
## TreatmentPBO                             
## TreatmentWater                           
## Trial2                                   
## Trial3                                   
## Mite.StrainTU-A:TreatmentDEM             
## Mite.StrainTU-A:TreatmentPBO             
## Mite.StrainTU-A:TreatmentWater           
## TreatmentDEM:Trial2                      
## TreatmentPBO:Trial2                      
## TreatmentWater:Trial2                    
## TreatmentDEM:Trial3                      
## TreatmentPBO:Trial3                      
## TreatmentWater:Trial3                    
## Mite.StrainTU-A:Trial2                *  
## Mite.StrainTU-A:Trial3                   
## Mite.StrainTU-A:TreatmentDEM:Trial2      
## Mite.StrainTU-A:TreatmentPBO:Trial2      
## Mite.StrainTU-A:TreatmentWater:Trial2    
## Mite.StrainTU-A:TreatmentDEM:Trial3      
## Mite.StrainTU-A:TreatmentPBO:Trial3      
## Mite.StrainTU-A:TreatmentWater:Trial3    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.383 on 120 degrees of freedom
## Multiple R-squared:  0.8557, Adjusted R-squared:  0.8281 
## F-statistic: 30.95 on 23 and 120 DF,  p-value: < 2.2e-16
Anova(m.0)
## Anova Table (Type II tests)
## 
## Response: Eggs.Mite
##                              Sum Sq  Df  F value    Pr(>F)    
## Mite.Strain                 3103.90   1 546.4799 < 2.2e-16 ***
## Treatment                     75.74   3   4.4448  0.005341 ** 
## Trial                        468.85   2  41.2735 2.287e-14 ***
## Mite.Strain:Treatment         76.63   3   4.4970  0.005001 ** 
## Treatment:Trial                8.69   6   0.2549  0.956459    
## Mite.Strain:Trial            285.03   2  25.0918 7.867e-10 ***
## Mite.Strain:Treatment:Trial   24.42   6   0.7166  0.636972    
## Residuals                    681.58 120                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Calculate effect size and display
result.anova<-Anova(m.0)
ss<-result.anova$"Sum Sq"    ##ss = sum of squares
pes<-ss/(ss+ss[length(ss)])  ##pes = partial e squared
pes[length(pes)]<-""
result.anova$"Part E Sq"<-pes
result.anova
## Anova Table (Type II tests)
## 
## Response: Eggs.Mite
##                              Sum Sq  Df  F value  Pr(>F) Part E Sq
## Mite.Strain                 3103.90   1 546.4799 0.00000   0.81995
## Treatment                     75.74   3   4.4448 0.00534   0.10001
## Trial                        468.85   2  41.2735 0.00000   0.40755
## Mite.Strain:Treatment         76.63   3   4.4970 0.00500   0.10106
## Treatment:Trial                8.69   6   0.2549 0.95646   0.01258
## Mite.Strain:Trial            285.03   2  25.0918 0.00000   0.29488
## Mite.Strain:Treatment:Trial   24.42   6   0.7166 0.63697   0.03459
## Residuals                    681.58 120
# plot interactions
interaction.plot(detox.inhibitor.data$Mite.Strain, detox.inhibitor.data$Treatment, detox.inhibitor.data$Eggs.Mite, type="l", leg.bty="o", leg.bg="grey95", lwd=2, ylab="Eggs.Mite", xlab="Mite.Strain", main="Mite.Strain:Treatment")

interaction.plot(detox.inhibitor.data$Treatment, detox.inhibitor.data$Trial, detox.inhibitor.data$Eggs.Mite, type="l", leg.bty="o", leg.bg="grey95", lwd=2, ylab="Eggs.Mite", xlab="Treatment", main="Treatment:Trial")

interaction.plot(detox.inhibitor.data$Mite.Strain, detox.inhibitor.data$Trial, detox.inhibitor.data$Eggs.Mite, type="l", leg.bty="o", leg.bg="grey95", lwd=2, ylab="Eggs.Mite", xlab="Mite.Strain", main="Mite.Strain:Trial")

# linear model without non-significant 3-way interaction and performing analysis within mite strains (want to reduce comparisons made in Tukey-Kramer post-hoc test)
m.TU <- lm(Eggs.Mite ~ Treatment + Trial + Treatment:Trial, data = detox.inhibitor.data.TU)
summary(m.TU)
## 
## Call:
## lm(formula = Eggs.Mite ~ Treatment + Trial + Treatment:Trial, 
##     data = detox.inhibitor.data.TU)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.6582 -0.8296  0.0258  0.7705  3.4658 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             3.7545     0.5641   6.655 9.66e-09 ***
## TreatmentDEM            1.1572     0.7978   1.450   0.1521    
## TreatmentPBO            1.5353     0.7978   1.924   0.0591 .  
## TreatmentWater          1.5166     0.7978   1.901   0.0621 .  
## Trial2                  2.1104     0.7978   2.645   0.0104 *  
## Trial3                 -0.9497     0.7978  -1.190   0.2386    
## TreatmentDEM:Trial2    -2.3323     1.1283  -2.067   0.0430 *  
## TreatmentPBO:Trial2    -1.3607     1.1283  -1.206   0.2326    
## TreatmentWater:Trial2  -2.5451     1.1283  -2.256   0.0277 *  
## TreatmentDEM:Trial3     0.1771     1.1283   0.157   0.8758    
## TreatmentPBO:Trial3     0.5909     1.1283   0.524   0.6024    
## TreatmentWater:Trial3  -0.6214     1.1283  -0.551   0.5838    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.382 on 60 degrees of freedom
## Multiple R-squared:  0.3378, Adjusted R-squared:  0.2164 
## F-statistic: 2.782 on 11 and 60 DF,  p-value: 0.005483
Anova(m.TU)
## Anova Table (Type II tests)
## 
## Response: Eggs.Mite
##                  Sum Sq Df F value   Pr(>F)   
## Treatment        15.366  3  2.6824 0.054720 . 
## Trial            26.242  2  6.8714 0.002056 **
## Treatment:Trial  16.834  6  1.4693 0.204142   
## Residuals       114.570 60                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Calculate effect size and display
result.anova<-Anova(m.TU)
ss<-result.anova$"Sum Sq"  ##ss = sum of squares
pes<-ss/(ss+ss[length(ss)])  ##pes = partial e squared
pes[length(pes)]<-""
result.anova$"Part E Sq"<-pes
result.anova
## Anova Table (Type II tests)
## 
## Response: Eggs.Mite
##                  Sum Sq Df F value   Pr(>F) Part E Sq
## Treatment        15.366  3  2.6824 0.054720   0.11826
## Trial            26.242  2  6.8714 0.002056   0.18636
## Treatment:Trial  16.834  6  1.4693 0.204142   0.12811
## Residuals       114.570 60
# perform post-hoc Tukey-Kramer test of contrasts
m.TU.Tukey <- aov(m.TU)
TukeyHSD(m.TU.Tukey, "Treatment")
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = m.TU)
## 
## $Treatment
##                  diff         lwr       upr     p adj
## DEM-DEF    0.43880350 -0.77838374 1.6559907 0.7765601
## PBO-DEF    1.27871053  0.06152328 2.4958978 0.0358903
## Water-DEF  0.46110764 -0.75607960 1.6782949 0.7493478
## PBO-DEM    0.83990703 -0.37728022 2.0570943 0.2725895
## Water-DEM  0.02230414 -1.19488311 1.2394914 0.9999588
## Water-PBO -0.81760289 -2.03479013 0.3995844 0.2953182
# plot interactions
interaction.plot(detox.inhibitor.data.TU$Treatment, detox.inhibitor.data.TU$Trial, detox.inhibitor.data.TU$Eggs.Mite, type="l", leg.bty="o", leg.bg="grey95", lwd=2, ylab="Eggs.Mite", xlab="Treatment", main="Treatment:Trial")

# linear model without non-significant 3-way interaction and performing analysis within mite strains (want to reduce comparisons made in Tukey-Kramer post-hoc test)
m.TA <- lm(Eggs.Mite ~ Treatment + Trial + Treatment:Trial, data = detox.inhibitor.data.TA)
Anova(m.TA)
## Anova Table (Type II tests)
## 
## Response: Eggs.Mite
##                 Sum Sq Df F value   Pr(>F)    
## Treatment       137.00  3  4.8323 0.004445 ** 
## Trial           727.64  2 38.4992 1.75e-11 ***
## Treatment:Trial  16.27  6  0.2870 0.940917    
## Residuals       567.01 60                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Calculate effect size and display
result.anova<-Anova(m.TA)
ss<-result.anova$"Sum Sq"  ##ss = sum of squares
pes<-ss/(ss+ss[length(ss)])  ##pes = partial e squared
pes[length(pes)]<-""
result.anova$"Part E Sq"<-pes
result.anova
## Anova Table (Type II tests)
## 
## Response: Eggs.Mite
##                 Sum Sq Df F value  Pr(>F) Part E Sq
## Treatment       137.00  3  4.8323 0.00445   0.19460
## Trial           727.64  2 38.4992 0.00000   0.56204
## Treatment:Trial  16.27  6  0.2870 0.94092   0.02790
## Residuals       567.01 60
# perform post-hoc Tukey-Kramer test of contrasts
m.TA.Tukey <- aov(m.TA)
TukeyHSD(m.TA.Tukey, "Treatment")
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = m.TA)
## 
## $Treatment
##                  diff        lwr      upr     p adj
## DEM-DEF    0.50660188 -2.2011923 3.214396 0.9600430
## PBO-DEF    0.43405907 -2.2737351 3.141853 0.9742297
## Water-DEF  3.46757969  0.7597855 6.175374 0.0067511
## PBO-DEM   -0.07254281 -2.7803370 2.635251 0.9998712
## Water-DEM  2.96097781  0.2531837 5.668772 0.0267892
## Water-PBO  3.03352062  0.3257265 5.741315 0.0222160
# plot interactions
interaction.plot(detox.inhibitor.data.TA$Treatment, detox.inhibitor.data.TA$Trial, detox.inhibitor.data.TA$Eggs.Mite, type="l", leg.bty="o", leg.bg="grey95", lwd=2, ylab="Eggs.Mite", xlab="Treatment", main="Treatment:Trial")

Validate models

# TU
detox.inhibitor.data.TU$m.TU.fit <- fitted(m.TU)    # fitted values
detox.inhibitor.data.TU$m.TU.res <- rstandard(m.TU) # Pearson residuals

# TA
detox.inhibitor.data.TA$m.TA.fit <- fitted(m.TA)    # fitted values
detox.inhibitor.data.TA$m.TA.res <- rstandard(m.TA) # Pearson residuals

Residual distribution / Overdispersion

We assumed normal residuals. This is the least important regression assumption but its can be tested with a qq plot.

# TU
ggplot(detox.inhibitor.data.TU, aes(sample = m.TU.res)) + geom_qq() +
  geom_abline(intercept = 0, slope = 1) + theme_classic()

# TA
ggplot(detox.inhibitor.data.TA, aes(sample = m.TA.res)) + geom_qq() +
  geom_abline(intercept = 0, slope = 1) + theme_classic()

Nice.

Residuals vs fitted values

Testing for:

Linearity - there should be no curvilinear pattern in the residuals.

Equal variance - the vertical spread of the residuals should be constant across all fitted values.

# TU
ggplot(detox.inhibitor.data.TU, aes(x = m.TU.fit, y = m.TU.res)) + 
  geom_point() +   geom_hline(yintercept = 0) +  geom_smooth() + theme_classic()
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'