Requires files: dat7.csv, lc.csv, Diversity_Analysis.r

library(readxl)
library(reshape2)
library(glmmTMB)
library(bbmle)
library(ggplot2)
library(plyr)
library(dplyr)
library(ggpubr)
library(kableExtra)
library(performance)
library(sjPlot)
library(runjags)
library(coda)
library(ggmcmc)
library(DHARMa)
dat7 <- read.csv("dat7.csv")
spp_counts <- aggregate(ct~spp, dat7, sum)
spp_occ <- aggregate(I(ct>0)~spp, dat7, mean)
dat8 <- dat7[which(dat7$spp%in%spp_counts$spp[which(spp_counts$ct>160)]),] 

Diversity analysis

Our hierarchical Bayesian multi-species occupancy model can be expressed as,

\[ R = n + \sum_{i=n+1}^{n+n_{aug}}{w_{(i)}}\\ w_{i} \sim \mathrm{Bernoulli}(\Omega)\\ Z_{j,i} \sim \mathrm{Bernoulli}(\psi_{j,i}\times w_{i})\\ X_{j,k,i} \sim \mathrm{Bernoulli}(p_{j,k,i} \times Z_{j,i}) \],

Where \(R\) is species richness, \(n\) is the number of of observed species, \(n_{aug}\) is the number of augmented (all-zero capture history) species added to the dataset, \(Z\) is the latent occupancy variable, and \(X\) is the data. Indexes refer to site (\(i\)), occasion (\(k\)), and species (\(i\)). The occurrence (\(\psi\)) and observation (\(p\)) processes were modeled as hierarchical random effects,

\[ \mathrm{logit}(p_{j,k,i}) = v_{i} \\ v_{i} \sim \mathrm{Normal}(\mu_v, \tau_v) \\ \mathrm{logit}(\psi_{j,i}) = u_{i} \\ u_{i} \sim \mathrm{Normal}(\mu_u, \tau_u) \].

We used weakly informative Gaussian priors for mean of \(v_{i}\) and \(u_{i}\) (mean = 0, sd = 2.25), and vague gamma priors for their precision (r = .1, lambda = 0.1). Omega was given a uniform (0, 1) prior. We also monitored the following derived diversity measures: \[ r_j = \sum_{i=1}^{n+n_{aug}}{Z_{j,i}} \\ \alpha = \mathrm{Mean}(r_j) \\ \beta = \frac{R}{\alpha}\\ \zeta = \sum_{i=1}^{n+n_{aug}}\prod_{j=1}^{J}Z_{j,i} \],

Where \(r_j\) is “local” species richness at site \(j\), \(\alpha\) is alpha-diversity, \(\beta\) is beta-diversity, and \(\zeta\) is zeta-diversity. We also monitored the number of sites each species occupied, \(Nsites_i = \sum_{j=1}^{J}{Z_{j,i}}\).

if (file.exists("runjagsfiles")){
   divout <- results.jags("runjagsfiles")
} else {
source("Diversity_Analysis.r")
   divout <- diversity_analysis()
}
## 
## Attaching package: 'reshape'
## The following objects are masked from 'package:tidyr':
## 
##     expand, smiths
## The following object is masked from 'package:dplyr':
## 
##     rename
## The following objects are masked from 'package:plyr':
## 
##     rename, round_any
## The following objects are masked from 'package:reshape2':
## 
##     colsplit, melt, recast
## Calling 3 simulations using the parallel method...
## Following the progress of chain 1 (the program will wait for all chains
## to finish before continuing):
## Welcome to JAGS 4.3.0 on Sat Sep 05 22:23:25 2020
## JAGS is free software and comes with ABSOLUTELY NO WARRANTY
## Loading module: basemod: ok
## Loading module: bugs: ok
## . . Reading data file data.txt
## . Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 7840
##    Unobserved stochastic nodes: 665
##    Total graph size: 9336
## . Reading parameter file inits1.txt
## . Initializing model
## . Adapting 1000
## -------------------------------------------------| 1000
## ++++++++++++++++++++++++++++++++++++++++++++++++++ 100%
## Adaptation successful
## . Updating 10000
## -------------------------------------------------| 10000
## ************************************************** 100%
## . . . . . . . . . . . . . . . . . Updating 20000
## -------------------------------------------------| 20000
## ************************************************** 100%
## . . . . Updating 0
## . Deleting model
## . 
## All chains have finished
## Simulation complete.  Reading coda files...
## Coda files loaded successfully
## Note: Summary statistics were not produced as there are >50 monitored
## variables
## [To override this behaviour see ?add.summary and ?runjags.options]
## FALSEFinished running the simulation
## JAGS files were saved to the 'runjagsfiles' folder in your current
## working directory
res.mcmc.list <- as.mcmc.list(divout)
res.ggs <- ggs(res.mcmc.list)
res.coda <- summary(divout)
## Calculating summary statistics...
## Calculating the Gelman-Rubin statistic for 102 variables....
res.df <- data.frame(res.coda)
library(xtable)
msomtab <- res.df[grep("Nsites\\[1\\]|Nsites\\[7\\]|alpha|beta|zeta|R|mu|sigma|omega",rownames(res.df)),c("Lower95", "Median", "Upper95", "psrf")]
msomtab$Parameter  <- rownames(msomtab)
msomtab<-msomtab[c(3,1,11,2,12,4,7,8,9,5,6),c("Parameter", "Median", "Lower95", "Upper95")]
kable(msomtab, format='html', digits=3, row.names=FALSE, escape=FALSE, caption="Table 3.")
Table 3.
Parameter Median Lower95 Upper95
omega 0.630 0.359 0.940
mu.u -0.737 -2.535 0.916
sigma.u 2.059 0.431 4.364
mu.v -2.265 -3.675 -1.149
sigma.v 1.625 0.339 3.481
R 12.000 10.000 18.000
alpha 4.433 3.000 7.233
beta 2.806 1.423 4.688
zeta 0.000 0.000 1.000
Nsites[1] 26.000 24.000 30.000
Nsites[7] 22.000 21.000 24.000
#land cover data
data1. <- dat7 %>%  arrange(td) %>% 
    mutate(Species =  sub("female", "", spp), Point = site, 
        Rep = as.numeric(as.factor(td)))
total.count. <- tapply(data1.$ct, data1.$Species, sum)
data1. <- subset(data1., Species %in% names(total.count.)[total.count.>0]) %>% 
    arrange(Point, Species)
uspecies <- as.character(unique(data1.$Species))
lc <- data1. %>% 
    group_by(Point) %>% 
    summarize(lc.ag=first(Agricultural), lc.ma=first(Mangrove),
        lc.rf=first(Rainforest), lc.sf=first(Scrub),
        lc.ur=first(Urban))
## `summarise()` ungrouping output (override with `.groups` argument)
rtab <- res.df[grep("r", rownames(res.df)), c("Lower95", "Median", "Upper95")] %>%
    mutate(param=rownames(.)) %>% 
    mutate(Point=as.numeric(substr(param,3,nchar(param)-1))) %>% 
    left_join(lc)
## Joining, by = "Point"
lcsiteplot <- function(lc.var="lc.sf"){
    rtab %>% ggplot(aes(x=eval(parse(text=lc.var)), y=Median)) + 
        geom_pointrange(aes(ymin=Lower95, ymax=Upper95))+ theme_classic() + xlim(0,1) + 
        scale_y_continuous(breaks = c(2,4,6,8,10))
}
#plot panels
dlpanel1 <- lcsiteplot("lc.sf") + ylab("N of Species") + xlab("Scrub")
dlpanel2 <- lcsiteplot("lc.ag") + ylab("") + xlab("Agriculture")
dlpanel3 <- lcsiteplot("lc.ma") + ylab("") + xlab("Mangrove")
dlpanel4 <- lcsiteplot("lc.rf") + ylab("N of Species") + xlab("Rainforest")
dlpanel5 <- lcsiteplot("lc.ur") + ylab("") + xlab("Urban")
dlpanel6 <- msomtab[grep("R|alpha|beta|zeta", rownames(msomtab)), 
        c("Lower95", "Median", "Upper95")] %>%
    mutate(param=rownames(.)) %>% filter(param!='R_exp') %>% 
    mutate(param=factor(param, levels=c("R", "alpha", "beta", "zeta"))) %>%
    arrange(-Median) %>%
    ggplot(aes(x=param, y=Median,)) + 
    geom_pointrange(aes(ymin=Lower95, ymax=Upper95)) + 
    ylab("") + xlab("") + theme_classic() + 
    scale_x_discrete(labels=c(parse(text="italic('R')"), expression(alpha), 
        expression(beta), expression(zeta)))
#plot
lcplotlist<- list(dlpanel1, dlpanel2, dlpanel3, dlpanel4, dlpanel5, dlpanel6)
LCDIVERSITY_PLOT <- ggarrange(plotlist=lcplotlist, labels=c("a)","b)","c)","d)","e)","f)"))
ggsave("Figure2.pdf", plot = LCDIVERSITY_PLOT, device = "pdf", width = 6.5, height = 4, dpi = 600)
LCDIVERSITY_PLOT
Figure 2.

Figure 2.

Relative abundance analysis

#Re-scale covariates prior to fitting the model
allcols <- ncol(dat8)
covlist <- c("precip_in", "Agricultural", "Rainforest", "Urban", "Scrub", "Mangrove")
for(i in 1:length(covlist)){
  dat8[,(allcols+i)] <- scale(dat8[,covlist[i]])
  names(dat8)[(allcols+i)] <- paste("c.",covlist[i],sep="")
}
dat8$pos <- numFactor(dat8$lon, dat8$lat) # for use in spatial autocorrelation test
dat8$group <- factor(rep(1, nrow(dat8)))  # for use in spatial autocorrelation test
# Mangrove mosquitoes and mangrove site location
dat8$mtrait <- ifelse(dat8$spp %in% c("De.ma.female","Ae.ta.female"), 1, 0)
dat8$msite <- as.numeric(dat8$landuse=="Mangrove")
dat8$c.Anthropogenic <- with(dat8, scale(Agricultural + Urban))
# Model hypotheses (suppress convergence warnings in nb1 and poisson models)
hlist <- list(
    H1 = ct ~ c.precip_in + spp:(c.Agricultural + c.Rainforest + c.Urban) + 
        spp*Mangrove + (1|site) + (c.precip_in|spp),
    H2 = ct ~ c.precip_in + spp:(c.Anthropogenic) + (1|site) + (c.precip_in|spp),
    H3 = ct ~ c.precip_in + spp:(c.Agricultural + c.Urban) + (1|site) + 
        (c.precip_in|spp),
    H4 = ct ~ c.precip_in + spp:(c.Anthropogenic) + mtrait*msite + (1|site) + 
        (c.precip_in|spp),
    H5 = ct ~ c.precip_in + spp:(c.Agricultural + c.Urban) + mtrait*msite + 
        (1|site) + (c.precip_in|spp),
    H6 = ct ~ c.precip_in + spp:(c.Anthropogenic + c.Scrub) + mtrait*msite + 
        (1|site) + (c.precip_in|spp),
    H7 = ct ~ c.precip_in + spp:(c.Agricultural + c.Urban + c.Scrub) + 
        mtrait*msite + (1|site) + (c.precip_in|spp),
    H8 = ct ~ c.precip_in + spp:(c.Anthropogenic + c.Rainforest) + mtrait*msite + 
        (1|site) + (c.precip_in|spp),
    H9 = ct ~ c.precip_in + spp:(c.Agricultural + c.Urban + c.Rainforest) + 
        mtrait:msite + (1|site) + (c.precip_in|spp),
    H10 = ct ~ c.precip_in + spp:(c.Anthropogenic + c.Scrub  + c.Rainforest) + 
        mtrait*msite + (1|site) + (c.precip_in|spp),
    H11 = ct ~ c.precip_in + spp:(c.Agricultural + c.Urban + c.Scrub + 
        c.Rainforest) + mtrait*msite + (1|site) + (c.precip_in|spp),
    H12 = ct ~ c.precip_in + spp:(c.Agricultural + c.Urban + c.Rainforest) + 
        mtrait*msite + (1|site) + (c.precip_in|spp)
)
# fit models
if (file.exists("models.RData")){
    load("models.RData")
} else {
    models <- lapply(hlist, function(x) glmmTMB(x, data=dat8, family=nbinom2))
    # re-fit with alternative error distributions - many throw convergence warnings
    models.nb1 <- lapply(models, function(x) update(x, .~., family=nbinom1))
    models.poi <- lapply(models, function(x) update(x, .~., family=poisson))
    #Save model list objects
save(models, models.nb1, models.poi, file="models.RData")
}
if (file.exists("compare_AIC.RData")){
    load("compare_AIC.RData")
} else {
    #!!!MUST!!! set 'sort = FALSE' for selection of top
    variable_selection <- AICtab(models, sort=FALSE, base = T)
    R.models <- sapply(models, r2)
    topmod <- models[[which(variable_selection$dAIC==0)]]
    ### Other families
    varsel_pois <- AICtab(models.poi, sort=FALSE, base = T)
    R.poi <- sapply(models.poi, r2)
    varsel_nb1 <- AICtab(models.nb1, sort=FALSE, base = T)
    R.nb1 <- sapply(models.nb1, r2)
    save(variable_selection, R.models, topmod, varsel_pois, varsel_nb1, 
        R.poi, R.nb1, file = "compare_AIC.RData")
}
vif_list <- lapply(models, check_collinearity)
check_collinearity(topmod)
## # Check for Multicollinearity
## 
## Low Correlation
## 
##           Parameter  VIF Increased SE
##         c.precip_in 1.01         1.00
##              mtrait 1.51         1.23
##               msite 2.53         1.59
##  spp:c.Agricultural 2.76         1.66
##         spp:c.Urban 2.04         1.43
##    spp:c.Rainforest 1.92         1.39
##        mtrait:msite 3.63         1.91
res <- simulateResiduals(topmod)
plot(res, quantreg = T)

res_site <- recalculateResiduals(res, group = dat8$site)
testSpatialAutocorrelation(res_site, x =  aggregate(dat8$lon, list(dat8$site), 
        mean)$x, y = aggregate(dat8$lat, list(dat8$site), mean)$x, plot = TRUE)

## 
##  DHARMa Moran's I test for spatial autocorrelation
## 
## data:  res_site
## observed = -0.036577, expected = -0.034483, sd = 0.066406, p-value =
## 0.9748
## alternative hypothesis: Spatial autocorrelation
res2 = recalculateResiduals(res, group = dat8$td)
testTemporalAutocorrelation(res2, time = as.Date(unique(dat8$td), format="%m/%d/%Y"))

## 
##  Durbin-Watson test
## 
## data:  simulationOutput$scaledResiduals ~ 1
## DW = 1.8622, p-value = 0.7911
## alternative hypothesis: true autocorrelation is not 0

Results

spp_summary <- left_join(spp_counts, spp_occ)
## Joining, by = "spp"
names(spp_summary)[3] <- "occ"

The best count model was H11, which can be expressed as,
\[ \begin{aligned} y_{(i,j,s)} &\sim \mathrm{NegBin}(\mu_{(i,j,s)}, \theta) \\ \mathrm{log}({\mu_{(i,j,s)}}) &= \beta_{0} + \alpha_{0(s)} + (\beta_{1} + \alpha_{1(s)})Precip_{(j)} + \beta_{2}m_{trait(s)} + \beta_{3} Mangrove_{(i)} + \beta_{4(s)}LocalAgriculture_{(i)} + \\ &\beta_{5(s)}LocalUrban_{(i)} + \beta_{6(s)}LocalRainforest_{n(i)} + \beta_{7}m_{trait(s)}Mangrove_{(i)} + \alpha_{2(i)} \\ \mathrm{log}(\theta) &= \eta_{(s)}, \end{aligned} \] where \(i\) indexes site, \(j\) indexes sampling occasion, and \(s\) indexes species. There are random intercepts for each species(\(\beta_{0} + \alpha_{0(s)}\)), a main effect of precipitation with species-specific random slopes (\((\beta_{1} + \alpha_{1(s)})Precip_{(j)}\)), a main effect of the Magrove breeding trait (\(\beta_{2(m_{trait(s)})}\)) multiplied by sites where points are located in Mangroves (\(Mangrove_{(i)}\)), species-specific main effects on proportional land cover variables (\(\beta_{3(s)}LocalAgriculture_{(i)} + \beta_{4(s)}LocalUrban_{(i)} + \beta_{5(s)}Local2^\circ Forest_{(i)} + \beta_{6(s)}LocalRainforest_{n(i)}\)), a site-level random intercept term (\(\alpha_{2(i)}\)), and a dispersion parameter, \(\eta_{(s)}\). The top model employs a negative binomial parameterization where the variance scales quadratically with the mean, \(Var(y_{(i,j,s)})=\mu_{(i,j,s)}(1+\frac{\mu_{(i,j,s)}}{\phi})\), where the predicted mean is \(\mu\), and \(\phi=e^{\eta}\) (Hardin & Hilbe 2007).

newdata = dat8
newdata$site <- NA
temp = predict(topmod, data = newdata, se.fit=TRUE, type="link")
newdata$predlink = temp$fit
newdata$predFE = exp(temp$fit)
newdata$predFE.min = exp(temp$fit-1.98*temp$se.fit)
newdata$predFE.max = exp(temp$fit+1.98*temp$se.fit)
real=ddply(dat8, ~spp+landuse+td, summarize, m=mean(ct), min = min(ct), max = max(ct))
spp_names <- function(string) {
    string <- gsub("\\.", "\\. ", string) #replace period with period & space
    tag <- sub('.*\\.', '.', string) #define tag
    string <- gsub(tag, ".", string) #remove tag at end of species name
    string <- paste("italic('",string,"')", sep="")
}
real <- transform(real, spp = factor(spp, levels=unique(real$spp), 
    labels=sapply(unique(real$spp), function(i) spp_names(i))))
real$landuse <- factor(real$landuse, 
    levels=c("Scrub", "Agricultural", "Mangrove", "Rainforest", "Urban"))
newdata2 <- transform(newdata, spp = factor(spp, levels=unique(newdata$spp), 
    labels=sapply(unique(newdata$spp), function(i) spp_names(i))))
newdata2 <- ddply(newdata2, ~spp+landuse+td, summarize,m=mean(predFE), 
    lo=exp(mean(predlink)-1.98*sd(predlink)), 
    hi=exp(mean(predlink)+1.98*sd(predlink)))
newdata2$landuse <- factor(newdata2$landuse, 
    levels = c("Scrub", "Agricultural", "Mangrove", "Rainforest", "Urban"))

Tables

tab_model(topmod, transform = NULL, wrap.labels = 50, title = "Table 4.")
Table 4.
  ct
Predictors Log-Mean CI p
(Intercept) -0.75 -1.33 – -0.18 0.010
c.precip_in 0.66 0.13 – 1.19 0.015
mtrait -1.40 -2.12 – -0.69 <0.001
msite 0.01 -1.29 – 1.31 0.990
spp [Ae.ae.female] : c.Agricultural 0.20 -0.38 – 0.77 0.503
spp [Ae.ta.female] : c.Agricultural -1.35 -1.98 – -0.71 <0.001
spp [Cu.qu.female] : c.Agricultural -1.07 -1.65 – -0.49 <0.001
spp [De.ma.female] : c.Agricultural -1.08 -1.80 – -0.35 0.003
spp [Ae.ae.female] : c.Urban 0.13 -0.42 – 0.68 0.644
spp [Ae.ta.female] : c.Urban -2.11 -2.91 – -1.30 <0.001
spp [Cu.qu.female] : c.Urban 0.32 -0.22 – 0.85 0.248
spp [De.ma.female] : c.Urban -0.97 -1.67 – -0.26 0.007
spp [Ae.ae.female] : c.Rainforest -0.43 -0.99 – 0.13 0.129
spp [Ae.ta.female] : c.Rainforest -1.21 -1.94 – -0.48 0.001
spp [Cu.qu.female] : c.Rainforest -1.39 -2.03 – -0.76 <0.001
spp [De.ma.female] : c.Rainforest -0.91 -1.64 – -0.18 0.014
mtrait * msite 2.32 0.83 – 3.81 0.002
Random Effects
σ2 2.60
τ00 site 0.54
τ00 spp 0.02
τ11 spp.c.precip_in 0.24
ρ01 spp -0.69
ICC 0.24
N site 30
N spp 4
Observations 1568
Marginal R2 / Conditional R2 0.601 / 0.695

Figures

predPalette <- c("#E69F00", "#F0E442", "#D55E00", "#009E73", "#888888")
gg_color <- function(n) {
  hues = seq(15, 375, length=n+1)
  hcl(h = hues, l = 65, c = 100)[1:n]
}
predPalette <- gg_color(5)
predcovers <- c("Scrub", "Agricultural", "Mangrove", "Rainforest", "Urban")
newdata2$name <- newdata2$spp
newdata2$name <- sub("Ae. ae.", "Aedes aegypti", newdata2$name)
newdata2$name <- sub("Ae. ta.", "Aedes taeniorhynchus", newdata2$name)
newdata2$name <- sub("Cu. qu.", "Culex quinquefasciatus", newdata2$name)
newdata2$name <- sub("De. ma.", "Deinocerites magnus", newdata2$name)
plotpreds <- function(name, ymin, ymax){newdata2[grep(name, newdata2$name),] %>%
    ggplot(aes(as.Date(td, "%m/%d/%Y"), m, colour=landuse)) + 
    ggtitle(parse(text=unique(grep(name,newdata2$name, value=T)))) + 
    geom_line(size = 1) + 
    theme_classic() + 
    theme(plot.title=element_text(hjust=0.5, size=11), 
        axis.text.x=element_text(angle=30, hjust=1)) + 
    ylab("") + xlab("") + 
  scale_x_date(date_labels = "%b %Y", 
    breaks=pretty(as.Date(newdata2$td, "%m/%d/%Y"), 5)) + 
  ylim(ymin,ymax)
}
panela <- plotpreds("Aedes aegypti", 0, 5)
panelb <- plotpreds("Aedes taeniorhynchus", 0, 200)
panelc <- plotpreds("Culex quinquefasciatus", 0, 15)
paneld <- plotpreds("Deinocerites magnus", 0, 15)

PREDICTIONS <- ggarrange(panela, panelb, panelc, paneld, ncol=2, nrow=2, 
        labels=c("a)","b)","c)","d)"), common.legend=TRUE, legend="right")

ggsave("Figure3.pdf", plot = PREDICTIONS, device = "pdf", width = 7.5, height = 6, dpi = 600)
PREDICTIONS
Figure 3.

Figure 3.

# Plot BLUPs
library(ggeffects)
precippred <- ggpredict(topmod, c("c.precip_in"), type = "fe")
blups <- ggpredict(topmod, c("c.precip_in", "spp"), type = "re")
pr <- ggplot(blups, aes(x = x , y = predicted, color = group)) + 
    geom_line() + 
    xlab("Precipitation (cm)") + ylab("Count") + theme_classic() + 
    theme(legend.title=element_blank(), 
        legend.justification=c(.05,.95), 
        legend.position=c(.05,.95)) + 
    scale_color_discrete(
        labels = c("Aedes aegypti", "Aedes taeniorhynchus", 
             "Culex quinquefasciatus", "Deinocerites magnus"),
        guide = guide_legend(label.theme = element_text(size = 11, angle = 0, face = "italic")))
ggsave("Figure4.pdf", plot = pr, device = "pdf", width = 3.5, height = 3.5, dpi = 600)
pr
Figure 4.

Figure 4.

Appendixes

# plot
cbPalette <- c("#F0E442", "#D55E00", "#0072B2", "#56B4E9", "#009E73", "#E69F00", "#888888")
lc <- read.csv("lc.csv")
LANDCOVER_BARPLOT <- ggplot() + 
    geom_bar(aes(y = pct, x = site, fill = LandCover), data = lc, stat="identity") +  
    scale_fill_manual(values=cbPalette, name = "Land Cover Class") + 
    labs(fill = "Land Cover Class") + ylab("Proportion of Land Cover") + 
    xlab("Site") + theme_classic() + scale_x_continuous(breaks = seq(0,30,by=2))
ggsave("FigureS1.pdf", plot = LANDCOVER_BARPLOT, device="pdf", width=6.5, height=4, dpi=600)
LANDCOVER_BARPLOT
Figure S1.

Figure S1.

real$name <- real$spp
real$name <- sub("Ae. ae.","Aedes aegypti",real$name)
real$name <- sub("Ae. ta.","Aedes taeniorhynchus",real$name)
real$name <- sub("Cu. qu.","Culex quinquefasciatus",real$name)
real$name <- sub("De. ma.","Deinocerites magnus",real$name)
FIGURE_PREDICTIONS2 <- newdata2 %>% mutate(td=as.Date(td,"%m/%d/%Y")) %>% 
    ggplot(aes(td, m, colour=landuse)) + geom_line() + 
  geom_line(aes(td, lo), lty=2) + 
  geom_line(aes(td, hi), lty=2) + 
  geom_point(data = real, aes(as.Date(td,"%m/%d/%Y"), m), size = 1) +
    facet_grid(name~landuse, scales = "free", labeller = label_parsed) +
    theme(legend.position = "none") + theme_classic() +
  ylab("Average mosquito relative abundance") +
  xlab("Month") + scale_x_date(date_labels = "%b %Y") + 
    theme(strip.background = element_blank(), axis.text.x=element_text(angle=30, hjust=1)) 
ggsave("FigureS2.pdf", plot = FIGURE_PREDICTIONS2, device="pdf", width=7.5, height=10, dpi=600)
FIGURE_PREDICTIONS2
Figure S2. a caption goes here.

Figure S2. a caption goes here.

sessionInfo()
## R version 4.0.2 (2020-06-22)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 14393)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_United States.1252 
## [2] LC_CTYPE=English_United States.1252   
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.1252    
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] ggeffects_0.15.1  xtable_1.8-4      reshape_0.8.8     DHARMa_0.3.2.0   
##  [5] ggmcmc_1.5.0      tidyr_1.1.2       coda_0.19-3       runjags_2.0.4-6  
##  [9] sjPlot_2.8.4      performance_0.4.8 kableExtra_1.2.1  ggpubr_0.4.0     
## [13] dplyr_1.0.2       plyr_1.8.6        ggplot2_3.3.2     bbmle_1.0.23.1   
## [17] glmmTMB_1.0.2.1   reshape2_1.4.4    readxl_1.3.1     
## 
## loaded via a namespace (and not attached):
##   [1] minqa_1.2.4         colorspace_1.4-1    ggsignif_0.6.0     
##   [4] ellipsis_0.3.1      rio_0.5.16          sjlabelled_1.1.6   
##   [7] estimability_1.3    parameters_0.8.2    rstudioapi_0.11    
##  [10] farver_2.0.3        mvtnorm_1.1-1       xml2_1.3.2         
##  [13] codetools_0.2-16    splines_4.0.2       doParallel_1.0.15  
##  [16] knitr_1.29          sjmisc_2.8.5        nloptr_1.2.2.2     
##  [19] broom_0.7.0         shiny_1.5.0         effectsize_0.3.2   
##  [22] compiler_4.0.2      httr_1.4.2          sjstats_0.18.0     
##  [25] emmeans_1.5.0       backports_1.1.9     fastmap_1.0.1      
##  [28] Matrix_1.2-18       later_1.1.0.1       htmltools_0.5.0    
##  [31] tools_4.0.2         gtable_0.3.0        glue_1.4.2         
##  [34] Rcpp_1.0.5          carData_3.0-4       cellranger_1.1.0   
##  [37] vctrs_0.3.4         ape_5.4-1           nlme_3.1-149       
##  [40] iterators_1.0.12    lmtest_0.9-37       insight_0.9.1      
##  [43] xfun_0.16           stringr_1.4.0       openxlsx_4.1.5     
##  [46] lme4_1.1-23         rvest_0.3.6         mime_0.9           
##  [49] lifecycle_0.2.0     statmod_1.4.34      rstatix_0.6.0      
##  [52] zoo_1.8-8           MASS_7.3-52         scales_1.1.1       
##  [55] promises_1.1.1      hms_0.5.3           parallel_4.0.2     
##  [58] TMB_1.7.18          RColorBrewer_1.1-2  yaml_2.2.1         
##  [61] curl_4.3            gridExtra_2.3       bdsmatrix_1.3-4    
##  [64] stringi_1.4.6       highr_0.8           bayestestR_0.7.2   
##  [67] gap_1.2.2           foreach_1.5.0       boot_1.3-25        
##  [70] zip_2.1.1           qgam_1.3.2          rlang_0.4.7        
##  [73] pkgconfig_2.0.3     evaluate_0.14       lattice_0.20-41    
##  [76] purrr_0.3.4         labeling_0.3        cowplot_1.0.0      
##  [79] tidyselect_1.1.0    GGally_2.0.0        magrittr_1.5       
##  [82] R6_2.4.1            generics_0.0.2      pillar_1.4.6       
##  [85] haven_2.3.1         foreign_0.8-80      withr_2.2.0        
##  [88] mgcv_1.8-33         abind_1.4-5         tibble_3.0.3       
##  [91] modelr_0.1.8        crayon_1.3.4        car_3.0-9          
##  [94] rmarkdown_2.3       grid_4.0.2          data.table_1.13.0  
##  [97] forcats_0.5.0       digest_0.6.25       webshot_0.5.2      
## [100] httpuv_1.5.4        numDeriv_2016.8-1.1 munsell_0.5.0      
## [103] viridisLite_0.3.0