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

The number of mites reamaining on Moneymaker after deposition of 50 TU or TU-A mites were counted on the leaflet 24h after application.

% retention was calculated.

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

dispersal.data <- read.csv("~/Lab Stuff/Adapted mites/Tomato/Dispersal Assay/R data/Dispersal R data.csv", header = TRUE)

# trial as a factor
dispersal.data$Trial <- factor(dispersal.data$Trial)

str(dispersal.data)
## 'data.frame':    32 obs. of  3 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 2 2 2 2 ...
##  $ Retention  : int  44 34 28 22 22 24 96 90 90 94 ...

Formulate hypothesis

H0: There will be no difference in dispersal of the mite strains.

HA: TU-A will be retained on the leaflet compared to TU mites.

Conduct data exploration

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

ggplot(dispersal.data, aes(x = Trial, y = Retention)) + geom_boxplot() + theme_classic()

ggplot(dispersal.data, aes(x = Mite.Strain, y = Retention)) + geom_boxplot() + theme_classic()

Only one outlier, probably real variability, I will keep it in the analysis and see if it causes problems in model validation.

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

Zero inflation in Y

No

Are categorical covariates balanced?

summary(dispersal.data)
##  Trial  Mite.Strain   Retention    
##  1:12   TU  :16     Min.   :12.00  
##  2:10   TU-A:16     1st Qu.:33.00  
##  3:10               Median :62.00  
##                     Mean   :61.12  
##                     3rd Qu.:92.50  
##                     Max.   :98.00

Not for Trial, but close

Yes for Mite.Strain

Apply model

# fit linear model and display model fit information and ANOVA table
m <- lm(Retention ~ Mite.Strain * Trial, data = dispersal.data)
summary(m)
## 
## Call:
## lm(formula = Retention ~ Mite.Strain * Trial, data = dispersal.data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -16.400  -2.700  -0.700   3.617  15.000 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              29.000      2.891  10.032 1.99e-10 ***
## Mite.StrainTU-A          63.333      4.088  15.491 1.21e-14 ***
## Trial2                   -0.600      4.288  -0.140    0.890    
## Trial3                    6.200      4.288   1.446    0.160    
## Mite.StrainTU-A:Trial2   -2.133      6.064  -0.352    0.728    
## Mite.StrainTU-A:Trial3   -6.133      6.064  -1.011    0.321    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.081 on 26 degrees of freedom
## Multiple R-squared:  0.9579, Adjusted R-squared:  0.9499 
## F-statistic: 118.4 on 5 and 26 DF,  p-value: < 2.2e-16
Anova(m)
## Anova Table (Type II tests)
## 
## Response: Retention
##                    Sum Sq Df  F value Pr(>F)    
## Mite.Strain       29524.5  1 588.7991 <2e-16 ***
## Trial               119.2  2   1.1889 0.3206    
## Mite.Strain:Trial    52.0  2   0.5188 0.6012    
## Residuals          1303.7 26                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Calculate effect size and display
result.anova<-Anova(m)
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: Retention
##                    Sum Sq Df  F value  Pr(>F) Part E Sq
## Mite.Strain       29524.5  1 588.7991 0.00000   0.95771
## Trial               119.2  2   1.1889 0.32057   0.08379
## Mite.Strain:Trial    52.0  2   0.5188 0.60124   0.03838
## Residuals          1303.7 26
# plot interactions

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

Validate model

dispersal.data$m.fit <- fitted(m)    # fitted values
dispersal.data$m.res <- rstandard(m) # 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.

ggplot(dispersal.data, aes(sample = m.res)) + geom_qq() +
  geom_abline(intercept = 0, slope = 1) + theme_classic()

Looks pretty good.

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 e constant across all fitted values.

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

Linearity - pretty good - confidence interval includes 0 for all values. Residuals farly equally dispersed.

Equal variance - no obvious problems

Residuals vs explanatory variables

Should be centered around 0, if not then model requires another explanatory variable(s), to account for observed variation.

ggplot(dispersal.data, aes(x = Mite.Strain, y = m.res)) + 
  geom_boxplot() + geom_hline(yintercept = 0) + theme_classic()

ggplot(dispersal.data, aes(x = Trial, y = m.res)) + 
  geom_boxplot() + geom_hline(yintercept = 0) + theme_classic()

ggplot(dispersal.data, aes(x = Mite.Strain:Trial, y = m.res)) + 
  geom_boxplot() + geom_hline(yintercept = 0) + theme_classic()

Looks good.

Model is valid and interpretation of ANOVA is good.

sessionInfo()
## R version 3.6.0 (2019-04-26)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 17763)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_Canada.1252  LC_CTYPE=English_Canada.1252   
## [3] LC_MONETARY=English_Canada.1252 LC_NUMERIC=C                   
## [5] LC_TIME=English_Canada.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] car_3.0-3     carData_3.0-2 ggplot2_3.1.1
## 
## loaded via a namespace (and not attached):
##  [1] zip_2.0.2         Rcpp_1.0.1        cellranger_1.1.0 
##  [4] pillar_1.4.1      compiler_3.6.0    plyr_1.8.4       
##  [7] forcats_0.4.0     tools_3.6.0       digest_0.6.19    
## [10] evaluate_0.14     tibble_2.1.1      gtable_0.3.0     
## [13] pkgconfig_2.0.2   rlang_0.3.4       openxlsx_4.1.0   
## [16] curl_3.3          yaml_2.2.0        haven_2.1.0      
## [19] xfun_0.7          rio_0.5.16        withr_2.1.2      
## [22] stringr_1.4.0     knitr_1.23        hms_0.4.2        
## [25] grid_3.6.0        data.table_1.12.2 readxl_1.3.1     
## [28] foreign_0.8-71    rmarkdown_1.13    magrittr_1.5     
## [31] scales_1.0.0      htmltools_0.3.6   abind_1.4-5      
## [34] colorspace_1.4-1  labeling_0.3      stringi_1.4.3    
## [37] lazyeval_0.2.2    munsell_0.5.0     crayon_1.3.4