“The impact of sound lure frequency and habitat type on male Aedes albopictus capture rates with the Male Aedes Sound Trap”
Preliminaries:
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)
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
#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.
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
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
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.
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).
#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.
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
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 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.
#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).
#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.
#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
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 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.
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).