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

Hormone were analysed in Moneymaker tomato leaflets that were fed on by no mites, TU mites and TU-A mites for 24 hours.

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

JA.data <- read.csv("~/Lab Stuff/Adapted mites/Tomato/Hormones/R data/JA R data.csv", header = TRUE)
# trial as a factor
JA.data$Trial <- factor(JA.data$Trial)
str(JA.data)
## 'data.frame':    13 obs. of  3 variables:
##  $ Trial    : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 2 2 2 ...
##  $ Treatment: Factor w/ 3 levels "control","TU",..: 1 3 3 3 2 2 2 1 1 3 ...
##  $ ng.g     : num  13.5 24.7 130.2 103.6 130.9 ...
JAI.data <- read.csv("~/Lab Stuff/Adapted mites/Tomato/Hormones/R data/JA-Ile R data.csv", header = TRUE)
# trial as a factor
JAI.data$Trial <- factor(JAI.data$Trial)
str(JAI.data)
## 'data.frame':    14 obs. of  3 variables:
##  $ Trial    : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 2 2 ...
##  $ Treatment: Factor w/ 3 levels "control","TU",..: 1 1 3 3 3 2 2 2 1 1 ...
##  $ ng.g     : num  7.67 7.1 6.19 21.61 13.02 ...
SA.data <- read.csv("~/Lab Stuff/Adapted mites/Tomato/Hormones/R data/SA R data.csv", header = TRUE)
# trial as a factor
SA.data$Trial <- factor(SA.data$Trial)
str(SA.data)
## 'data.frame':    15 obs. of  3 variables:
##  $ Trial    : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 2 2 ...
##  $ Treatment: Factor w/ 3 levels "control","TU",..: 1 1 3 3 3 2 2 2 1 1 ...
##  $ ng.g     : num  1568 516 573 2043 2028 ...
ABA.data <- read.csv("~/Lab Stuff/Adapted mites/Tomato/Hormones/R data/ABA R data.csv", header = TRUE)
# trial as a factor
ABA.data$Trial <- factor(ABA.data$Trial)
str(ABA.data)
## 'data.frame':    14 obs. of  3 variables:
##  $ Trial    : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 2 2 2 ...
##  $ Treatment: Factor w/ 3 levels "control","TU",..: 1 1 3 3 2 2 2 1 1 3 ...
##  $ ng.g     : num  2463 2229 2204 2387 2302 ...

Formulate hypothesis

H0: There will be no difference in hormone levels between treatments.

HA: Hormone levels will differ by treatment.

Conduct data exploration

Outliers in the response variable (ng.g) within explanatory variables (Trial, Treatment).

ggplot(JA.data, aes(x = Trial, y = ng.g)) + geom_boxplot() + theme_classic()

ggplot(JA.data, aes(x = Treatment, y = ng.g)) + geom_boxplot() + theme_classic()

ggplot(JAI.data, aes(x = Trial, y = ng.g)) + geom_boxplot() + theme_classic()

ggplot(JAI.data, aes(x = Treatment, y = ng.g)) + geom_boxplot() + theme_classic()

ggplot(SA.data, aes(x = Trial, y = ng.g)) + geom_boxplot() + theme_classic()

ggplot(SA.data, aes(x = Treatment, y = ng.g)) + geom_boxplot() + theme_classic()

ggplot(ABA.data, aes(x = Trial, y = ng.g)) + geom_boxplot() + theme_classic()

ggplot(ABA.data, aes(x = Treatment, y = ng.g)) + geom_boxplot() + theme_classic()

JA outliers removed: Control 58.18414322 and TU 35.95505618.

JA-Ile outlier removed: TU 9.550561798.

ABA outlier removed: TU-A 1450.

Collinearity of the explanatory variables

Des 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 Treatment will be performed to test for reproducibility.

Zero inflation in Y

No

Are categorical covariates balanced?

summary(JA.data)
##  Trial   Treatment      ng.g        
##  1:7   control:3   Min.   :  6.881  
##  2:6   TU     :5   1st Qu.: 17.372  
##        TU-A   :5   Median :103.646  
##                    Mean   : 78.538  
##                    3rd Qu.:121.582  
##                    Max.   :139.447
summary(JAI.data)
##  Trial   Treatment      ng.g       
##  1:8   control:4   Min.   : 6.186  
##  2:6   TU     :5   1st Qu.: 7.245  
##        TU-A   :5   Median :14.323  
##                    Mean   :13.462  
##                    3rd Qu.:17.817  
##                    Max.   :21.611
summary(SA.data)
##  Trial   Treatment      ng.g       
##  1:8   control:4   Min.   : 515.6  
##  2:7   TU     :6   1st Qu.: 855.4  
##        TU-A   :5   Median :1419.1  
##                    Mean   :1320.7  
##                    3rd Qu.:1750.8  
##                    Max.   :2123.2
summary(ABA.data)
##  Trial   Treatment      ng.g     
##  1:7   control:4   Min.   :1794  
##  2:7   TU     :6   1st Qu.:2153  
##        TU-A   :4   Median :2265  
##                    Mean   :2365  
##                    3rd Qu.:2581  
##                    Max.   :3195

No, but close. Sample size is a little low.

Apply model

for JA

# fit linear model and display model fit information and ANOVA table
m.JA <- lm(ng.g ~ Treatment * Trial, data = JA.data)
summary(m.JA)
## 
## Call:
## lm(formula = ng.g ~ Treatment * Trial, data = JA.data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -61.439 -11.428   1.762  11.428  43.975 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)  
## (Intercept)            13.494     38.263   0.353    0.735  
## TreatmentTU           106.326     44.182   2.407    0.047 *
## TreatmentTU-A          72.687     44.182   1.645    0.144  
## Trial2                 -1.368     46.863  -0.029    0.978  
## TreatmentTU:Trial2      9.568     58.448   0.164    0.875  
## TreatmentTU-A:Trial2  -30.214     58.448  -0.517    0.621  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 38.26 on 7 degrees of freedom
## Multiple R-squared:  0.7041, Adjusted R-squared:  0.4927 
## F-statistic: 3.331 on 5 and 7 DF,  p-value: 0.07403
Anova(m.JA)
## Anova Table (Type II tests)
## 
## Response: ng.g
##                  Sum Sq Df F value  Pr(>F)  
## Treatment       21309.6  2  7.2775 0.01952 *
## Trial             273.7  1  0.1869 0.67848  
## Treatment:Trial  1005.1  2  0.3433 0.72075  
## Residuals       10248.5  7                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Calculate effect size and display
result.anova<-Anova(m.JA)
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: ng.g
##                  Sum Sq Df F value  Pr(>F) Part E Sq
## Treatment       21309.6  2  7.2775 0.01952   0.67525
## Trial             273.7  1  0.1869 0.67848   0.02601
## Treatment:Trial  1005.1  2  0.3433 0.72075   0.08932
## Residuals       10248.5  7
# perform post-hoc Tukey-Kramer test of contrasts
TukeyHSD(aov(ng.g ~ Treatment * Trial, data = JA.data))
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = ng.g ~ Treatment * Trial, data = JA.data)
## 
## $Treatment
##                   diff        lwr       upr     p adj
## TU-control   110.51766   28.22262 192.81271 0.0132957
## TU-A-control  60.96679  -21.32826 143.26183 0.1425884
## TU-A-TU      -49.55088 -120.82048  21.71872 0.1710765
## 
## $Trial
##          diff       lwr      upr     p adj
## 2-1 -8.967202 -59.30446 41.37006 0.6862128
## 
## $`Treatment:Trial`
##                            diff        lwr       upr     p adj
## TU:1-control:1       106.325762  -61.10101 273.75253 0.2691222
## TU-A:1-control:1      72.687439  -94.73933 240.11421 0.5988447
## control:2-control:1   -1.368151 -178.95106 176.21476 1.0000000
## TU:2-control:1       114.525264  -63.05765 292.10817 0.2575654
## TU-A:2-control:1      41.105554 -136.47736 218.68846 0.9407941
## TU-A:1-TU:1          -33.638322 -152.02693  84.75028 0.8764110
## control:2-TU:1      -107.693912 -240.05640  24.66857 0.1184562
## TU:2-TU:1              8.199502 -124.16298 140.56199 0.9998507
## TU-A:2-TU:1          -65.220207 -197.58269  67.14228 0.4861780
## control:2-TU-A:1     -74.055590 -206.41808  58.30690 0.3727813
## TU:2-TU-A:1           41.837825  -90.52466 174.20031 0.8257826
## TU-A:2-TU-A:1        -31.581885 -163.94437 100.78060 0.9336074
## TU:2-control:2       115.893415  -29.10242 260.88925 0.1266568
## TU-A:2-control:2      42.473705 -102.52213 187.46954 0.8632513
## TU-A:2-TU:2          -73.419710 -218.41555  71.57613 0.4615268
# plot interactions
interaction.plot(JA.data$Treatment, JA.data$Trial, JA.data$ng.g, type="l", leg.bty="o", leg.bg="grey95", lwd=2, ylab="ng.g", xlab="Treatment Strain", main="Treatment:Trial")

for JAI

# fit linear model and display model fit information and ANOVA table
m.JAI <- lm(ng.g ~ Treatment * Trial, data = JAI.data)
summary(m.JAI)
## 
## Call:
## lm(formula = ng.g ~ Treatment * Trial, data = JAI.data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.4202 -1.5552 -0.2696  1.5552  8.0052 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)  
## (Intercept)            7.3875     3.4551   2.138   0.0650 .
## TreatmentTU           10.3216     4.4605   2.314   0.0494 *
## TreatmentTU-A          6.2183     4.4605   1.394   0.2008  
## Trial2                 0.5846     4.8862   0.120   0.9077  
## TreatmentTU:Trial2     1.3533     6.6159   0.205   0.8430  
## TreatmentTU-A:Trial2  -1.9345     6.6159  -0.292   0.7774  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.886 on 8 degrees of freedom
## Multiple R-squared:  0.5836, Adjusted R-squared:  0.3233 
## F-statistic: 2.242 on 5 and 8 DF,  p-value: 0.1483
Anova(m.JAI)
## Anova Table (Type II tests)
## 
## Response: ng.g
##                  Sum Sq Df F value  Pr(>F)  
## Treatment       260.820  2  5.4622 0.03194 *
## Trial             0.490  1  0.0205 0.88966  
## Treatment:Trial   6.546  2  0.1371 0.87390  
## Residuals       190.999  8                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Calculate effect size and display
result.anova<-Anova(m.JAI)
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: ng.g
##                  Sum Sq Df F value  Pr(>F) Part E Sq
## Treatment       260.820  2  5.4622 0.03194   0.57727
## Trial             0.490  1  0.0205 0.88966   0.00256
## Treatment:Trial   6.546  2  0.1371 0.87390   0.03313
## Residuals       190.999  8
# perform post-hoc Tukey-Kramer test of contrasts
TukeyHSD(aov(ng.g ~ Treatment * Trial, data = JAI.data))
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = ng.g ~ Treatment * Trial, data = JAI.data)
## 
## $Treatment
##                   diff        lwr       upr     p adj
## TU-control   10.804454   1.438443 20.170466 0.0263632
## TU-A-control  5.386069  -3.979943 14.752081 0.2836177
## TU-A-TU      -5.418385 -14.248746  3.411976 0.2449779
## 
## $Trial
##          diff       lwr     upr     p adj
## 2-1 0.3763395 -5.708851 6.46153 0.8901206
## 
## $`Treatment:Trial`
##                           diff        lwr       upr     p adj
## TU:1-control:1      10.3215881  -5.975737 26.618913 0.2886223
## TU-A:1-control:1     6.2183473 -10.078978 22.515672 0.7303290
## control:2-control:1  0.5846355 -17.268189 18.437460 0.9999949
## TU:2-control:1      12.2595483  -5.593277 30.112373 0.2264182
## TU-A:2-control:1     4.8684466 -12.984378 22.721271 0.9068433
## TU-A:1-TU:1         -4.1032408 -18.680011 10.473530 0.8957360
## control:2-TU:1      -9.7369526 -26.034277  6.560372 0.3377418
## TU:2-TU:1            1.9379602 -14.359365 18.235285 0.9972474
## TU-A:2-TU:1         -5.4531415 -21.750466 10.844183 0.8155276
## control:2-TU-A:1    -5.6337117 -21.931037 10.663613 0.7963767
## TU:2-TU-A:1          6.0412011 -10.256124 22.338526 0.7508868
## TU-A:2-TU-A:1       -1.3499007 -17.647225 14.947424 0.9995088
## TU:2-control:2      11.6749128  -6.177912 29.527738 0.2630767
## TU-A:2-control:2     4.2838111 -13.569014 22.136636 0.9421277
## TU-A:2-TU:2         -7.3911017 -25.243927 10.461723 0.6672027
# plot interactions
interaction.plot(JAI.data$Treatment, JAI.data$Trial, JAI.data$ng.g, type="l", leg.bty="o", leg.bg="grey95", lwd=2, ylab="ng.g", xlab="Treatment Strain", main="Treatment:Trial")

for SA

# fit linear model and display model fit information and ANOVA table
m.SA <- lm(ng.g ~ Treatment * Trial, data = SA.data)
summary(m.SA)
## 
## Call:
## lm(formula = ng.g ~ Treatment * Trial, data = SA.data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -974.81 -271.99  -46.02  424.78  526.07 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)  
## (Intercept)            1041.7      403.5   2.581   0.0296 *
## TreatmentTU             713.2      521.0   1.369   0.2042  
## TreatmentTU-A           506.3      521.0   0.972   0.3565  
## Trial2                 -186.3      570.7  -0.326   0.7516  
## TreatmentTU:Trial2     -186.2      736.8  -0.253   0.8062  
## TreatmentTU-A:Trial2   -382.0      772.7  -0.494   0.6329  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 570.7 on 9 degrees of freedom
## Multiple R-squared:  0.3464, Adjusted R-squared:  -0.01678 
## F-statistic: 0.9538 on 5 and 9 DF,  p-value: 0.4925
Anova(m.SA)
## Anova Table (Type II tests)
## 
## Response: ng.g
##                  Sum Sq Df F value Pr(>F)
## Treatment        927865  2  1.4244 0.2901
## Trial            550234  1  1.6894 0.2260
## Treatment:Trial   80028  2  0.1229 0.8858
## Residuals       2931242  9
# Calculate effect size and display
result.anova<-Anova(m.SA)
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: ng.g
##                  Sum Sq Df F value  Pr(>F) Part E Sq
## Treatment        927865  2  1.4244 0.29010  0.240435
## Trial            550234  1  1.6894 0.22598  0.158046
## Treatment:Trial   80028  2  0.1229 0.88585  0.026576
## Residuals       2931242  9
# plot interactions
interaction.plot(SA.data$Treatment, SA.data$Trial, SA.data$ng.g, type="l", leg.bty="o", leg.bg="grey95", lwd=2, ylab="ng.g", xlab="Treatment Strain", main="Treatment:Trial")

for ABA

# fit linear model and display model fit information and ANOVA table
m.ABA <- lm(ng.g ~ Treatment * Trial, data = ABA.data)
summary(m.ABA)
## 
## Call:
## lm(formula = ng.g ~ Treatment * Trial, data = ABA.data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -523.34 -212.42   14.76  223.61  474.26 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           2345.80     245.45   9.557 1.19e-05 ***
## TreatmentTU           -307.67     316.87  -0.971    0.360    
## TreatmentTU-A          -50.15     347.11  -0.144    0.889    
## Trial2                  15.94     347.11   0.046    0.965    
## TreatmentTU:Trial2     667.02     448.12   1.488    0.175    
## TreatmentTU-A:Trial2   104.65     490.89   0.213    0.837    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 347.1 on 8 degrees of freedom
## Multiple R-squared:  0.4264, Adjusted R-squared:  0.06791 
## F-statistic: 1.189 on 5 and 8 DF,  p-value: 0.3931
Anova(m.ABA)
## Anova Table (Type II tests)
## 
## Response: ng.g
##                 Sum Sq Df F value Pr(>F)
## Treatment         2111  2  0.0088 0.9913
## Trial           385096  1  3.1962 0.1116
## Treatment:Trial 329345  2  1.3667 0.3086
## Residuals       963898  8
# Calculate effect size and display
result.anova<-Anova(m.ABA)
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: ng.g
##                 Sum Sq Df F value  Pr(>F) Part E Sq
## Treatment         2111  2  0.0088 0.99129  0.002185
## Trial           385096  1  3.1962 0.11162  0.285469
## Treatment:Trial 329345  2  1.3667 0.30861  0.254666
## Residuals       963898  8
# plot interactions
interaction.plot(ABA.data$Treatment, ABA.data$Trial, ABA.data$ng.g, type="l", leg.bty="o", leg.bg="grey95", lwd=2, ylab="ng.g", xlab="Treatment Strain", main="Treatment:Trial")

Validate models

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

# JA-Ile
JAI.data$m.JAI.fit <- fitted(m.JAI)    # fitted values
JAI.data$m.JAI.res <- rstandard(m.JAI) # Pearson residuals

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

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