---
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")
```
