2. Punctual scale model
2.1 Input data and model settings
At the punctual scale, the model only requires a temperature time series, recorded by e.g. a weather station, provided as a numerical matrix (in degree Celsius). For the purpose of this tutorial, we simulate a 1-year long temperature time series and, for the sake of brevity, we will not discuss the chunks of code already presented in other tutorials.
# [1] "en_GB.UTF-8"ndays <- 365*1 #length of the time series in days
set.seed(123)
sim_temp <- create_sims(n_reps = 1,
n = ndays,
central = 16,
sd = 2,
exposure_type = "continuous",
exposure_trend = "cos1", exposure_amp = -1.0,
average_outcome = 12,
outcome_trend = "cos1",
outcome_amp = 0.8,
rr = 1.0055)
# Model settings
## Define the day of introduction (July 1st is day 1)
str <- "2000-07-01"
## Define the end-day of life cycle (August 1st is the last day)
endr <- "2000-08-01"
## Define the number of eggs to be introduced
ie <- 1000
## Define the number of model iterations
it <- 1 # The higher the number of simulations the better
## Define the number of liters for the larval density-dependent mortality
habitat_liters <- 1
## Define latitude and longitude for the diapause process
myLat <- 42
myLon <- 7
## Define the number of parallel processes (for sequential iterations set nc=1)
cl <- 1
## convert float temperatures to integer
df_temp <- data.frame("Date" = sim_temp[[1]]$date, "temp" = sim_temp[[1]]$x)
w <- t(as.integer(df_temp$temp*1000)[format(as.Date(str)+1,"%j"):format(as.Date(endr)+1,"%j")])2.2 Run the model
It is crucial to run the model the model specifying the argument compressed.output = FALSE. This will return the number of simulated individuals for each sub-compartments.
simout <- dynamAedes.m(species="albopictus",
scale="ws",
jhwv=habitat_liters,
temps.matrix=w,
startd=str,
endd=endr,
n.clusters=cl,
iter=it,
intro.eggs=ie,
compressed.output=FALSE,
lat=myLat,
long=myLon,
verbose=FALSE,
seeding=TRUE)2.3 Analyse the results
A first summary of simulations can be obtained with:
summary(simout)# Summary of dynamAedes simulations:
# ----------------------------------
# Species: Aedes albopictus
# Scale: WEATHER STATION
# Start Date: 2000-07-01
# End Date: 2000-08-01
# Number of Iterations: 1
# Introduced Stage: egg
# Number Introduced: 1000
# Is Output Compressed?: No
# Water in the System: 1 L
# Min days with population: 30
# Max days with population: 30The simout object is a S4 object where the outputs of the model and related details are saved in different slot. For example, the number of model iterations is saved in:
simout@n_iterations# [1] 1The model output, i.e. the number of simulated individuals, is stored in simout@simulation. For the uncompressed model, simout@simulation is a list where the first level stores the simulation of different iteration, while the second corresponds to the simulated days in the corresponding iteration. If we inspect the first iteration, we observe that the model has computed length(simout[[1]]) days, since we have started the simulation on the 1st of July and ended on the 1st of August.
length(simout@simulation[[1]])# [1] 30The third level corresponds to an array reporting, for a given iteration and a given day, the number of individuals belonging to each compartment (rows) for each sub-compartment (the third dimension of the array, noted as sc1-…-scN in the print). As example, if we inspect the 1st and the 15th day within the first iteration, we obtain a matrix having:
class(simout@simulation[[1]][[1]])# [1] "array"simout@simulation[[1]][[1]]# , , sc1
#
# [,1]
# egg 0
# juvenile 355
# adult 0
# diapause_egg 0
#
# , , sc2
#
# [,1]
# egg 0
# juvenile 0
# adult 0
# diapause_egg 0
#
# , , sc3
#
# [,1]
# egg 0
# juvenile 0
# adult 0
# diapause_egg 0
#
# , , sc4
#
# [,1]
# egg 645
# juvenile 0
# adult 0
# diapause_egg 0
#
# , , sc5
#
# [,1]
# egg 0
# juvenile 0
# adult 0
# diapause_egg 0
#
# , , sc6
#
# [,1]
# egg 0
# juvenile 0
# adult 0
# diapause_egg 0simout@simulation[[1]][[15]]# , , sc1
#
# [,1]
# egg 0
# juvenile 82
# adult 5
# diapause_egg 0
#
# , , sc2
#
# [,1]
# egg 48
# juvenile 17
# adult 2
# diapause_egg 0
#
# , , sc3
#
# [,1]
# egg 0
# juvenile 5
# adult 0
# diapause_egg 0
#
# , , sc4
#
# [,1]
# egg 80
# juvenile 1
# adult 1
# diapause_egg 0
#
# , , sc5
#
# [,1]
# egg 0
# juvenile 2
# adult 0
# diapause_egg 0
#
# , , sc6
#
# [,1]
# egg 0
# juvenile 1
# adult 0
# diapause_egg 02.4 Derive abundance 95% CI for a specific sub-compartment
We can use the auxiliary functions of the package to analyse the results. We now compute the interquantile range abundance for the host-seeking sub-compartment of the simulated population using the function adci.
# Retrieve the maximum number of simulated days
dd <- max(simout)
# Compute the inter-quartile of abundances along the iterations
breaks <- c(0.25,0.50,0.75)
ed <- 1:dd
hs <- adci(simout, eval_date=ed, breaks=breaks,
stage="Adults",
sub_stage = "Host-seeking" )
head(hs)# X25. X50. X75. day stage sub_stage
# 1 0 0 0 1 Adult Host-seeking
# 2 0 0 0 2 Adult Host-seeking
# 3 0 0 0 3 Adult Host-seeking
# 4 0 0 0 4 Adult Host-seeking
# 5 0 0 0 5 Adult Host-seeking
# 6 0 0 0 6 Adult Host-seekingtail(hs)# X25. X50. X75. day stage sub_stage
# 25 1 1 1 25 Adult Host-seeking
# 26 1 1 1 26 Adult Host-seeking
# 27 4 4 4 27 Adult Host-seeking
# 28 1 1 1 28 Adult Host-seeking
# 29 1 1 1 29 Adult Host-seeking
# 30 1 1 1 30 Adult Host-seekingWe can now simply plot it.
ggplot(hs, aes(x=day, y=X50., group=factor(stage), col=factor(stage))) +
ggtitle("Host-seeking Ae. albopictus Interquantile range abundance") +
geom_ribbon(aes(ymin=X25., ymax=X75., fill=factor(stage)),
col="white",
alpha=0.2,
outline.type="full") +
geom_line(linewidth=0.8) +
ylim(0,10)+
labs(x="Date", y="Interquantile range abundance", col="Stage", fill="Stage") +
theme_classic() +
theme(legend.position="bottom",
text = element_text(size=16),
strip.text = element_text(face = "italic"))
