This vignette exemplifies the use of the devRate package i) to fit a
development rate models to empirical datasets with eight species using
six nonlinear models, and ii) to compare models based on goodness of fit
and the trade-off between the model’s goodness of fit and its structural
complexity.
It reproduces some of the results of the study by Shi et al. 2016.
Datasets
Table 2. Datasets of temperature-dependent
development rates used by Shi et al. 2016 and method used to retrieve
the source dataset.
| Helicoverpa armigera |
Wu et al. (2009) |
IPEC manual |
| Kampimodromus aberrans |
Broufas et al. (2007) |
Table 3 |
| Toxorhyynchites brevipalpis |
Trpis (1972) |
Table 1 |
| Bactrocera dorsalis |
Messenger and Flitters (1958) |
Table 1 |
| Aedes aegypti |
Gilpin and McClelland (1979) |
article not found |
| Bemisia tabaci (B-biotype) |
Xiang et al. (2007) |
article not found |
| Lipaphis erysimi |
Liu and Meng (2000) |
Table 1 |
| Myzus persicae |
Liu and Meng (1999) |
Table 1 |
| Epilachna varivestis |
Shirai and Yara (2001) |
Table 4 |
| Drosophila buzzatii |
de Jong (2010) |
Web Plot Digitizer |
Gilpin and McClelland (1979), and Xiang et al. (2007) articles were
not found and were excluded from this analysis.
The Wu et al. (2009) dataset was retrieved from the IPEC package
manual (article not found).
wuDS <- data.frame(
temp = seq(from = 15, to = 37, by = 1),
devRate = 1/c(41.24, 37.16, 32.47, 26.22, 22.71, 19.01, 16.79, 15.63, 14.27, 12.48, 11.3,
10.56, 9.69, 9.14, 8.24, 8.02, 7.43, 7.27, 7.35, 7.49, 7.63, 7.9, 10.03))
The Broufas et al. (2007) dataset was retrieved from Table 3 where
the mean development period is provided in hours and transformed in the
inverse of days.
broufasDS <- data.frame(
temp = c(15.0, 17.0, 20.0, 22.0, 25.0, 27.0, 30.0, 33.0, 35.0),
devRate = 1/(c(595.4, 463.0, 260.5, 222.5, 161.8, 153.1, 136.0, 147.8, 182.1)/24))
The Trpis et al. (2007) dataset was retrieved from Table 1 where the
mean egg development time is provided in hours and transformed in the
inverse of days.
trpisDS <- data.frame(
temp = c(14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0, 21.0, 22.0, 23.0, 24.0, 25.0, 26.0,
27.0, 28.0, 29.0, 30.0, 31.0, 32.0),
devRate = 1/(c(209.0, 174.0, 165.0, 125.0, 102.0, 90.0, 76.0, 62.0, 55.0, 50.0, 48.0,
44.0, 41.0, 39.0, 38.0, 37.5, 37.0, 38.0, 39.0)/24))
The Messenger and Flitters (1958) dataset was retrieved from Table 1
where the temperature is in Fahrenheit and transformed into Celsius and
the mean egg development time is provided in hours and transformed in
the inverse of days.
messengerDS <- data.frame(
temp = (c(55.0, 56.0, 57.0, 58.0, 60.0, 62.5, 65.0, 67.5, 70.0, 75.0, 80.0, 85.0, 87.5,
90.0, 92.5, 95.0, 96.0, 97.0, 97.5) - 32)/1.8,
devRate = 1/(c(263.0, 232.0, 170.5, 148.0, 121.3, 95.5, 74.0, 62.5, 51.5, 38.0, 30.5,
27.0, 25.0, 24.0, 23.5, 25.0, 26.5, 29.3, 34.3)/24))
The Liu and Meng (2000) dataset was retrieved from Table 1 (constant
temperatures), and for the altae form (this may be different from the
dataset used by Shi et al. 2016).
liu1DS <- data.frame(
temp = c(8.3, 11.3, 14.3, 17.3, 20.1, 22.3, 24.7, 26.5, 28.0, 30.0, 33.0, 35.1),
devRate = 1/c(47.5, 26.1, 15.8, 11.2, 8.3, 6.7, 6.2, 5.7, 5.4, 5.1, 5.6, 7.1))
The Liu and Meng (1999) dataset was retrieved from Table 1 for the
altae form (this may be different from the dataset used by Shi et
al. 2016).
liu2DS <- data.frame(
temp = c(6.2, 8.3, 11.3, 14.3, 17.3, 20.1, 22.3, 24.7, 26.5, 28.0, 30.0),
devRate = 1/c(50.4, 33.9, 21.5, 14.3, 11.2, 8.0, 6.8, 6.2, 5.8, 6.1, 6.5))
The Shirai and Yara (2001) dataset was retrieved from Table 4 (this
may be different from the dataset used by Shi et al. 2016).
shiraiDS <- data.frame(
temp = c(12.5, 15.0, 17.5, 20.0, 22.5, 25.0, 27.5, 30.0, 32.5),
devRate = 1/c(87.0, 55.5, 41.9, 32.4, 26.9, 21.7, 20.1, 20.6, 20.8))
The de Jong (2010) dataset was retrieved using Web Plot Digitizer on
Figure 1B (this should be different from the dataset used by Shi et
al. 2016).
deJongDS <- data.frame(
temp = c(15.1, 16.1, 17.1, 18.0, 21.7, 27.9, 30.1, 31.8),
devRate = c(0.0253, 0.0291, 0.0370, 0.0381, 0.0710, 0.1079, 0.1098, 0.0972))
Visualization
The fitted models can be plotted using the devRatePlot function from
the devRate package, or alternatively they can be plotted using the
predict function from the stats package. In the following code,
listDS[[4]] refers to the fourth dataset from Messenger and Flitters
(1958) and listNLS[[i]][[4]] refers to the corresponding model fit.
plot(x = listDS[[4]][,1],
y = listDS[[4]][,2],
ylim = c(0, 1), ylab = "Development rate",
xlim = c(5, 40), xlab = "Temperature", type = "n")
for(i in 1:6){
points(x = seq(from = 0, to = 50, by = 0.1),
y = predict(listNLS[[i]][[4]], newdata = list(T = seq(from = 0, to = 50, by = 0.1))),
type = 'l', lwd = 2, col = i)
}
points(x = listDS[[4]][,1], y = listDS[[4]][,2], pch = 16, cex = 1.5)
legend("topleft", col = 1:6, lwd = 2, legend = c("Briere-1", "Briere-2", "Lactin", "Perf-2", "Beta", "Ratkowsky"))
Figure 1. The six models fitted to the Messenger and
Flitters (1958) dataset.
plot(x = listDS[[3]][,1],
y = listDS[[3]][,2],
ylim = c(0, 0.7), ylab = "Development rate",
xlim = c(5, 40), xlab = "Temperature", type = "n")
for(i in 1:6){
points(x = seq(from = 0, to = 50, by = 0.1),
y = predict(listNLS[[i]][[3]], newdata = list(T = seq(from = 0, to = 50, by = 0.1))),
type = 'l', lwd = 2, col = i)
}
points(x = listDS[[3]][,1], y = listDS[[3]][,2], pch = 16, cex = 1.5)
legend("topleft", col = 1:6, lwd = 2, legend = c("Briere-1", "Briere-2", "Lactin", "Perf-2", "Beta", "Ratkowsky"))
Figure 2. The six models fitted to the Trpis (1972)
dataset.
getPlot <- function(datasetNumber, ...){
plot(x = listDS[[datasetNumber]][,1],
y = listDS[[datasetNumber]][,2],
ylab = "Development rate",
xlab = "Temperature", type = "n", ...)
for(i in 1:6){
points(x = seq(from = 0, to = 50, by = 0.1),
y = predict(listNLS[[i]][[datasetNumber]], newdata = list(T = seq(from = 0, to = 50, by = 0.1))),
type = 'l', lwd = 2, col = i)
}
points(x = listDS[[datasetNumber]][,1], y = listDS[[datasetNumber]][,2], pch = 16, cex = 1.5)
legend("topleft", col = 1:6, lwd = 2,
legend = c("Briere-1", "Briere-2", "Lactin", "Perf-2", "Beta", "Ratkowsky"))
}
par(mfrow = c(3, 2), mar = c(4, 4, 1, 1))
getPlot(datasetNumber = 1, ylim = c(0, 0.20), xlim = c(5, 40))
text(x = 40, y = 0.20, labels = "A", cex = 2, pos = 1)
getPlot(datasetNumber = 2, ylim = c(0, 0.25), xlim = c(5, 40))
text(x = 40, y = 0.25, labels = "B", cex = 2, pos = 1)
getPlot(datasetNumber = 5, ylim = c(0, 0.20), xlim = c(0, 40))
text(x = 40, y = 0.20, labels = "C", cex = 2, pos = 1)
getPlot(datasetNumber = 6, ylim = c(0, 0.20), xlim = c(0, 35))
text(x = 35, y = 0.20, labels = "D", cex = 2, pos = 1)
getPlot(datasetNumber = 7, ylim = c(0, 0.08), xlim = c(5, 40))
text(x = 40, y = 0.08, labels = "E", cex = 2, pos = 1)
getPlot(datasetNumber = 8, ylim = c(0, 0.15), xlim = c(5, 35))
text(x = 35, y = 0.15, labels = "F", cex = 2, pos = 1)
Figure 3. The six models fitted to the other
datasets: A: Wu et al. (2009), B: Broufas et al. (2007), C: Liu and Meng
(2000), D: Liu and Meng (1999), E: Shirai and Yara (2001), and F: de
Jong (2010).