2. Simulating phenology
Based on the relationship between development and temperature, and on a temperature time series, it is then possible to run simulations of the pest’s phenology.
2.1. Temperature time series
Any temperature time series can be used to simulate pest phenology. In the following simulations, the time step is expressed in days. If the temperature time series has a time step of 1h, a time step of 1/24 should be specified. In this example, we’ll simulate a theoretical temperature time series using a random generation for the normal distribution with mean equal to 25 and standard deviation equal to 2 (90 days with one temperature data every 6 hours ; 90x4 temperature values).
temp <- stats::rnorm(n = 90*4, mean = 25, sd = 2)
par(mar = c(4, 4, 0, 0))
plot(
x = temp, type = "o",
xlab = "Time (6 hours time step)", ylab = "Temperature",
ylim = c(15, 35)
)
2.2. Model selection and computation
Now we can run simulations, selecting our temperature time series and the models of our choice.
2.2.1. Simulation using a single model
# This model was designed for eggs, larvae and pupae, using simple linear
# regression models (following Campbell 1974).
forecastX <- devRateIBMparam(
tempTS = temp, # the temperature time series described above
timeStepTS = 1/4, # we have 4 temperature values per day (every 6 hours)
eq = list("campbell_74", "campbell_74", "campbell_74"), # the model eq for
# the three life stages
myParam = list(
my_model[["model"]]$pupa, # starting with pupae
my_model[["model"]]$egg, # then adults and eggs
my_model[["model"]]$larva # then larvae and the cycle goes on
),
adultLifeStage = 1, # adult life stage after the first model step (pupae
# to adults)
timeLayEggs = 10, # adults longevity fixed to 10 days
numInd = 10, # 10 individuals for the simulation
stocha = 0.05 # a bit of intra-population variability
)The result of the simulation is a list object, with the first element being the predictions of the various stages for all individuals (phenology). The data correspond to the number of time steps expressed in the time series value (6h time step here). The second element corresponds to the model parameters used, and the third element to the temperature series used. We can plot the simualtion result using the first element.
sim01 <- forecastX[[1]]
hist(
-999, # empty histogram
xlim = c(0, max(sim01)), ylim = c(0, nrow(sim01)),
xlab = "Time", ylab = "Number of individuals", axes = FALSE, main = ""
)
axis(1)
axis(2, las = 1)
cycleLabels <- rep(c("adults", "larvae", "pupae"), 10)
trash <- lapply(1:ncol(sim01), function(i){
hist(sim01[,i], add = TRUE)
text(
x = mean(sim01[,i]),
y = max(table(sim01[,i])),
labels = colnames(sim01)[i],
pos = 2
)
text(
x = mean(sim01[,i]),
y = max(table(sim01[,i]))+1,
labels = cycleLabels[i],
pos = 2
)
})
2.2.2. Simulation using different models
set.seed(1234)
my_model01 <- ha_bartekova2006(plotfig = FALSE)
my_model02 <- ha_mironidis2008_nls(plotfig = FALSE)
my_model03 <- ha_foley1981(plotfig = FALSE)
forecastX <- devRateIBMparam(
tempTS = temp, # the temperature time series described above
timeStepTS = 1/4, # we have 4 temperature values per day (every 6 hours)
eq = rep(list(
"campbell_74", "campbell_74", "lactin2_95"
), 3), # the model eq for
# the three life stages
myParam = list(
my_model03[["model"]]$diapausingpupae, # starting with diapausing pupae (Foley 1981)
my_model01[["model"]]$egg, # then adults and eggs (Bartekova 2006)
my_model02[["model"]]$larva, # then larvae (Mironidis 2008)
my_model03[["model"]]$nondiapausingpupae, # pupae G2 (Foley 1981)
my_model01[["model"]]$egg,
my_model02[["model"]]$larva,
my_model03[["model"]]$nondiapausingpupae,
my_model01[["model"]]$egg,
my_model02[["model"]]$larva
),
adultLifeStage = 1, # adult life stage after the first model step (pupae
# to adults)
timeLayEggs = 10, # adults longevity fixed to 10 days
numInd = 10, # 10 individuals for the simulation
stocha = 0.05 # a bit of intra-population variability
)sim02 <- forecastX[[1]]
hist(
-999, # empty histogram
xlim = c(0, max(sim02)), ylim = c(0, nrow(sim02)),
xlab = "Time", ylab = "Number of individuals", axes = FALSE, main = ""
)
axis(1)
axis(2, las = 1)
cycleLabels <- rep(c("adults", "larvae", "pupae"), 10)
trash <- lapply(1:ncol(sim02), function(i){
hist(sim02[,i], add = TRUE)
text(
x = mean(sim02[,i]),
y = max(table(sim02[,i])),
labels = colnames(sim02)[i],
pos = 2
)
text(
x = mean(sim02[,i]),
y = max(table(sim02[,i]))+1,
labels = cycleLabels[i],
pos = 2
)
})
2.3. Simulation results and comparison
## g1s1 g1s2 g1s3 g2s1 g2s2
## Min. 52.00 87.00 187.0 241.0 276.0
## 1st Qu. 52.25 87.25 188.0 242.0 277.0
## Median 53.00 88.00 188.0 242.0 277.5
## Mean 52.80 87.90 188.3 242.1 277.5
## 3rd Qu. 53.00 88.00 189.0 242.0 278.0
## Max. 54.00 89.00 189.0 243.0 279.0## g1s1 g1s2 g1s3 g1s4 g1s5 g1s6
## Min. 53.00 88.00 152 213.0 239.0 303.0
## 1st Qu. 53.00 88.25 153 214.0 240.0 304.0
## Median 53.00 89.00 153 214.5 240.5 304.5
## Mean 53.40 88.90 153 214.3 240.4 304.6
## 3rd Qu. 53.75 89.00 153 215.0 241.0 305.0
## Max. 55.00 90.00 154 215.0 241.0 306.0A more detailed comparison would not be meaningful here because the time series is theoretical and the choice of models used is arbitrary. Nevertheless, these simulations highlight notable differences in phenology for H. armigera. A comparison could be made based on meteorological data from sites where insect trapping is carried out, in order to compare simulated data with observed data and choose the most appropriate model from among the many possible combinations.