# Function for Implementing the Accumulated Developmental Progress Method Using Mean Daily Temperatures Estimates the starting date ($S$, in day-of-year) and the parameters of a developmental rate model in the accumulated developmental progress (ADP) method using mean daily air temperatures (Wagner et al., 1984; Shi et al., 2017a, b). ```r ADP( S.arr, expr, ini.val, Year1, Time, Year2, DOY, Temp, DOY.ul = 120, fig.opt = TRUE, control = list(), verbose = TRUE ) ``` ## Arguments - `S.arr`: the candidate starting dates for thermal accumulation (in day-of-year) - `expr`: a user-defined model that is used in the accumulated developmental progress (ADP) method - `ini.val`: a vector or a list that saves the initial values of the parameters in `expr` - `Year1`: the vector of the years in which a particular phenological event was recorded - `Time`: the vector of the occurrence times (in day-of-year) of a particular phenological event across many years - `Year2`: the vector of the years recording the climate data corresponding to the occurrence times - `DOY`: the vector of the dates (in day-of-year) for which climate data exist - `Temp`: the mean daily air temperature data (in ${}^{\circ}$C) corresponding to `DOY` - `DOY.ul`: the upper limit of `DOY` used to predict the occurrence time - `fig.opt`: an optional argument to draw the figures associated with the temperature-dependent developmental rate curve, the mean daily temperatures versus years, and a comparison between the predicted and observed occurrence times - `control`: the list of control parameters for using the `optim` function in the `stats` package - `verbose`: an optional argument allowing users to suppress the printing of computation progress ## Details It is better not to set too many candidate starting dates, as doing so will be time-consuming. If `expr` is selected as Arrhenius' equation, `S.arr` can be selected as the `S` obtained from the output of the `ADTS` function. Here, `expr` can be other nonlinear temperature-dependent developmental rate functions (see Shi et al. [2017b] for details). Further, `expr` can be any an arbitrary user-defined temperature-dependent developmental rate function, e.g., a function named `myfun`, but it needs to take the form of `myfun <- function(P, x){...}`, where `P` is the vector of the model parameter(s), and `x` is the vector of the predictor variable, i.e., the temperature variable. $\qquad$The function does not require that `Year1` is the same as `unique(Year2)`, and the intersection of the two vectors of years will be kept. The unused years that have phenological records but lack climate data will be showed in `unused.years` in the returned list. $\qquad$The numerical value of `DOY.ul` should be greater than or equal to the maximum `Time`. $\qquad$ Let $r$ represent the temperature-dependent developmental rate, i.e., the reciprocal of the developmental duration required for completing a particular phenological event, at a constant temperature. In the accumulated developmental progress (ADP) method, when the annual accumulated developmental progress (AADP) reaches 100%, the phenological event is predicted to occurr for each year. Let $\mathrm{AADP}_{i}$ denote the AADP of the $i$th year, which equals $$ \mathrm{AADP}_{i} = \sum_{j=S}^{E_{i}}r_{ij}\left(\mathrm{\mathbf{P}}; T_{ij}\right), $$ where $S$ represents the starting date (in day-of-year), $E_{i}$ represents the ending date (in day-of-year), i.e., the occurrence time of a pariticular phenological event in the $i$th year, $\mathrm{\mathbf{P}}$ is the vector of the model parameters in `expr`, and $T_{ij}$ represents the mean daily temperature of the $j$th day of the $i$th year (in ${}^{\circ}$C or K). In theory, $\mathrm{AADP}_{i} = 100\%$, i.e., the AADP values of different years are a constant 100%. However, in practice, there is a certain deviation of $\mathrm{AADP}_{i}$ from 100%. The following approach is used to determine the predicted occurrence time. When $\sum_{j=S}^{F}r_{ij} = 100\%$ (where $F \geq S$), it follows that $F$ is the predicted occurrence time; when $\sum_{j=S}^{F}r_{ij} < 100\%$ and $\sum_{j=S}^{F+1}r_{ij} > 100\%$, the trapezoid method (Ring and Harris, 1983) is used to determine the predicted occurrence time. Let $\hat{E}_{i}$ represent the predicted occurrence time of the $i$th year. Assume that there are $n$-year phenological records. When the starting date $S$ and the temperature-dependent developmental rate model are known, the model parameters can be estimated using the Nelder-Mead optiminization method (Nelder and Mead, 1965) to minimize the root-mean-square error (RMSE) between the observed and predicted occurrence times, i.e., $$ \mathrm{\mathbf{\hat{P}}} = \mathrm{arg}\,\mathop{\mathrm{min}}\limits_{\mathrm{\mathbf{P}}}\left\{\mathrm{RMSE}\right\} = \mathrm{arg}\,\mathop{\mathrm{min}}\limits_{\mathrm{\mathbf{P}}}\sqrt{\frac{\sum_{i=1}^{n}\left(E_{i}-\hat{E}_i\right){}^{2}}{n}}. $$ Because $S$ is not determined, a group of candidate $S$ values (in day-of-year) need to be provided. Assume that there are $m$ candidate $S$ values, i.e., c("$S_{1}, S_{2}, S_{3}, \n$", "$\\cdots, S_{m}$"). For each $S_{q}$ (where $q$ ranges between 1 and $m$), we can obtain a vector of the estimated model parameters, $\mathrm{\mathbf{\hat{P}}}_{q}$, by minimizing $\mathrm{RMSE}_{q}$ using the Nelder-Mead optiminization method. Then we finally selected $\mathrm{\mathbf{\hat{P}}}$ associated with c("$\\mathrm{min}\\left\\{\\mathrm{RMSE}_{1}, \n$", "$\\mathrm{RMSE}_{2}, \\mathrm{RMSE}_{3}, \\cdots, \\mathrm{RMSE}_{m}\\right\\}$") as the target parameter vector. ## Returns - **TDDR**: the temperature-dependent developmental rate matrix consisting of the year, day-of-year, mean daily temperature and developmental rate columns - **MAT**: a matrix consisting of the candidate starting dates and the estimates of candidate model parameters with the corresponding RMSEs - **Dev.accum**: the calculated annual accumulated developmental progresses in different years - **Year**: The overlapping years between `Year1` and `Year2` - **Time**: The observed occurrence times (day-of-year) in the overlapping years between `Year1` and `Year2` - **Time.pred**: the predicted occurrence times in different years - **S**: the determined starting date (day-of-year) - **par**: the estimates of model parameters - **RMSE**: the RMSE (in days) between the observed and predicted occurrence times - **unused.years**: the years that have phenological records but lack climate data ## References Nelder, J.A., Mead, R. (1965) A simplex method for function minimization. **Computer Journal** 7, 308$-$313. tools:::Rd_expr_doi("10.1093/comjnl/7.4.308") Ring, D.R., Harris, M.K. (1983) Predicting pecan nut casebearer (Lepidoptera: Pyralidae) activity at College Station, Texas. **Environmental Entomology** 12, 482$-$486. tools:::Rd_expr_doi("10.1093/ee/12.2.482") Shi, P., Chen, Z., Reddy, G.V.P., Hui, C., Huang, J., Xiao, M. (2017a) Timing of cherry tree blooming: Contrasting effects of rising winter low temperatures and early spring temperatures. **Agricultural and Forest Meteorology** 240$-$241, 78$-$89. tools:::Rd_expr_doi("10.1016/j.agrformet.2017.04.001") Shi, P., Fan, M., Reddy, G.V.P. (2017b) Comparison of thermal performance equations in describing temperature-dependent developmental rates of insects: (III) Phenological applications. **Annals of the Entomological Society of America** 110, 558$-$564. tools:::Rd_expr_doi("10.1093/aesa/sax063") Wagner, T.L., Wu, H.-I., Sharpe, P.J.H., Shcoolfield, R.M., Coulson, R.N. (1984) Modelling insect development rates: a literature review and application of a biophysical model. **Annals of the Entomological Society of America** 77, 208$-$225. tools:::Rd_expr_doi("10.1093/aesa/77.2.208") ## Author(s) Peijian Shi pjshi@njfu.edu.cn , Zhenghong Chen chenzh64@126.com , Jing Tan jmjwyb@163.com , Brady K. Quinn Brady.Quinn@dfo-mpo.gc.ca . ## Note The entire mean daily temperature data set for the spring of each year should be provided. In `TDDR`, the first column of `Year` saves the years, the second column of `DOY` saves the day-of-year values, the third column of `Temperature` saves the mean daily air temperatures calculated between the starting date to the occurrence times, and the fourth column of `Rate` saves the calculated developmental rates corresponding to the mean daily temperatures. ## See Also `predADP` ## Examples ```r data(apricotFFD) data(BJDAT) X1 <- apricotFFD X2 <- BJDAT Year1.val <- X1$Year Time.val <- X1$Time Year2.val <- X2$Year DOY.val <- X2$DOY Temp.val <- X2$MDT DOY.ul.val <- 120 S.arr0 <- 47 #### Defines a re-parameterized Arrhenius' equation ########################### Arrhenius.eqn <- function(P, x){ B <- P[1] Ea <- P[2] R <- 1.987 * 10^(-3) x <- x + 273.15 10^12*exp(B-Ea/(R*x)) } ############################################################################## #### Provides the initial values of the parameter of Arrhenius' equation ##### ini.val0 <- list( B = 20, Ea = 14 ) ############################################################################## res5 <- ADP( S.arr = S.arr0, expr = Arrhenius.eqn, ini.val = ini.val0, Year1 = Year1.val, Time = Time.val, Year2 = Year2.val, DOY = DOY.val, Temp = Temp.val, DOY.ul = DOY.ul.val, fig.opt = TRUE, control = list(trace = FALSE, reltol = 1e-12, maxit = 5000), verbose = TRUE ) res5 TDDR <- res5$TDDR T <- TDDR$Temperature r <- TDDR$Rate Y <- res5$Year DP <- res5$Dev.accum dev.new() par1 <- par(family="serif") par2 <- par(mar=c(5, 5, 2, 2)) par3 <- par(mgp=c(3, 1, 0)) Ind <- sort(T, index.return=TRUE)$ix T1 <- T[Ind] r1 <- r[Ind] plot( T1, r1, cex.lab = 1.5, cex.axis = 1.5, pch = 1, cex = 1.5, col = 2, type = "l", xlab = expression(paste("Mean daily temperature (", degree, "C)", sep = "")), ylab = expression(paste("Calculated developmental rate (", {day}^{"-1"}, ")", sep = "")) ) par(par1) par(par2) par(par3) dev.new() par1 <- par(family="serif") par2 <- par(mar=c(5, 5, 2, 2)) par3 <- par(mgp=c(3, 1, 0)) plot( Y, DP * 100, xlab = "Year", ylab = "Accumulated developmental progress (%)", ylim = c(50, 150), cex.lab=1.5, cex.axis = 1.5, cex = 1.5 ) abline( h = 1 * 100, lwd = 1, col = 4, lty = 2 ) par(par1) par(par2) par(par3) # graphics.off() ```

ADP function