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)]),]
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.")
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
#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
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"))
tab_model(topmod, transform = NULL, wrap.labels = 50, title = "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 |
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
# 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
# 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
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
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