Model evaluation

In devRate: Quantify the Relationship Between Development Rate and Temperature in Ectotherms

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^[Shi, P. J., Reddy, G. V., Chen, L., & Ge, F. (2015). Comparison of thermal performance equations in describing temperature-dependent developmental rates of insects:(I) empirical models. Annals of the Entomological Society of America, 109(2), 211-215.].

require("devRate")

Nonlinear models

Table 1. Models used by Shi et al. 2016

| Model name | Model name in devRate | |----------------|--------------------------| | Briere-1 |briere1_99 | | Briere-2 |briere2_99 | | Lactin |lactin2_95 | | Performance-2 |perf2_11 | | Beta |beta_16 | | Ratkowsky |ratkowsky_83 |

Datasets

Table 2. Datasets of temperature-dependent development rates used by Shi et al. 2016 and method used to retrieve the source dataset.

| Sample species | Source | Method | |------------------------------|-------------------------------|-------------------| |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))

Fitting the six models to the eight empirical datasets

[1] Briere-1 model

Starting values for the NLS fit were chosen on the basis of available information on the package database with the briere1_99 object.

The devRateModel function can be applied to each dataset, or a list containing all the datasets can be built to ease the process. In the first case:

wuDS_briere1_99 <- devRateModel(eq = briere1_99, 
  dfData = wuDS,
  startValues = list(aa = 0.01, Tmin = 10, Tmax = 40))

broufasDS_briere1_99 <- devRateModel(eq = briere1_99, 
  dfData = broufasDS,
  startValues = list(aa = 0.01, Tmin = 10, Tmax = 40))

trpisDS_briere1_99 <- devRateModel(eq = briere1_99, 
  dfData = trpisDS,
  startValues = list(aa = 0.01, Tmin = 10, Tmax = 40))

messengerDS_briere1_99 <- devRateModel(eq = briere1_99, 
  dfData = messengerDS,
  startValues = list(aa = 0.01, Tmin = 10, Tmax = 40))

liu1DS_briere1_99 <- devRateModel(eq = briere1_99, 
  dfData = liu1DS,
  startValues = list(aa = 0.01, Tmin = 10, Tmax = 40))

liu2DS_briere1_99 <- devRateModel(eq = briere1_99, 
  dfData = liu2DS,
  startValues = list(aa = 0.01, Tmin = 10, Tmax = 40))

shiraiDS_briere1_99 <- devRateModel(eq = briere1_99, 
  dfData = shiraiDS,
  startValues = list(aa = 0.01, Tmin = 10, Tmax = 40))

deJongDS_briere1_99 <- devRateModel(eq = briere1_99, 
  dfData = deJongDS,
  startValues = list(aa = 0.01, Tmin = 10, Tmax = 40))

The results can then be stored in a single list.

briere1NLS <- list(wuDS_briere1_99, broufasDS_briere1_99, trpisDS_briere1_99, 
                   messengerDS_briere1_99, liu1DS_briere1_99, liu2DS_briere1_99, 
                   shiraiDS_briere1_99, deJongDS_briere1_99)

Or alternatively:

listDS <- list(wuDS, broufasDS, trpisDS, messengerDS, liu1DS, liu2DS, shiraiDS, deJongDS)
briere1NLS <- lapply(listDS, function(myDataSet){
  devRateModel(eq = briere1_99, 
    dfData = myDataSet,
    startValues = list(aa = 0.01, Tmin = 10, Tmax = 40)
  )
})

[2] Briere-2 model

Likewise, the briere2_99 model can be fitted to the dataset each one at a time, or using the lapply function.

broufasDS_briere2_99 <- devRateModel(eq = briere2_99, 
  dfData = broufasDS,
  startValues = list(aa = 0.01, Tmin = 10, Tmax = 40, bb = 2))
### and so on ...

briere2NLS <- lapply(listDS, function(myDataSet){
  devRateModel(eq = briere2_99, 
    dfData = myDataSet,
    startValues = list(aa = 0.01, Tmin = 10, Tmax = 40, bb = 2)
  )
})

[3] Lactin model

broufasDS_lactin2_95 <- devRateModel(eq = lactin2_95, 
  dfData = broufasDS,
  startValues = list(aa = 0.03, Tmax = 30, deltaT = 5.0, bb = -1.5))
### and so on ...

lactin2NLS <- lapply(listDS, function(myDataSet){
  devRateModel(eq = lactin2_95,
    dfData = myDataSet,
    startValues = list(aa = 0.03, Tmax = 30, deltaT = 5.0, bb = -1.5)
  )
})

[4] Performance-2 Model

perf2NLS <- lapply(listDS, function(myDataSet){
  devRateModel(eq = perf2_11, 
    dfData = myDataSet,
    startValues = list(cc = 0.02, T1 = 10, k = 0.5, T2 = 35)
  )
})

[5] Beta model

betaNLS <- lapply(listDS, function(myDataSet){
  devRateModel(eq = beta_16, 
    dfData = myDataSet,
    startValues = list(rm = 0.2, T1 = 5, T2 = 40, Tm = 30)
  )
})

[6] Ratkowsky model

ratkowskyNLS <- lapply(listDS, function(myDataSet){
  devRateModel(eq = ratkowsky_83, 
    dfData = myDataSet,
    startValues = list(cc = 0.02, T1 = 5, T2 = 40, k = 0.2)
  )
})

Organizing models' results

Model results can be grouped into a single object to ease the data processing.

listNLS <- list(briere1NLS, briere2NLS, lactin2NLS, perf2NLS, betaNLS, ratkowskyNLS)

Model evaluation

Shi et al. 2016 used the RSS, R2, R2 adjusted, RMSE, and AIC corrected in their study.

[1] RSS

The RSS can be computed from the nls objects returned by the devRateModel function, and equation (7) from Shi et al. 2016.

\begin{equation} RSS = \sum_{i=1}^n (obs_i - pred_i)^2 \end{equation}

RSS <- t(sapply(1:6, function(myModel){
  sapply(1:8, function(myDS){
    return(sum((listDS[[myDS]][,2] - predict(listNLS[[myModel]][[myDS]]))^2))
  })
}))
rownames(RSS) <- c("Briere-1", "Briere-2", "Lactin", "Perf-2", "Beta", "Ratkowsky")
colnames(RSS) <- c("wu", "broufas", "trpis", "messenger", "liu1", "liu2", "shirai", "deJong")

Table 3. RSS results

knitr::kable(RSS) 

[2] R squared

The R2 can be computed using equation (8) from Shi et al. 2016.

\begin{equation} R^2 = 1 - \frac{\sum_{i=1}^n (obs_i - pred_i)^2}{\sum_{i=1}^n (obs_i - meanObs_i)^2} \end{equation}

R2 <- t(sapply(1:6, function(myModel){
  sapply(1:8, function(myDS){
    return(1 - sum((listDS[[myDS]][,2] - predict(listNLS[[myModel]][[myDS]]))^2) / 
      sum((listDS[[myDS]][,2] - mean(listDS[[myDS]][,2]))^2))
  })
}))
rownames(R2) <- c("Briere-1", "Briere-2", "Lactin", "Perf-2", "Beta", "Ratkowsky")
colnames(R2) <- c("wu", "broufas", "trpis", "messenger", "liu1", "liu2", "shirai", "deJong")

Table 4. R2 results

knitr::kable(R2) 

[3] R squared adjusted

The R2 adjusted can be computed using equation (9) from Shi et al. 2016.

\begin{equation} R^2_{adj} = 1 - \frac{n-1}{n-p}(1-R^2) \end{equation}

R2adj <- t(sapply(1:6, function(myModel){
  p <- 1 + length(coef(listNLS[[myModel]][[1]]))
  sapply(1:8, function(myDS){
    n <- length(listDS[[myDS]][,2])
    Rsq <- 1 - sum((listDS[[myDS]][,2] - predict(listNLS[[myModel]][[myDS]]))^2) / 
      sum((listDS[[myDS]][,2] - mean(listDS[[myDS]][,2]))^2)
    return(1 - (n - 1)/(n - p) * (1 - Rsq))
  })
}))
rownames(R2adj) <- c("Briere-1", "Briere-2", "Lactin", "Perf-2", "Beta", "Ratkowsky")
colnames(R2adj) <- c("wu", "broufas", "trpis", "messenger", "liu1", "liu2", "shirai", "deJong")

Table 5. R2 adjusted results

knitr::kable(R2adj) 

[4] RMSE

The RMSE can be computed using equation (10) from Shi et al. 2016.

\begin{equation} RMSE = \sqrt{RSS/(n-p+1)} \end{equation}

RMSE <- t(sapply(1:6, function(myModel){
  p <- 1 + length(coef(listNLS[[myModel]][[1]]))
  sapply(1:8, function(myDS){
    n <- length(listDS[[myDS]][,2])
    RSS <- sum((listDS[[myDS]][,2] - predict(listNLS[[myModel]][[myDS]]))^2)
    return(sqrt(RSS / (n - p + 1)))
  })
}))
rownames(RMSE) <- c("Briere-1", "Briere-2", "Lactin", "Perf-2", "Beta", "Ratkowsky")
colnames(RMSE) <- c("wu", "broufas", "trpis", "messenger", "liu1", "liu2", "shirai", "deJong")

Table 6. RMSE results

knitr::kable(RMSE) 

[5] AIC corrected

The AICc can be computed using equations (11) and (12) from Shi et al. 2016.

\begin{equation} AIC_c = -2L+2pn/(n-p-1) \end{equation}

\begin{equation} L = -\frac{n}{2}ln\frac{RSS}{n} \end{equation}

AICc <- t(sapply(1:6, function(myModel){
  p <- 1 + length(coef(listNLS[[myModel]][[1]]))
  sapply(1:8, function(myDS){
    n <- length(listDS[[myDS]][,2])
    RSS <- sum((listDS[[myDS]][,2] - predict(listNLS[[myModel]][[myDS]]))^2)
    L <- -n/2 * log(RSS / n)
    return(-2 * L + 2 * p * n / (n - p - 1))
  })
}))
# number of parameters + 1
p <- sapply(1:6, function(myModel){
  p <- 1 + length(coef(listNLS[[myModel]][[1]]))
  return(p)
})
# sample size
n <- sapply(1:8, function(myDS){
  n <- length(listDS[[myDS]][,2])
  return(n)
})
rownames(AICc) <- paste0(c("Briere-1", "Briere-2", "Lactin", 
      "Perf-2", "Beta", "Ratkowsky"), " (", p, ")")
colnames(AICc) <- paste0(c("wu", "broufas", "trpis", "messenger", 
      "liu1", "liu2", "shirai", "deJong"), " (", n, ")")

Table 7. AICc results with the number of parameters + 1 in parenthesis next to the model names and the sample size in parenthesis next to the sample names.

knitr::kable(AICc) 

[6] AIC

The AIC can also be computed with the AIC function in the stats package.

AIC <- t(sapply(1:6, function(myModel){
  sapply(1:8, function(myDS){
    return(AIC(listNLS[[myModel]][[myDS]]))
  })
}))
rownames(AIC) <- c("Briere-1", "Briere-2", "Lactin", "Perf-2", "Beta", "Ratkowsky")
colnames(AIC) <- c("wu", "broufas", "trpis", "messenger", "liu1", "liu2", "shirai", "deJong")

Table 8. AIC results

knitr::kable(AIC) 

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).

Conclusion

The first four datasets gave the same results as in Shi et al. 2016, and the last four different results, probably because the datasets retrieved from the source may differ (e.g., altae or apterae aphid data in Liu et al. 1999; 2000). This vignette exemplifies the use of the devRate package to compare models. For interpretations, see the extensive literature on model evaluation with RSS, R2, R2 adjusted, RMSE, AIC, AIC corrected available online, and the study by Shi et al. 2016.

devRate documentation built on Aug. 24, 2023, 9:07 a.m.