“The impact of sound lure frequency and habitat type on male Aedes albopictus capture rates with the Male Aedes Sound Trap”

Preliminaries:

  1. Make sure all required R packages are installed and loaded.
  2. Set an appropriate working directory where the dataset “Swan_et_al_2020_data.xlsx” is stored.

Required packages:

library(qtl)
library(ggplot2)
library(pscl)
library(boot)
library(psych) 
library(stats) 
library(MASS)
library(car)
library(lmtest)
library(lme4)
library(tidyverse)
library(plyr)
library(tibble)
library(readxl)
library(aods3)
library(emmeans)

Data import

Read in the xlsx file, renaming it “data_all”.

Code used to analyse experiments 1 (“freq_01”), 2 (“freq_02”) and 3 (“habitat_01”).

data_all <- "C:/Users/jc486631/OneDrive - James Cook University/Chapter 2 - MAST/Manuscript/co-authors comments/Further revisions/July revisions/submission/data/Swan_et_al_2020_data.xlsx"

#lists the 3 datasheets - "freq_01"    "freq_02"    "habitat_01"

excel_sheets(path = data_all)
## [1] "freq_01"    "freq_02"    "habitat_01"

Experiments 1 (“freq_01”) data

  1. to determine the most effective frequency (450, 500, 550 & 600 Hz) for capturing male Ae. albopictus.
#we will now read in exp 1 - "freq_01" data and call it "t1"

t1 <- read_excel(path = data_all, sheet = "freq_01")
str(t1)
## Classes 'tbl_df', 'tbl' and 'data.frame':    96 obs. of  7 variables:
##  $ square   : chr  "A" "B" "C" "A" ...
##  $ position : chr  "four" "three" "three" "one" ...
##  $ station  : chr  "Afour" "Bthree" "Cthree" "Aone" ...
##  $ replicate: chr  "eight" "eight" "eight" "five" ...
##  $ trap     : chr  "a450" "a450" "a450" "a450" ...
##  $ nday     : chr  "eight" "eight" "eight" "five" ...
##  $ daily_sum: num  0 2 0 0 2 3 0 3 0 1 ...

A few summary statistics about each trap is given below:

t1_stats <- describeBy(t1, t1$trap)
t1_stats
## 
##  Descriptive statistics by group 
## group: a450
##            vars  n mean   sd median trimmed  mad min  max range skew kurtosis
## square*       1 24  NaN   NA     NA     NaN   NA Inf -Inf  -Inf   NA       NA
## position*     2 24  NaN   NA     NA     NaN   NA Inf -Inf  -Inf   NA       NA
## station*      3 24  NaN   NA     NA     NaN   NA Inf -Inf  -Inf   NA       NA
## replicate*    4 24  NaN   NA     NA     NaN   NA Inf -Inf  -Inf   NA       NA
## trap*         5 24  NaN   NA     NA     NaN   NA Inf -Inf  -Inf   NA       NA
## nday*         6 24  NaN   NA     NA     NaN   NA Inf -Inf  -Inf   NA       NA
## daily_sum     7 24 0.96 1.08      1    0.85 1.48   0    3     3 0.67    -0.99
##              se
## square*      NA
## position*    NA
## station*     NA
## replicate*   NA
## trap*        NA
## nday*        NA
## daily_sum  0.22
## ------------------------------------------------------------ 
## group: b500
##            vars  n mean   sd median trimmed  mad min  max range skew kurtosis
## square*       1 24  NaN   NA     NA     NaN   NA Inf -Inf  -Inf   NA       NA
## position*     2 24  NaN   NA     NA     NaN   NA Inf -Inf  -Inf   NA       NA
## station*      3 24  NaN   NA     NA     NaN   NA Inf -Inf  -Inf   NA       NA
## replicate*    4 24  NaN   NA     NA     NaN   NA Inf -Inf  -Inf   NA       NA
## trap*         5 24  NaN   NA     NA     NaN   NA Inf -Inf  -Inf   NA       NA
## nday*         6 24  NaN   NA     NA     NaN   NA Inf -Inf  -Inf   NA       NA
## daily_sum     7 24 3.08 3.32      2    2.55 2.22   0   13    13  1.4     1.41
##              se
## square*      NA
## position*    NA
## station*     NA
## replicate*   NA
## trap*        NA
## nday*        NA
## daily_sum  0.68
## ------------------------------------------------------------ 
## group: c550
##            vars  n mean   sd median trimmed  mad min  max range skew kurtosis
## square*       1 24  NaN   NA     NA     NaN   NA Inf -Inf  -Inf   NA       NA
## position*     2 24  NaN   NA     NA     NaN   NA Inf -Inf  -Inf   NA       NA
## station*      3 24  NaN   NA     NA     NaN   NA Inf -Inf  -Inf   NA       NA
## replicate*    4 24  NaN   NA     NA     NaN   NA Inf -Inf  -Inf   NA       NA
## trap*         5 24  NaN   NA     NA     NaN   NA Inf -Inf  -Inf   NA       NA
## nday*         6 24  NaN   NA     NA     NaN   NA Inf -Inf  -Inf   NA       NA
## daily_sum     7 24 4.25 4.22      2     3.8 2.97   0   14    14 0.72    -0.76
##              se
## square*      NA
## position*    NA
## station*     NA
## replicate*   NA
## trap*        NA
## nday*        NA
## daily_sum  0.86
## ------------------------------------------------------------ 
## group: d600
##            vars  n mean  sd median trimmed  mad min  max range skew kurtosis
## square*       1 24  NaN  NA     NA     NaN   NA Inf -Inf  -Inf   NA       NA
## position*     2 24  NaN  NA     NA     NaN   NA Inf -Inf  -Inf   NA       NA
## station*      3 24  NaN  NA     NA     NaN   NA Inf -Inf  -Inf   NA       NA
## replicate*    4 24  NaN  NA     NA     NaN   NA Inf -Inf  -Inf   NA       NA
## trap*         5 24  NaN  NA     NA     NaN   NA Inf -Inf  -Inf   NA       NA
## nday*         6 24  NaN  NA     NA     NaN   NA Inf -Inf  -Inf   NA       NA
## daily_sum     7 24 4.75 4.1      4     4.5 4.45   0   12    12 0.43    -1.23
##              se
## square*      NA
## position*    NA
## station*     NA
## replicate*   NA
## trap*        NA
## nday*        NA
## daily_sum  0.84
#the average daily catch of male Aedes albopictus between 450 - 500, 550 and 600 appears to differ considerably.
#Statistical investigations using generalised linear mixed models (GLMM) will need to occur.  

GLMM - Model selection: Poisson “freq_01”

Experimental design - 1 fixed effect (trap) + 2 random factors (1|nday) + (1|station)

alboglmm.p_t1 <- glmer(daily_sum ~ trap + (1|nday) + (1|station), family = poisson, data=t1)
summary(alboglmm.p_t1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: daily_sum ~ trap + (1 | nday) + (1 | station)
##    Data: t1
## 
##      AIC      BIC   logLik deviance df.resid 
##    422.8    438.2   -205.4    410.8       90 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3772 -0.9410 -0.1999  0.4644  3.8610 
## 
## Random effects:
##  Groups  Name        Variance Std.Dev.
##  station (Intercept) 0.5366   0.7326  
##  nday    (Intercept) 0.1133   0.3366  
## Number of obs: 96, groups:  station, 12; nday, 8
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.3723     0.3240  -1.149    0.251    
## trapb500      1.2106     0.2382   5.082 3.74e-07 ***
## trapc550      1.4838     0.2287   6.489 8.65e-11 ***
## trapd600      1.6971     0.2300   7.377 1.61e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr) trp500 trp550
## trapb500 -0.560              
## trapc550 -0.577  0.782       
## trapd600 -0.593  0.784  0.813
Anova(alboglmm.p_t1)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: daily_sum
##       Chisq Df Pr(>Chisq)    
## trap 57.702  3   1.82e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Testing for overdispersion with the poisson model

overdis.P_t1<-gof(alboglmm.p_t1)
overdis.P_t1
## D  = 136.6373, df = 90, P(>D) = 0.001114409
## X2 = 123.4842, df = 90, P(>X2) = 0.01106384
sum(residuals(alboglmm.p_t1,"pearson")^2)
## [1] 123.4842

GLMM. Model selection: negative binomial “freq_01”

alboglmm.nb_t1 <- glmer.nb(daily_sum ~ trap + (1|nday) + (1|station), data=t1)
summary(alboglmm.nb_t1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Negative Binomial(4.2387)  ( log )
## Formula: daily_sum ~ trap + (1 | nday) + (1 | station)
##    Data: t1
## 
##      AIC      BIC   logLik deviance df.resid 
##    410.7    428.6   -198.3    396.7       89 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.49325 -0.69618 -0.09611  0.39240  3.00351 
## 
## Random effects:
##  Groups  Name        Variance Std.Dev.
##  station (Intercept) 0.5397   0.7346  
##  nday    (Intercept) 0.1292   0.3595  
## Number of obs: 96, groups:  station, 12; nday, 8
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.4013     0.3529  -1.137    0.255    
## trapb500      1.2206     0.2919   4.182 2.89e-05 ***
## trapc550      1.4831     0.2838   5.227 1.73e-07 ***
## trapd600      1.7445     0.2842   6.137 8.38e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr) trp500 trp550
## trapb500 -0.578              
## trapc550 -0.590  0.696       
## trapd600 -0.611  0.707  0.726
Anova(alboglmm.nb_t1)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: daily_sum
##       Chisq Df Pr(>Chisq)    
## trap 39.405  3  1.424e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Testing for overdispersion with the negative binomial model

overdis.NB_t1 <- gof(alboglmm.nb_t1)
overdis.NB_t1
## D  = 87.2148, df = 89, P(>D) = 0.5337226
## X2 = 70.3108, df = 89, P(>X2) = 0.9281349
sum(residuals(alboglmm.nb_t1,"pearson")^2)
## [1] 70.31078

Comparing models

sum(residuals(alboglmm.p_t1,"pearson")^2)
## [1] 123.4842
sum(residuals(alboglmm.nb_t1,"pearson")^2)
## [1] 70.31078
#The negative binomial model (70) is less overdispersed than the poisson model (123). 

#comparing the AIC weights between models
AIC(alboglmm.p_t1 ,alboglmm.nb_t1)
##                df      AIC
## alboglmm.p_t1   6 422.8062
## alboglmm.nb_t1  7 410.6734
#alboglmm.nb has a lower AIC (410), compared with alboglmm.p_t1 (422) and is less overdispersed than alboglmm.p_t1. 

“freq_01” Tukey post hoc test for significance - negative binomial model

tukeyTrap_t1<-emmeans(alboglmm.nb_t1, list(pairwise ~ trap), adjust = "tukey")
confint(tukeyTrap_t1, level = 0.95)
## $`emmeans of trap`
##  trap emmean    SE  df asymp.LCL asymp.UCL
##  a450 -0.401 0.353 Inf   -1.2803     0.478
##  b500  0.819 0.301 Inf    0.0691     1.569
##  c550  1.082 0.295 Inf    0.3475     1.816
##  d600  1.343 0.288 Inf    0.6265     2.060
## 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 4 estimates 
## 
## $`pairwise differences of trap`
##  contrast    estimate    SE  df asymp.LCL asymp.UCL
##  a450 - b500   -1.221 0.292 Inf    -1.971   -0.4707
##  a450 - c550   -1.483 0.284 Inf    -2.212   -0.7541
##  a450 - d600   -1.744 0.284 Inf    -2.475   -1.0143
##  b500 - c550   -0.262 0.225 Inf    -0.840    0.3147
##  b500 - d600   -0.524 0.221 Inf    -1.091    0.0432
##  c550 - d600   -0.261 0.210 Inf    -0.802    0.2789
## 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## Conf-level adjustment: tukey method for comparing a family of 4 estimates
tukeyTrap_t1
## $`emmeans of trap`
##  trap emmean    SE  df asymp.LCL asymp.UCL
##  a450 -0.401 0.353 Inf   -1.2803     0.478
##  b500  0.819 0.301 Inf    0.0691     1.569
##  c550  1.082 0.295 Inf    0.3475     1.816
##  d600  1.343 0.288 Inf    0.6265     2.060
## 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 4 estimates 
## 
## $`pairwise differences of trap`
##  contrast    estimate    SE  df z.ratio p.value
##  a450 - b500   -1.221 0.292 Inf -4.182  0.0002 
##  a450 - c550   -1.483 0.284 Inf -5.227  <.0001 
##  a450 - d600   -1.744 0.284 Inf -6.137  <.0001 
##  b500 - c550   -0.262 0.225 Inf -1.168  0.6469 
##  b500 - d600   -0.524 0.221 Inf -2.373  0.0823 
##  c550 - d600   -0.261 0.210 Inf -1.243  0.5994 
## 
## Results are given on the log (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 4 estimates

Tukey post hoc shows that traps with sound lures set to 450 Hz captured significantly less male Ae. albopictus than traps with sound lures set to 500, 550 and 600 Hz. No significant differences in average daily catch of male Ae. albopictus were found between 500, 550 and 600 Hz.

This result is displayed in Figure 4A of the Swan et al. 2020 manuscript (graph produced in GraphPad Prism 8).

Experiment 2 - “freq_02”

  1. to determine the most effective frequency (450, 600, 650 & 700 Hz) for capturing male Ae. albopictus.
#lists the 3 datasheets - "freq_01"    "freq_02"    "habitat_01"
excel_sheets(path = data_all)
## [1] "freq_01"    "freq_02"    "habitat_01"
#select "freq_02" and we will call it "t2"

t2 <- read_excel(path = data_all, sheet = "freq_02")
str(t2)
## Classes 'tbl_df', 'tbl' and 'data.frame':    144 obs. of  7 variables:
##  $ square   : chr  "A" "B" "C" "A" ...
##  $ position : chr  "two" "two" "two" "two" ...
##  $ station  : chr  "Atwo" "Btwo" "Ctwo" "Atwo" ...
##  $ replicate: chr  "eight" "eight" "eight" "eleven" ...
##  $ trap     : chr  "c650" "c650" "c650" "d700" ...
##  $ nday     : chr  "eight" "eight" "eight" "eleven" ...
##  $ daily_sum: num  0 1 1 0 1 0 0 0 23 3 ...
t2_stats <- describeBy(t2, t2$trap)
t2_stats
## 
##  Descriptive statistics by group 
## group: a450
##            vars  n mean   sd median trimmed mad min  max range skew kurtosis
## square*       1 36  NaN   NA     NA     NaN  NA Inf -Inf  -Inf   NA       NA
## position*     2 36  NaN   NA     NA     NaN  NA Inf -Inf  -Inf   NA       NA
## station*      3 36  NaN   NA     NA     NaN  NA Inf -Inf  -Inf   NA       NA
## replicate*    4 36  NaN   NA     NA     NaN  NA Inf -Inf  -Inf   NA       NA
## trap*         5 36  NaN   NA     NA     NaN  NA Inf -Inf  -Inf   NA       NA
## nday*         6 36  NaN   NA     NA     NaN  NA Inf -Inf  -Inf   NA       NA
## daily_sum     7 36 1.03 1.86      0    0.67   0   0    9     9  2.5     7.18
##              se
## square*      NA
## position*    NA
## station*     NA
## replicate*   NA
## trap*        NA
## nday*        NA
## daily_sum  0.31
## ------------------------------------------------------------ 
## group: b600
##            vars  n mean   sd median trimmed  mad min  max range skew kurtosis
## square*       1 36  NaN   NA     NA     NaN   NA Inf -Inf  -Inf   NA       NA
## position*     2 36  NaN   NA     NA     NaN   NA Inf -Inf  -Inf   NA       NA
## station*      3 36  NaN   NA     NA     NaN   NA Inf -Inf  -Inf   NA       NA
## replicate*    4 36  NaN   NA     NA     NaN   NA Inf -Inf  -Inf   NA       NA
## trap*         5 36  NaN   NA     NA     NaN   NA Inf -Inf  -Inf   NA       NA
## nday*         6 36  NaN   NA     NA     NaN   NA Inf -Inf  -Inf   NA       NA
## daily_sum     7 36 3.19 5.53      1    1.97 1.48   0   23    23 2.07     3.64
##              se
## square*      NA
## position*    NA
## station*     NA
## replicate*   NA
## trap*        NA
## nday*        NA
## daily_sum  0.92
## ------------------------------------------------------------ 
## group: c650
##            vars  n mean   sd median trimmed  mad min  max range skew kurtosis
## square*       1 36  NaN   NA     NA     NaN   NA Inf -Inf  -Inf   NA       NA
## position*     2 36  NaN   NA     NA     NaN   NA Inf -Inf  -Inf   NA       NA
## station*      3 36  NaN   NA     NA     NaN   NA Inf -Inf  -Inf   NA       NA
## replicate*    4 36  NaN   NA     NA     NaN   NA Inf -Inf  -Inf   NA       NA
## trap*         5 36  NaN   NA     NA     NaN   NA Inf -Inf  -Inf   NA       NA
## nday*         6 36  NaN   NA     NA     NaN   NA Inf -Inf  -Inf   NA       NA
## daily_sum     7 36 3.72 5.95    1.5    2.47 2.22   0   28    28 2.68     7.31
##              se
## square*      NA
## position*    NA
## station*     NA
## replicate*   NA
## trap*        NA
## nday*        NA
## daily_sum  0.99
## ------------------------------------------------------------ 
## group: d700
##            vars  n mean   sd median trimmed mad min  max range skew kurtosis
## square*       1 36  NaN   NA     NA     NaN  NA Inf -Inf  -Inf   NA       NA
## position*     2 36  NaN   NA     NA     NaN  NA Inf -Inf  -Inf   NA       NA
## station*      3 36  NaN   NA     NA     NaN  NA Inf -Inf  -Inf   NA       NA
## replicate*    4 36  NaN   NA     NA     NaN  NA Inf -Inf  -Inf   NA       NA
## trap*         5 36  NaN   NA     NA     NaN  NA Inf -Inf  -Inf   NA       NA
## nday*         6 36  NaN   NA     NA     NaN  NA Inf -Inf  -Inf   NA       NA
## daily_sum     7 36 2.06 3.11      0    1.53   0   0   10    10 1.34     0.49
##              se
## square*      NA
## position*    NA
## station*     NA
## replicate*   NA
## trap*        NA
## nday*        NA
## daily_sum  0.52
#the average daily catch of male Aedes albopictus between 450 - 600, 650 and 700 appears to differ considerably.
#less so between 600 - 700 Hz.
#Statistical investigations using generalised linear mixed models (GLMM) will need to occur.  

GLMM. Model selection: Poisson “freq_02”

Experimental design - 1 fixed effect (trap) + 2 random factors (1|nday) + (1|station)

alboglmm.p_t2 <- glmer(daily_sum ~ trap + (1|nday)+(1|station),family = poisson, data=t2)
summary(alboglmm.p_t2)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: daily_sum ~ trap + (1 | nday) + (1 | station)
##    Data: t2
## 
##      AIC      BIC   logLik deviance df.resid 
##    572.2    590.0   -280.1    560.2      138 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1883 -0.7925 -0.3955  0.4930  7.6257 
## 
## Random effects:
##  Groups  Name        Variance Std.Dev.
##  nday    (Intercept) 1.3026   1.1413  
##  station (Intercept) 0.9665   0.9831  
## Number of obs: 144, groups:  nday, 12; station, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.8312     0.4735  -1.755 0.079205 .  
## trapb600      1.1347     0.1868   6.074 1.25e-09 ***
## trapc650      1.2832     0.1837   6.987 2.81e-12 ***
## trapd700      0.6933     0.1989   3.485 0.000492 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr) trp600 trp650
## trapb600 -0.298              
## trapc650 -0.303  0.769       
## trapd700 -0.280  0.708  0.722
Anova(alboglmm.p_t2)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: daily_sum
##      Chisq Df Pr(>Chisq)    
## trap  57.9  3  1.651e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Test for overdispersion
overdis.P_t2<-gof(alboglmm.p_t2)
overdis.P_t2
## D  = 235.6984, df = 138, P(>D) = 4.338784e-07
## X2 = 268.7795, df = 138, P(>X2) = 1.86077e-10
sum(residuals(alboglmm.p_t2,"pearson")^2)
## [1] 268.7795

GLMM. Model selection: negative binomial “freq_02”

Experimental design - 1 fixed effect (trap) + 2 random factors (1|nday) + (1|station)

alboglmm.nb_t2 <- glmer.nb(daily_sum ~ trap + (1|nday) + (1|station), data=t2)
summary(alboglmm.nb_t2)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Negative Binomial(1.5181)  ( log )
## Formula: daily_sum ~ trap + (1 | nday) + (1 | station)
##    Data: t2
## 
##      AIC      BIC   logLik deviance df.resid 
##    524.4    545.2   -255.2    510.4      137 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.0033 -0.6208 -0.3595  0.4671  4.2343 
## 
## Random effects:
##  Groups  Name        Variance Std.Dev.
##  nday    (Intercept) 1.1006   1.0491  
##  station (Intercept) 0.7689   0.8769  
## Number of obs: 144, groups:  nday, 12; station, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.7757     0.4755  -1.631 0.102810    
## trapb600      1.0507     0.3071   3.421 0.000624 ***
## trapc650      1.4379     0.3112   4.621 3.82e-06 ***
## trapd700      0.6998     0.3163   2.213 0.026922 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr) trp600 trp650
## trapb600 -0.408              
## trapc650 -0.429  0.619       
## trapd700 -0.395  0.607  0.597
Anova(alboglmm.nb_t2)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: daily_sum
##      Chisq Df Pr(>Chisq)    
## trap 22.91  3  4.216e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#testing overdispersion
overdis.NB_t2 <- gof(alboglmm.nb_t2)
overdis.NB_t2
## D  = 117.6902, df = 137, P(>D) = 0.8821049
## X2 = 112.9936, df = 137, P(>X2) = 0.9337036
sum(residuals(alboglmm.nb_t2,"pearson")^2)
## [1] 112.9936

Comparing models

#comparing overdispersion

sum(residuals(alboglmm.nb_t2,"pearson")^2)
## [1] 112.9936
sum(residuals(alboglmm.p_t2,"pearson")^2)
## [1] 268.7795
#alboglmm.nb_t2 is far less overdispersed (112) than alboglmm.p_t2 (268)

#compare AIC weights
AIC(alboglmm.nb_t2, alboglmm.p_t2)
##                df      AIC
## alboglmm.nb_t2  7 524.3878
## alboglmm.p_t2   6 572.2269
#alboglmm.nb_t2 has both a lower AIC value (524) to alboglmm.p_t2 (572) and is less overdispersed than the poisson model. 

“freq_02” Tukey post hoc test for significance - negative binomial model

#Tukey post hoc test for significance.

tukeytrap_02<-emmeans(alboglmm.nb_t2, list(pairwise ~ trap), adjust = "tukey")
tukeytrap_02
## $`emmeans of trap`
##  trap emmean    SE  df asymp.LCL asymp.UCL
##  a450 -0.776 0.476 Inf    -1.960     0.409
##  b600  0.275 0.449 Inf    -0.843     1.393
##  c650  0.662 0.443 Inf    -0.441     1.765
##  d700 -0.076 0.455 Inf    -1.210     1.058
## 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 4 estimates 
## 
## $`pairwise differences of trap`
##  contrast    estimate    SE  df z.ratio p.value
##  a450 - b600   -1.051 0.307 Inf -3.421  0.0035 
##  a450 - c650   -1.438 0.311 Inf -4.621  <.0001 
##  a450 - d700   -0.700 0.316 Inf -2.213  0.1197 
##  b600 - c650   -0.387 0.270 Inf -1.435  0.4775 
##  b600 - d700    0.351 0.276 Inf  1.270  0.5820 
##  c650 - d700    0.738 0.282 Inf  2.619  0.0437 
## 
## Results are given on the log (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 4 estimates

Tukey post hoc test shows that traps with sound lures set to 450 Hz captured significantly less male Ae. albopictus than traps with sound lures set to 600, 650 and 700 Hz. Traps with sound lures set to 700 Hz captured significantly less male Ae. albopictus compared to traps with sound lures set to 650 Hz.

This result is displayed in Figure 4B of the Swan et al. 2020 manuscript (graph produced in GraphPad Prism 8).

Experiment 3 - habitat_01

#lists the 3 datasheets - "freq_01"    "freq_02"    "habitat_01"
excel_sheets(path = data_all)
## [1] "freq_01"    "freq_02"    "habitat_01"
#select "freq_02" and we will call it "t2"

t3 <- read_excel(path = data_all, sheet = "habitat_01")

t3_stats <- describeBy(t3, t3$habitat_type)
t3_stats
## 
##  Descriptive statistics by group 
## group: house
##               vars  n mean   sd median trimmed  mad min  max range skew
## station*         1 28  NaN   NA     NA     NaN   NA Inf -Inf  -Inf   NA
## nday             2 28 4.00 2.04      4    4.00 2.97   1    7     6 0.00
## habitat_type*    3 28  NaN   NA     NA     NaN   NA Inf -Inf  -Inf   NA
## daily_sum        4 28 0.32 0.86      0    0.12 0.00   0    4     4 3.04
##               kurtosis   se
## station*            NA   NA
## nday             -1.37 0.38
## habitat_type*       NA   NA
## daily_sum         9.36 0.16
## ------------------------------------------------------------ 
## group: woodland
##               vars  n mean   sd median trimmed  mad min  max range skew
## station*         1 28  NaN   NA     NA     NaN   NA Inf -Inf  -Inf   NA
## nday             2 28 4.00 2.04      4     4.0 2.97   1    7     6 0.00
## habitat_type*    3 28  NaN   NA     NA     NaN   NA Inf -Inf  -Inf   NA
## daily_sum        4 28 7.29 6.87      5     6.5 4.45   0   24    24 1.36
##               kurtosis   se
## station*            NA   NA
## nday             -1.37 0.38
## habitat_type*       NA   NA
## daily_sum         1.01 1.30
## ------------------------------------------------------------ 
## group: woodland edge
##               vars  n mean   sd median trimmed  mad min  max range skew
## station*         1 28  NaN   NA     NA     NaN   NA Inf -Inf  -Inf   NA
## nday             2 28    4 2.04    4.0    4.00 2.97   1    7     6 0.00
## habitat_type*    3 28  NaN   NA     NA     NaN   NA Inf -Inf  -Inf   NA
## daily_sum        4 28    7 8.54    3.5    5.54 4.45   0   40    40 2.23
##               kurtosis   se
## station*            NA   NA
## nday             -1.37 0.38
## habitat_type*       NA   NA
## daily_sum         5.59 1.61
#It appears as if there is considerable differences between the daily catch of male Aedes albopictus in House - woodland and woodland edge habitats.
#Statistical investigations using generalised linear mixed models (GLMM) will need to occur.  

GLMM, Poisson model “habitat_01”

#Poisson 
alboglmm.p_t3 <- glmer(daily_sum ~ habitat_type + (1|nday) + (1|station),family = poisson, data=t3)
summary(alboglmm.p_t3)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: daily_sum ~ habitat_type + (1 | nday) + (1 | station)
##    Data: t3
## 
##      AIC      BIC   logLik deviance df.resid 
##    383.7    395.9   -186.9    373.7       79 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2172 -0.7873 -0.3679  0.4199  8.9043 
## 
## Random effects:
##  Groups  Name        Variance Std.Dev.
##  nday    (Intercept) 0.5566   0.7460  
##  station (Intercept) 0.4391   0.6626  
## Number of obs: 84, groups:  nday, 7; station, 4
## 
## Fixed effects:
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                -1.5257     0.5493  -2.777  0.00548 ** 
## habitat_typewoodland        3.1209     0.3387   9.215  < 2e-16 ***
## habitat_typewoodland edge   3.0809     0.3390   9.089  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) hbtt_t
## hbtt_typwdl -0.590       
## hbtt_typwde -0.590  0.957
Anova(alboglmm.p_t3)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: daily_sum
##               Chisq Df Pr(>Chisq)    
## habitat_type 85.795  2  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Test for overdispersion
overdis.P_t3<-gof(alboglmm.p_t3)
overdis.P_t3
## D  = 138.7412, df = 79, P(>D) = 3.816759e-05
## X2 = 197.7626, df = 79, P(>X2) = 3.677099e-12
sum(residuals(alboglmm.p_t3,"pearson")^2)
## [1] 197.7626

GLMM, negative binomial model “habitat_01”

alboglmm.nb_t3 <- glmer.nb(daily_sum ~ habitat_type + (1|station) + (1|nday), data=t3)
summary(alboglmm.nb_t3)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Negative Binomial(4.933)  ( log )
## Formula: daily_sum ~ habitat_type + (1 | station) + (1 | nday)
##    Data: t3
## 
##      AIC      BIC   logLik deviance df.resid 
##    356.8    371.3   -172.4    344.8       78 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.3242 -0.5881 -0.3500  0.3510  7.8204 
## 
## Random effects:
##  Groups  Name        Variance Std.Dev.
##  nday    (Intercept) 0.5321   0.7294  
##  station (Intercept) 0.4320   0.6573  
## Number of obs: 84, groups:  nday, 7; station, 4
## 
## Fixed effects:
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                -1.4984     0.5578  -2.686  0.00722 ** 
## habitat_typewoodland        3.1238     0.3706   8.428  < 2e-16 ***
## habitat_typewoodland edge   3.0260     0.3706   8.164 3.24e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) hbtt_t
## hbtt_typwdl -0.601       
## hbtt_typwde -0.599  0.894
Anova(alboglmm.nb_t3)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: daily_sum
##               Chisq Df Pr(>Chisq)    
## habitat_type 73.005  2  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#overdispersion
overdis.NB_t3 <- gof(alboglmm.nb_t3)
overdis.NB_t3
## D  = 74.502, df = 78, P(>D) = 0.5912466
## X2 = 115.1485, df = 78, P(>X2) = 0.003991458
sum(residuals(alboglmm.nb_t3,"pearson")^2)
## [1] 115.1485

Comparing models

#comparing overdispersion
sum(residuals(alboglmm.p_t3,"pearson")^2)
## [1] 197.7626
sum(residuals(alboglmm.nb_t3,"pearson")^2)
## [1] 115.1485
#alboglmm.nb_t3 (115) is less overdispered than alboglmm.p_t3 (197).

#comparing the AIC weights between models
AIC(alboglmm.p_t3 ,alboglmm.nb_t3)
##                df      AIC
## alboglmm.p_t3   5 383.7172
## alboglmm.nb_t3  6 356.7513
#alboglmm.nb_t3 has a lower AIC (356) than alboglmm.p_t3 (383) and is less overdispersed. 

Tukey post hoc test for significance - negative binomial model

tukeyTrap_t3<-emmeans(alboglmm.nb_t3, list(pairwise ~ habitat_type), adjust = "tukey")
tukeyTrap_t3
## $`emmeans of habitat_type`
##  habitat_type  emmean    SE  df asymp.LCL asymp.UCL
##  house          -1.50 0.558 Inf    -2.830    -0.167
##  woodland        1.63 0.447 Inf     0.558     2.693
##  woodland edge   1.53 0.448 Inf     0.457     2.598
## 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of habitat_type`
##  contrast                 estimate    SE  df z.ratio p.value
##  house - woodland          -3.1238 0.371 Inf -8.428  <.0001 
##  house - woodland edge     -3.0260 0.371 Inf -8.164  <.0001 
##  woodland - woodland edge   0.0979 0.171 Inf  0.574  0.8342 
## 
## Results are given on the log (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 3 estimates

Stations at house habitats captured significantly less male Ae. albopictus compared to stations set at both woodland edge and woodland habitats. There was no significant difference between average daily catch of male Ae. albopictus between woodland edge and woodland stations.

This result is displayed in Figure 5 of the Swan et al. 2020 manuscript (graph produced in GraphPad Prism 8).