Statistical Analysis

Analysis of mosquito trap data for the manuscript “Evaluation of modified autocidal gravid ovitraps for the control of container-breeding mosquitoes in Saint Augustine, FL”

Adam Rivers
USDA-ARS Genomics and Bioinformatics Unit

Overview

In this experiment three treatment and control areas were defined in St. Augustine, FL. Mosquito monitoring traps were placed in treatment and control areas at equal densities. In the treatment areas AGO traps were placed to test their effectiveness atreducing mosquito populations. After a 4 week period with no treatment, AGO kill traps were added and each area was monitored for 20 weeks.

Factors	comments
datetime	(24 levels, week of measurement)
type.y	(2 levels, the treatment of interest)
city_pair	(3 levels, geographic pair replicate)

All data are count data of mosquito populations

Analysis

Load the data frame.

library('tidyverse')
library('lubridate')
library('glmmTMB')

# Load raw data filtering NA's 
data <-as_tibble(read.csv("all_trap_data.csv.gz")) %>% filter(!is.na(count)) %>% filter(!is.na(organism))

# break dates into weeks
datetime <- mdy_hm(data$datetime)
data$week <- factor(cut(datetime, breaks="week"),ordered=TRUE)

Analysis of BG trap data

We will first focus on BG traps and female Aedes aegypti mosquitoes in BG traps

data_filtered <- filter(data, organism == "aedes_aegypti_female_count" & trap_type =="bg")

Visualize the count data

Data are highly skewed and there appear to be excess zeros. We will use the negative binomial with and withou zero inflation.

mean(data_filtered$count)

[1] 3.560185

sd(data_filtered$count)

[1] 6.065848

ggplot(data=data_filtered, aes(count)) + geom_histogram(bins=45)+ xlab("Mosquito count") + ylab("Number of traps")
ggsave("bg_hist.pdf")

Saving 7.29 x 4.51 in image

Explore the distribution of female Aedes aegypti over time in BG traps

ggplot(aes(x=week, y=count, group=treatment_type, color=treatment_type), data=data_filtered,) +
  geom_jitter(alpha=0.5) +
  theme(axis.text.x=element_text(angle =- 90, vjust = 0.5), legend.title=element_blank()) +
  ylab("Mosquito count") +
  xlab("Week")
ggsave("Aedes_aegypti_bg_weekly.pdf")

Saving 7.29 x 4.51 in image

Model 1

Estimate the number of mosquitoes in a trap as a function of the fixed effect of the treatment (treatment_type) and the random effect of the week.We will evaluate controlling for auto-correlation in time series data by adding auto-correlation structure of order 1 using week.

To use the Zero inflated model we will use the following R package:

Brooks ME, Kristensen K, van Benthem KJ, Magnusson A, Berg CW, Nielsen A, Skaug HJ, Maechler M, Bolker BM (2017). “glmmTMB Balances Speed and Flexibility Among Packages for Zero-inflated Generalized Linear Mixed Modeling.” The R Journal, 9(2), 378–400. https://journal.r-project.org/archive/2017/RJ-2017-066/index.html.

Model Comparison 1

modela <- glmmTMB(count ~ treatment_type + ar1(week-1|treatment_type), data = data_filtered, family = nbinom2)
modelb <- glmmTMB(count ~ treatment_type + ar1(week-1|treatment_type), data = data_filtered, zi=~factor(treatment_type), family = nbinom2)
modelc <- glmmTMB(count ~ treatment_type + (1|week), data = data_filtered, family = nbinom2)
modeld <- glmmTMB(count ~ treatment_type + (1|week), data = data_filtered, zi=~factor(treatment_type), family = nbinom2)
modele <- glmmTMB(count ~ treatment_type, data = data_filtered, family = nbinom2)
modelf <- glmmTMB(count ~ treatment_type, data = data_filtered, zi=~factor(treatment_type), family = nbinom2)

anova(modela,modelb,modelc,modeld,modele,modelf)

Data: data_filtered
Models:
modele: count ~ treatment_type, zi=~0, disp=~1
modelc: count ~ treatment_type + (1 | week), zi=~0, disp=~1
modela: count ~ treatment_type + ar1(week - 1 | treatment_type), zi=~0, disp=~1
modelf: count ~ treatment_type, zi=~factor(treatment_type), disp=~1
modeld: count ~ treatment_type + (1 | week), zi=~factor(treatment_type), disp=~1
modelb: count ~ treatment_type + ar1(week - 1 | treatment_type), zi=~factor(treatment_type), disp=~1
       Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
modele  3 1931.0 1943.2 -962.49   1925.0                             
modelc  4 1926.5 1942.8 -959.27   1918.5 6.4486      1   0.011104 *  
modela  5 1930.3 1950.7 -960.16   1920.3 0.0000      1   1.000000    
modelf  5 1929.6 1949.9 -959.78   1919.6 0.7598      0  < 2.2e-16 ***
modeld  6 1923.6 1948.0 -955.79   1911.6 7.9705      1   0.004755 ** 
modelb  7 1927.2 1955.7 -956.59   1913.2 0.0000      1   1.000000    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The model with the lowest BIC and AIC is model 1d. It is zero inflated but not autoregressive.

summary(modeld)

 Family: nbinom2  ( log )
Formula:          count ~ treatment_type + (1 | week)
Zero inflation:         ~factor(treatment_type)
Data: data_filtered

     AIC      BIC   logLik deviance df.resid 
  1923.6   1948.0   -955.8   1911.6      426 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 week   (Intercept) 0.2051   0.4529  
Number of obs: 432, groups:  week, 24

Overdispersion parameter for nbinom2 family (): 0.77 

Conditional model:
                        Estimate Std. Error z value Pr(>|z|)    
(Intercept)               1.1421     0.1638   6.974 3.09e-12 ***
treatment_typetreatment   0.5162     0.1671   3.089  0.00201 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Zero-inflation model:
                                Estimate Std. Error z value Pr(>|z|)   
(Intercept)                      -1.4080     0.5023  -2.803  0.00506 **
factor(treatment_type)treatment   0.3795     0.4294   0.884  0.37684   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

library('DHARMa')
modeld_simres <- simulateResiduals(modeld)
plotResiduals(modeld_simres, quantreg = T)

Model 1d is the best supported model. It indicates a non-significant positive association with traps and mosquitoes overall in the ZINB model and a significant positive relation in the conditional model. This suggests the traps were not effective at reducing female Aedes aegypti.

Model 2

Estimate the number of mosquitoes in a trap as a function of the fixed effect of the treatment (treatment_type) and the random effect of the week. Use NB as the family since the data is count data and may be over-dispersed. Account for auto-correlation in time series data by adding autocorrelation structure of order 1 using datetime. This is different from model 1 because the 3 site pairs are broken out as a factor.

model2 <- glmmTMB(count ~ treatment_type*site_pair + (1|week), data = data_filtered, zi=~treatment_type*site_pair, family = nbinom2)


summary(model1)

 Family: nbinom2  ( log )
Formula:          count ~ treatment_type * site_pair + (1 | week)
Zero inflation:         ~treatment_type * site_pair
Data: data_filtered

     AIC      BIC   logLik deviance df.resid 
  1725.8   1782.8   -848.9   1697.8      418 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 week   (Intercept) 0.3517   0.593   
Number of obs: 432, groups:  week, 24

Overdispersion parameter for nbinom2 family ():  1.8 

Conditional model:
                                          Estimate Std. Error z value Pr(>|z|)    
(Intercept)                                1.46075    0.17350   8.419  < 2e-16 ***
treatment_typetreatment                    0.71316    0.16152   4.415 1.01e-05 ***
site_pairsasnorth                         -1.60595    0.30612  -5.246 1.55e-07 ***
site_pairsassouth                         -0.01161    0.19607  -0.059    0.953    
treatment_typetreatment:site_pairsasnorth  0.32859    0.34823   0.944    0.345    
treatment_typetreatment:site_pairsassouth -2.37506    0.44165  -5.378 7.54e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Zero-inflation model:
                                          Estimate Std. Error z value Pr(>|z|)   
(Intercept)                                -2.4210     0.7613  -3.180  0.00147 **
treatment_typetreatment                    -1.1338     1.6662  -0.680  0.49621   
site_pairsasnorth                           1.9447     0.9161   2.123  0.03377 * 
site_pairsassouth                           1.6613     0.7997   2.078  0.03775 * 
treatment_typetreatment:site_pairsasnorth  -0.5154     1.8879  -0.273  0.78484   
treatment_typetreatment:site_pairsassouth   2.4422     1.7565   1.390  0.16441   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

model2_simres <- simulateResiduals(model2)
plotResiduals(modeld_simres, quantreg = T)

There is a positive association with traps and mosquitoes in downtown and a negative relationship in the north and south. This model suggests that there is not a strong consistient relationship between traps and mosquitoes.

SAGO trap data

We will first focus on BG traps and female Aedes aegypti mosquitoes in BG traps

data_filtered2 <- filter(data, organism == "aedes_aegypti_female_count" & trap_type =="sago")

Visualize the count data.

Data are highly skewed and there appear to be many excess zeros. I will evaluate zero inflated methods.

mean(data_filtered2$count)

[1] 0.06468531

sd(data_filtered2$count)

[1] 0.2761478

ggplot(data=data_filtered2, aes(count)) + geom_histogram(bins=6)+ xlab("Mosquito count") + ylab("Number of traps")
ggsave("sago_hist.pdf")

Saving 7.29 x 4.51 in image

There are many Zeros in the data. A zero inflated Negative Binomial or zero inflated Poisson model may be the most appropriate.

Explore the distribution of female Aedes aegypti over time in SAGO traps

ggplot( aes(x=week, y=count, group=treatment_type, color=treatment_type), data=data_filtered2,) + geom_jitter(alpha=0.5)+ theme(axis.text.x=element_text(angle =- 90, vjust = 0.5), legend.title=element_blank()) + ylab("Mosquito count") + xlab("Week")

ggsave("Aedes_aegypti_sago_weekly.pdf")

Saving 7.29 x 4.51 in image

Zero inflated regression Models

To use the Zero inflated model we will use the follwoing R package:

Brooks ME, Kristensen K, van Benthem KJ, Magnusson A, Berg CW, Nielsen A, Skaug HJ, Maechler M, Bolker BM (2017). “glmmTMB Balances Speed and Flexibility Among Packages for Zero-inflated Generalized Linear Mixed Modeling.” The R Journal, 9(2), 378–400. https://journal.r-project.org/archive/2017/RJ-2017-066/index.html.

Model 3

Glmm’s negative binomial link function.

library(glmmTMB)
model3a <- glmmTMB(count ~ treatment_type, data = data_filtered2, family = nbinom2)
model3b <- glmmTMB(count ~ treatment_type , data = data_filtered2, zi=~treatment_type, family = nbinom2)
model3c <- glmmTMB(count ~ treatment_type + ar1(week-1|treatment_type), data = data_filtered2, family = nbinom2)
model3d <- glmmTMB(count ~ treatment_type + ar1(week-1|treatment_type), data = data_filtered2, zi=~treatment_type, family = nbinom2)
model3e <- glmmTMB(count ~ treatment_type + (1|week), data = data_filtered2, zi=~treatment_type, family = nbinom2)
anova(model3a, model3b, model3c, model3d, model3e)

Data: data_filtered2
Models:
model3a: count ~ treatment_type, zi=~0, disp=~1
model3b: count ~ treatment_type, zi=~treatment_type, disp=~1
model3c: count ~ treatment_type + ar1(week - 1 | treatment_type), zi=~0, disp=~1
model3e: count ~ treatment_type + (1 | week), zi=~treatment_type, disp=~1
model3d: count ~ treatment_type + ar1(week - 1 | treatment_type), zi=~treatment_type, disp=~1
        Df    AIC    BIC  logLik deviance   Chisq Chi Df Pr(>Chisq)    
model3a  3 1671.2 1689.7 -832.61   1665.2                              
model3b  5 1674.5 1705.2 -832.26   1664.5  0.6972      2   0.705681    
model3c  5 1633.1 1663.8 -811.56   1623.1 41.4041      0  < 2.2e-16 ***
model3e  6 1625.4 1662.2 -806.69   1613.4  9.7352      1   0.001808 ** 
model3d  7 1636.8 1679.8 -811.39   1622.8  0.0000      1   1.000000    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Model 3e has the lowest AIC/BIC scores.

summary(model3e)

 Family: nbinom2  ( log )
Formula:          count ~ treatment_type + (1 | week)
Zero inflation:         ~treatment_type
Data: data_filtered2

     AIC      BIC   logLik deviance df.resid 
  1625.4   1662.2   -806.7   1613.4     3426 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 week   (Intercept) 0.5593   0.7478  
Number of obs: 3432, groups:  week, 24

Overdispersion parameter for nbinom2 family (): 1.22 

Conditional model:
                        Estimate Std. Error z value Pr(>|z|)    
(Intercept)              -2.5649     0.4784  -5.362 8.24e-08 ***
treatment_typetreatment  -0.6047     0.4627  -1.307    0.191    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Zero-inflation model:
                        Estimate Std. Error z value Pr(>|z|)
(Intercept)               -1.102      1.759  -0.626    0.531
treatment_typetreatment  -15.141   3476.414  -0.004    0.997

model3e_simres <- simulateResiduals(model3e)
plotResiduals(modeld_simres, quantreg = T)

There is not a significant effect of treatment on counts of female Aedes aegypti in SAGO traps