--- title: "Using devRate package to fit development rate models to an empirical dataset" author: "Francois Rebaudo, using data from Crespo-Perez et al. 2011 " date: "`r Sys.Date()`" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{main} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- This vignette exemplifies the use of the devRate package to fit a development rate model to an empirical dataset (laboratory experiments at air constant temperatures) and to build a simple phenological model. The preliminary step is to perform an experiment where the arthropod study model is reared at different constant air temperatures, and to monitor the time at which individuals change of life stage (e.g., from eggs to larva). The monitoring of life stages is commonly expressed in development rate, corresponding to the inverse of the development time needed to reach the next life stage. Development time is usually expressed in days. For example, assuming that an individual needs 10 days to develop from eggs to larva at 15 degrees Celsius, its development rate would be 1/10 = 0.1. In this vignette we illustrate the different functions using a dataset from the literature. The dataset for _T. solanivora_ was retrieved from Crespo-Perez et al. 2011^[Crespo-Pérez, V., Rebaudo, F., Silvain, J.-F. and Dangles, O. (2011) Modeling invasive species spread in complex landscapes: the case of potato moth in Ecuador. Landscape Ecology, 26, 1447–1461.], using Web Plot Digitizer^[Rohatgi, A. (2015) WebPlotDigitalizer: HTML5 Based Online Tool to Extract Numerical Data from Plot Images.]. The dataset is also included in the package in the exTropicalMoth object. This object is composed of the laboratory data and results of the modeling (see below). In this vignette we present how to organize your own dataset from scratch. ## Organizing the dataset In this example various individuals of _T. solanivora_ were reared at different temperatures and the development rates were monitored for three life stages (eggs, larvae, pupae): * the eggs were reared at 10, 13, 15, 15.5, 16, 17, 20, 25, 30, and 35 degrees * the larva were reared at 10, 13, 15, 15.5, 17, 20, 25, 30, and 35 degrees * the pupa were reared at 10, 13, 15, 15.5, 16, 17, 20, 25, 30, and 35 degrees. The dataset resulting from the rearing experiments looks like the following table (for one life stage) and represents the only input needed by the package devRate. | Temperature | Development Rate | |------------:|-----------------:| |10.0 |0.031 | |10.0 |0.039 | |13.0 |0.072 | |... |... | If the dataset is stored into a file, it can be read directly from that file (e.g., using _read.table()_). In this example, we copy-pasted the experimental results to create vectors and data frames that are the required input types for the package. ```{r} ### Experimental temperatures and development rate for the eggs expeTempEggs <- c(10.0, 10.0, 13.0, 15.0, 15.0, 15.5, 16.0, 16.0, 17.0, 20.0, 20.0, 25.0, 25.0, 30.0, 30.0, 35.0) expeDevEggs <- c(0.031, 0.039, 0.072, 0.047, 0.059, 0.066, 0.083, 0.1, 0.1, 0.1, 0.143, 0.171, 0.2, 0.2, 0.18, 0.001) dfDevEggs <- data.frame(expeTempEggs, expeDevEggs) ### Experimental temperatures and development rate for the larva expeTempLarva <- c(10.0, 10.0, 10.0, 13.0, 15.0, 15.5, 15.5, 15.5, 17.0, 20.0, 25.0, 25.0, 30.0, 35.0) expeDevLarva <- c(0.01, 0.014, 0.019, 0.034, 0.024, 0.029, 0.034, 0.039, 0.067, 0.05, 0.076, 0.056, 0.0003, 0.0002) dfDevLarva <- data.frame(expeTempLarva, expeDevLarva) ### Experimental temperatures and development rate for the pupa expeTempPupa <- c(10.0, 10.0, 10.0, 13.0, 15.0, 15.0, 15.5, 15.5, 16.0, 16.0, 17.0, 20.0, 20.0, 25.0, 25.0, 30.0, 35.0) expeDevPupa <- c(0.001, 0.008, 0.012, 0.044, 0.017, 0.044, 0.039, 0.037, 0.034, 0.051, 0.051, 0.08, 0.092, 0.102, 0.073, 0.005, 0.0002) dfDevPupa <- data.frame(expeTempPupa, expeDevPupa) ### Same dataset included in the package in the form of matrices library("devRate") data(exTropicalMoth) str(exTropicalMoth[[1]]) ``` ## Finding models in the literature Before attempting to fit any model to the empirical data, the devRate function "devRateFind" search the database for previous articles fitting models to the organism, either by Order, Family, or species. The most used models are presented at the top of the data.frame. ```{r} devRateFind(orderSP = "Lepidoptera") ``` Here we can see that at the time of this vignette, campbell_74 model was used 108 times to characterize the relationship between development rate and temperature for the Lepidoptera Order, and the taylor_81 model 22 times. ```{r} devRateFind(familySP = "Gelechiidae") ``` If we focus on the Family (Gelechiidae), campbell_74 has been used 12 times, lactin2_95 7 times, logan6_76 4 times, lactin1_95 4 times, briere1_99 4 times, damos_08 4 times, analytis_77 3 times, taylor_81 2 times, and so on (this may change as the database is growing). ```{r} devRateFind(species = "Tecia solanivora") ``` Unfortunately, the species _Tecia solanivora_ is not in the package database at this time, so that we have to find a closely related species. A rapid search in the literature indicates that _T. solanivora_ is often associated with two tuber moths: _Symmetrischema tangolias_ and _Phthorimaea operculella_, both being of the Gelechiidae family. The devRateFind function can be used to browse the database for these two Gelechiidae species. ```{r} devRateFind(species = "Symmetrischema tangolias") devRateFind(species = "Phthorimaea operculella") ``` The taylor_81 model was used for _S. tangolias_ and among others, the lactin1_95 model was used for _P. operculella_. The "devRateInfo" function provides additional information on these models and on parameter estimations. ```{r} devRateInfo(eq = taylor_81) devRateInfo(eq = lactin1_95) ``` Here we can see for example that _S. tangolias_ study by Sporleder et al. 2016 has used the taylor_81 model, and that _P. operculella_ study by Golizadeh et al. 2012 has used the lactin1_95 model. To return only the information on the closely related species, a specific query can be performed on the model. ```{r} taylor_81$startVal[taylor_81$startVal["genSp"] == "Symmetrischema tangolias",] lactin1_95$startVal[lactin1_95$startVal["genSp"] == "Phthorimaea operculella",] ``` Information from the database can be plotted using the "devRatePlotInfo" function if we want to have a first look on the taylor_81 or lactin1_95 models. ```{r, fig.width = 6, fig.height = 6} devRatePlotInfo (eq = taylor_81, sortBy = "ordersp", ylim = c(0, 0.20), xlim = c(0, 50)) devRatePlotInfo (eq = lactin1_95, sortBy = "ordersp", ylim = c(0, 1.00), xlim = c(0, 50)) ``` If there is no a priori model selection (e.g., guided by a specific research question), the taylor_81 and/or lactin1_95 models can be used as candidate models. ## Fitting models to empirical datasets Now that we have candidate models and starting parameter estimates from closely related species, we can start the fitting process with our empirical data (NLS fit). The empirical data can be fitted to any model in the database with the "devRateModel" function, which returns an object of class "nls". ```{r} ### using the vectors from section "Organizing the dataset" ############################################################ ### for the taylor_81 model mEggs01 <- devRateModel(eq = taylor_81, temp = expeTempEggs, devRate = expeDevEggs, startValues = list(Rm = 0.05, Tm = 30, To = 5)) mLarva01 <- devRateModel(eq = taylor_81, temp = expeTempLarva, devRate = expeDevLarva, startValues = list(Rm = 0.05, Tm = 25, To = 10)) mPupa01 <- devRateModel(eq = taylor_81, temp = expeTempPupa, devRate = expeDevPupa, startValues = list(Rm = 0.1, Tm = 30, To = 10)) ### for the lactin1_95 model mEggs01b <- devRateModel(eq = lactin1_95, temp = expeTempEggs, devRate = expeDevEggs, startValues = list(aa = 0.177, Tmax = 36.586, deltaT = 5.631)) # mLarva01b <- devRateModel(eq = lactin1_95, # temp = expeTempLarva, # devRate = expeDevLarva, # startValues = list(aa = 0.169, Tmax = 37.914, deltaT = 5.912)) ### The algorithm has not found a solution after 50 iterations ### One possibility is to increase the maximum number of iterations ### using the "control" argument (see ?nls() for more details). mLarva01b <- devRateModel(eq = lactin1_95, temp = expeTempLarva, devRate = expeDevLarva, startValues = list(aa = 0.169, Tmax = 37.914, deltaT = 5.912), control = list(maxiter = 500)) mPupa01b <- devRateModel(eq = lactin1_95, temp = expeTempPupa, devRate = expeDevPupa, startValues = list(aa = 0.193, Tmax = 36.291, deltaT = 5.18), control = list(maxiter = 500)) ### using the data frames from section "Organizing the dataset" ############################################################ mEggs02 <- devRateModel(eq = taylor_81, dfData = dfDevEggs, startValues = list(Rm = 0.05, Tm = 30, To = 5)) mLarva02 <- devRateModel(eq = taylor_81, dfData = dfDevLarva, startValues = list(Rm = 0.05, Tm = 25, To = 10)) mPupa02 <- devRateModel(eq = taylor_81, dfData = dfDevPupa, startValues = list(Rm = 0.1, Tm = 30, To = 10)) ### using the dataset included in the package (only for taylor_81 model) ############################################################ mEggs <- devRateModel(eq = taylor_81, temp = exTropicalMoth$raw$eggs[,1], devRate = exTropicalMoth$raw$eggs[,2], startValues = list(Rm = 0.05, Tm = 30, To = 5)) mLarva <- devRateModel(eq = taylor_81, temp = exTropicalMoth$raw$larva[,1], devRate = exTropicalMoth$raw$larva[,2], startValues = list(Rm = 0.05, Tm = 25, To = 10)) mPupa <- devRateModel(eq = taylor_81, temp = exTropicalMoth$raw$pupa[,1], devRate = exTropicalMoth$raw$pupa[,2], startValues = list(Rm = 0.1, Tm = 30, To = 10)) ``` The summary of the "devRateModel" can be obtained with the "devRatePrint" function. ```{r} resultNLS <- devRatePrint(myNLS = mLarva) resultNLSb <- devRatePrint(myNLS = mLarva01b) ``` Empirical data can be plotted against the model using the "devRatePlot" function. ```{r, fig.width = 7, fig.height = 5} par(mfrow = c(1, 2), mar = c(4, 4, 0, 0)) devRatePlot(eq = taylor_81, nlsDR = mEggs, pch = 16, ylim = c(0, 0.2)) devRatePlot(eq = lactin1_95, nlsDR = mEggs01b, pch = 16, ylim = c(0, 0.2)) devRatePlot(eq = taylor_81, nlsDR = mLarva, pch = 16, ylim = c(0, 0.1)) devRatePlot(eq = lactin1_95, nlsDR = mLarva01b, pch = 16, ylim = c(0, 0.1)) devRatePlot(eq = taylor_81, nlsDR = mPupa, pch = 16, ylim = c(0, 0.15)) devRatePlot(eq = lactin1_95, nlsDR = mPupa01b, pch = 16, ylim = c(0, 0.15)) ``` ## Models comparison Models can be compared using the AIC or any function compatible with an "nls" object, such as BIC or logLik (see the R manual for the use and interpretation of these functions, outside of the scope of this vignette). ```{r} ### Models for the larva life stage c(AIC(mLarva), AIC(mLarva01b)) c(BIC(mLarva), BIC(mLarva01b)) c(logLik(mLarva), logLik(mLarva01b)) ``` ## Forecasting phenologies from empirical temperature From the "nls" object obtained above, we can build a simple phenology model using temperature time series (e.g., to forecast the organism outbreaks). In this example the temperature dataset is built from a Normal distribution of mean 15 and a standard deviation of 1, with a time step of one day. The development models used are those previously fitted with the Taylor model for the three stages. We assumed that the average time for female adults to lay eggs was of 1 day. We simulated 50 individuals, with a stochasticity in development rate centered on the development rate, with a standard deviation of 0.015 (Normal distribution). ```{r} forecastTsolanivora <- devRateIBM( tempTS = rnorm(n = 100, mean = 15, sd = 1), timeStepTS = 1, models = list(mEggs, mLarva, mPupa), numInd = 50, stocha = 0.015, timeLayEggs = 1) print(forecastTsolanivora) ``` Results can be plotted using the "devRateIBMPlot" function. From this simple model we can observe that if eggs of the first generation are laid at time t = 0, adults should be seen at t = [65:85], if the temperature time series correspond to the temperatures experienced by the organism and in the absence of other limiting factors. ```{r, fig.width = 6, fig.height = 6} devRateIBMPlot(ibm = forecastTsolanivora, typeG = "density") ```