PhenoFlex Model constitutes a framework for modeling the spring phenology of deciduous trees, with a particular focus on fruit and nut trees. It uses the structure of the Dynamic Model for chill accumulation and the Growing Degree Hours model for heat accumulation. In order to relate spring phenology to chill and heat, and to account for varying theories about the relationship between these two agroclimatic phenomena,
PhenoFlex includes a flexible transition function that defines responsiveness to heat (during ecodormancy) as a function of accumulated chill (during endodormancy). Unlike most previous applications of these models,
PhenoFlex can flexibly fit the parameters of both models to observed spring phenology dates. The model can thus accommodate variation among the temperature responses of different species and cultivars.
PhenoFlex Model is implemented in the
chillR package, which also contains functions to fit the model parameters to observed spring phenology data. This vignette demonstrates the use of the PhenoFlex Model to predict spring phenology dates, as well as the procedure for fitting model parameters to data.
We demonstrate the
PhenoFlex functions using the example of cherry bloom data recorded at Campus Klein-Altendorf, the experimental station of the University of Bonn. This dataset, along with long-term records of daily minimum and maximum temperatures, is included in the
chillR package (
KA_weather, respectively). Since the PhenoFlex model requires hourly temperatures, we use the
stack_hourly_temps function from
chillR to interpolate the daily data.
library(chillR) library(ggplot2) data(KA_weather) data(KA_bloom) <- stack_hourly_temps(KA_weather, latitude=50.4)hourtemps
PhenoFlex, the chilling requirement is denoted \(y_c\) and the heat requirement \(z_c\). We first illustrate how the model can be used to predict spring phenology events, based on user-specified \(y_c\) and \(z_c\) values and using the usual parameters for the Dynamic Model and the Growing Degree Hours model (the default values in
PhenoFlex). To run the analysis for one year only, we select all data for the 2009 season. The 2009 season is the dormancy season that ends in 2009. The number of months to consider for the season is specified in the
genSeason function by the
mrange parameter. Since the default (August to June) is appropriate here, this is not specified below.
<- 40 yc <- 190 zc <- genSeason(hourtemps, years=c(2009)) iSeason <- PhenoFlex(temp=hourtemps$hourtemps$Temp[iSeason[]], res times=c(1: length(hourtemps$hourtemps$Temp[iSeason[]])), zc=zc, stopatzc=TRUE, yc=yc, basic_output=FALSE)
The PhenoFlex function generates a list that tracks the accumulation of x (precursor to the dormancy-breaking factor), y (the dormancy-breaking factor; or Chill Portion), z (heat accumulation) and xs (the ratio of the formation to the destruction rate of x) over time. It also returns a bloomindex element, which points to the row in the temps input table that corresponds to estimated budbreak (or whatever the phenological stage of interest is).
<- res$bloomindex DBreakDay <-hourtemps$hourtemps[iSeason[],] seasontemps"x"]<-res$x seasontemps[,"y"]<-res$y seasontemps[,"z"]<-res$z seasontemps[,<-add_date(seasontemps) seasontemps ggplot(data=seasontemps[1:DBreakDay,],aes(x=Date,y=y)) + geom_line(col="blue",lwd=1.5) + theme_bw(base_size=20) + geom_hline(yintercept=yc,lty=2) + labs(title="Chill (y) accumulation")
ggplot(data=seasontemps[1:DBreakDay,],aes(x=Date,y=z)) + geom_line(col="red",lwd=1.5) + theme_bw(base_size=20) + geom_hline(yintercept=zc,lty=2) + labs(title="Heat (z) accumulation")
PhenoFlex model can be fitted to phenological data, provided that sufficient observations are available for the cultivar of interest. For this kind of model, parameters are usually determined using an empirical solver. Solvers identify the best combination of model parameters by trying out different values and gradually adjusting these parameters, until the objective function - a measure of how far predictions are from the observed data - does not decrease further. In this framework, we fit the model using a generalized simulated annealing algorithm. In contrast to other algorithms, simulated annealing can deal with discrete data (we calculate the residual sum of squares between observed and predicted bloom days as objective function). Since simulated annealing is a stochastic solver, there is a risk of not finding the overall best parameter combination (global minimum of residual sum of squares). The generalisation of the basic solver reduces this risk. Still, the search should be repeated several times with different initial parameter values and random seeds. This iterative procedure requires lots of model runs, which can take substantial time and/or computing power. Here, we only demonstrate the functionality using a maximum of 10 iterations of the simulated annealing procedure. In order to achieve reliable parameters, we recommend at least 1000 iterations. We also recommend using many observations of a cultivar’s spring phenology (many more than in this example), and ideally to include winter seasons spanning a wide range of temperature conditions.
For fitting the
PhenoFlex model to phenological data, we have to generate a list of seasons as follows. We are only using 6 seasons in this example, but a real-life application should be based on at least 15 to 20 seasons or even more, comprising variable temperature profiles across years.
<- genSeasonList(hourtemps$hourtemps, years=c(2003:2008))SeasonList
Parameters are then fitted using the generalized simulated annealing algorithm, which is called by the
phenologyFitter function and requires several inputs:
PhenoFlexmodel has more than one parameter, we have to wrap it into another function (
PhenoFlex_GDHwrapper), which requires a temperature dataset
xand a vector containing all the parameters of the
PhenoFlexmodel. Within the wrapper function, this parameter vector is unpacked, with each value assigned to the respective parameter.
The order, in which the
PhenoFlex parameters have to be provided, is given in the description of the
PhenoFlex_GDHwrapper function as follows:
PhenoFlex- default value: 0.5
Note that the
PhenoFlex_GDHwrapper can only fit the version of the
PhenoFlex model that uses a heat model based on the GDH concept. For use of the Gaussian heat model that can also be considered by
PhenoFlex, use the
For this example, we initiate the procedure with the default values, and we set upper and lower bounds as follows:
<- c(40, 190, 0.5, 25, 3372.8, 9900.3, 6319.5, 5.939917e13, 4, 36, 4, 1.60) par <- c(41, 200, 1.0, 30, 4000.0, 10000.0, 7000.0, 6.e13, 10, 40, 10, 50.00) upper <- c(38, 180, 0.1, 0 , 3000.0, 9000.0, 6000.0, 5.e13, 0, 0, 0, 0.05)lower
Now, we can run the fitter:
<- phenologyFitter(par.guess=par, Fit_res modelfn = PhenoFlex_GDHwrapper, bloomJDays=KA_bloom$pheno[which(KA_bloom$Year %in% c(2003:2008))], SeasonList=SeasonList, lower=lower, upper=upper, control=list(smooth=FALSE, verbose=FALSE, maxit=10, nb.stop.improvement=5))
## Loading required namespace: parallel
Note that the list
control regulates the behavior of the
GenSA algorithm, the solving engine we are using here. Note further that
smooth=FALSE must be set for the
PhenoFlex model, since the objective function, the residual sum of squares between predicted and observed bloom days is calculated based on discrete data (days) and is thus not smooth.
verbose controls whether messages from the algorithm are shown.
maxit defines the maximum number of iterations of the algorithm and should have a value of at least 1000.
nb.stop.improvement is the number of search steps of the algorithm, after which the solver is stopped when there is no further improvement of the fit.
In the example above, the values are chosen in a way that allows the fitter to finish quickly. Thus, the result may not be particularly meaningful. More reasonable parameters are, for instance, the default values of
=list(smooth=FALSE, verbose=TRUE, maxit=1000, controlnb.stop.improvement=250)
phenologyFitter can also be used with other models. For instance, in order to fit the
StepChill model, we could add
to the argument list of
phenologyFitter and adapt the parameter vectors
The result of a fitting procedure can be summarized as follows
Phenology Fitter Final RSS: 114.4583 RMSE : 4.367653 R^2 : 3.28061 data versus predicted: data predicted delta 1 107 105.0417 1.958333 2 110 106.1250 3.875000 3 105 100.0417 4.958333 4 116 117.2500 -1.250000 5 105 105.1250 -0.125000 6 115 123.3333 -8.333333
data being the days where bloom or any other phenological event was observed,
predicted the day that was predicted by the model and
delta the difference between
predicted, i.e. the residual error. These results can also be visualized using the
Model prediction are usually uncertain, and the possible errors this may produce should be expressed. To estimate these errors, we use a bootstrapping technique, which involves the following steps:
PhenoFlexframework, make predictions for the validation dataset and record the results for each year
PhenoFlexconcept will eventually be provided in a peer-reviewed paper
In short, the residuals are bootstrapped and then the fit is repeated. This procedure is performed
<- bootstrap.phenologyFit(Fit_res, boot.R=10, Fit_res.boot control=list(smooth=FALSE, verbose=TRUE, maxit=10, nb.stop.improvement=5), lower=lower, upper=upper, seed=1726354)
## ========It: 2, obj value: 202.8385417 ## ================================================It: 1, obj value: 96.80208333 ## ========It: 1, obj value: 29.62152778 ## ================
Like for the result of
phenologyFitter there is also a
summary and a
plot function for
par Err q16 q84 1 3.951783e+01 7.892125e-01 3.863460e+01 4.017894e+01 2 1.807182e+02 4.716312e+00 1.807182e+02 1.898302e+02 3 5.818040e-01 1.693573e-01 5.818040e-01 8.276560e-01 4 2.286709e+01 1.596126e+00 2.286709e+01 2.286709e+01 5 3.289466e+03 1.631851e+01 3.289466e+03 3.289466e+03 6 9.822295e+03 3.614660e+01 9.822295e+03 9.822295e+03 7 6.675057e+03 1.773010e+02 6.534836e+03 6.748974e+03 8 5.939918e+13 6.158073e+07 5.939917e+13 5.939925e+13 9 4.797002e+00 1.855391e+00 4.146226e+00 5.048300e+00 10 3.175800e+01 3.997594e+00 3.021162e+01 3.391548e+01 11 6.896034e+00 1.101434e+00 6.896034e+00 6.896034e+00 12 1.047828e+01 9.509628e+00 1.047828e+01 2.880682e+01
par being the estimated parameter values,
Err the standard deviation of the bootstrap distribution and
q84 the 16. and 84. percentiles, respectively.
This plot now shows observed bloom dates, as well as bloom dates predicted with the
PhenoFlex model, including an estimate of prediction uncertainty (shown as error bars).