Contents

• Introduction
• Mathematical Concepts for Disease Forecasting
• Introduction to Plant Disease Forecasting
• Case Study #1: Simple Generic Infection Model
• Case Study #2: Development and Field Evaluation of a Model
• Case Study #3: Modeling Population Dynamics Nematodes
• Final Conclusions
• Suggested Quiz and Study Questions
• References

Esker, P.D., A.H. Sparks, L. Campbell, Z. Guo, M. Rouse, S.D. Silwal, S. Tolos, B. Van Allen, and K.A. Garrett, 2008. Ecology and Epidemiology in R: Disease Forecasting. The Plant Health Instructor. DOI:10.1094/PHI-A-2008-0129-01.

Ecology and Epidemiology in R: Disease Forecasting

Mathematical Concepts for Disease Forecasting

Application of multiple regression for disease forecasting

In a simple linear regression, a response variable y is regressed on a single predictor (explanatory) variable x. In the context of disease forecasting, y might denote disease severity with x denoting a predictor like mean temperature. The model is y=β₀+β₁x+ε where ε is iid Normal (0,>σ²); see Sparks et al. (2008) for an introduction to linear regression with examples in R. Multiple linear regression generalizes this methodology to allow for multiple predictor variables, such as mean temperature and mean relative humidity. Notation for a model including m predictors is y=β₀+β₁x₁+β₂x₂+...+β_mx_m+ε.
One common multiple linear regression model is the polynomial model:
y=β₀+β₁x₁+β₂x₁²+β₃x₂+β₄x₂²+β₅x₁x₂+ε.
The estimate of the parameter can be calculated by redefining variables:
w₁=x₁, w₂=x₁², w₃=x₂, w₄=x₂², w₅=x₁x₂.
and performing a multiple linear regression model using
y=β₀+β₁w₁+β₂w₂+β₃w₃+β₄w₄+β₅w₅+ε.
This now is an ordinary multiple linear regression model using w as predictor variables.

The goal of multiple linear regression may be accurate prediction, or the regression coefficients may also be of direct interest themselves (Maindonald and Braun 2003).

Example of Multiple Linear Regression using R

The function used to fit multiple regression models is the same as for simple linear regression models:
lm(formula,data)

The following example, using simulated data, is an extension of the example provided for simple linear regression (Sparks et al. 2008). Here the response variable is disease severity, and the predictor variables are temperature, relative humidity and rainfall amount. A review of the use of the R programming environment may be helpful (Garrett et al. 2007).

disease<-c(1.9,3.1,3.3,4.8,5.3,6.1,6.4,7.6,9.8,12.4) #response variable
temperature<-c(62,61,65,65,80,80,83,70,90,85) #predictor variable
rel.humidity<-c(55,57,53,52,32,34,30,45,25,36) #predictor variable
rainfall<-c(1,1,0.5,0.9,0.1,0.2,0.8,0.4,0,0.2) #in inches
severity.mult <- as.data.frame(cbind(disease,temperature,rel.humidity,rainfall))

A scatterplot of all variables may be constructed using the following command:

pairs(severity.mult)

Output

This produces the following scatterplot matrix illustrating the relationship between each pair of variables:

Click to enlarge.

The patterns illustrated in the graph indicate an approximately straight-line relationship between each pair of variables. A scatterplot is a useful tool to examine not only how the potential predictor variables possibly relate to the response, but to also provide an exploratory method to determine if predictor variables are associated with one another (i.e., multicollinearity).

Next we fit a multiple linear regression model to the data with disease as a response variable, and temperature, rel.humidity, and rainfall as predictor variables.

##Multiple regression model
severity.mult.lm <- lm(disease~temperature+rel.humidity+rainfall,data=severity.mult)
summary(severity.mult.lm) #summary of the linear model

These commands produce the following output:

Output

Residuals:
    Min      1Q  Median      3Q     Max 
-1.5630 -0.7119 -0.2576  0.3774  2.9775 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept)  -74.9798    25.7954  -2.907   0.0271 *
temperature    0.7954     0.2342   3.396   0.0146 *
rel.humidity   0.5423     0.2037   2.663   0.0374 *
rainfall      -1.2061     2.1027  -0.574   0.5871  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 1.563 on 6 degrees of freedom
Multiple R-Squared: 0.8415,     Adjusted R-squared: 0.7623 
F-statistic: 10.62 on 3 and 6 DF,  p-value: 0.008172

In the output from this first analysis, the p-value for the full model is 0.008172, small enough to suggest that at least one of the predictor variables may be useful for prediction. The adjusted R² is 0.7622 which indicates that the predictor variables explain 76% of the variance in the response variable. The next step is to address the apparent correlation between the predictor variables.

Multicollinearity

When using multiple linear regression, multicollinearity, or correlation among predictor variables, is an important consideration. When multicollinearity is a problem, standard errors for the estimates will tend to be large and parameter estimates may be illogical. One can calculate a variance inflation factor (VIF) for each independent variable, where a VIF of five or greater indicates the existence of multicollinearity. A function vif() is available in R version 2.4.1. and the car library (requires downloading of the package from CRAN). In our example, the VIF values were 22.9 for temperature, 21.3 for relative humidity, and 2.4 for rainfall. The scatterplots had shown that temperature and relative humidity were probably strongly related, which may explain such large VIF values. For more conclusions on multicollinearity and how to handle it, see Heiburger and Holland (2003). In the example given above, there is some indication of multicollinearity. One way to try to eliminate multicollinearity problems is by doing variable selection to select a subset of predictor variables that are not strongly correlated.

One method to select appropriate predictor variables, ones that are useful for prediction and not strongly correlated with each other, is backward elimination (Neter et al. 1996). In this approach, all predictor variables are included in the starting models, and then particular predictor variables are removed to iteratively search for the best fit. The Akaike information criterion (AIC) can be used to select among models for an optimal combination of parsimony (limiting the model to the smallest number of parameters needed) and good fit. A lower AIC indicates a better model.

step(severity.mult.lm)

Output

step(severity.mult.lm)
Start:  AIC= 11.83 
 disease ~ temperature + rel.humidity + rainfall 

               Df Sum of Sq    RSS    AIC
- rainfall      1     0.804 15.467 10.361
                      14.662 11.827
- rel.humidity  1    17.326 31.988 17.628
- temperature   1    28.189 42.852 20.552

Step:  AIC= 10.36 
 disease ~ temperature + rel.humidity 

               Df Sum of Sq    RSS    AIC
                      15.467 10.361
- rel.humidity  1    17.425 32.892 15.906
- temperature   1    33.031 48.497 19.789

Call:
lm(formula = disease ~ temperature + rel.humidity, data = severity.mult)

Coefficients:
 (Intercept)   temperature  rel.humidity  
    -78.2751        0.8308        0.5438  

This output indicates the steps taken to select the model. The first model tested includes all three predictor variables and the AIC for that model is calculated. Then the AIC for each model that excludes a single predictor variable is calculated. Excluding rainfall results in a lower AIC, so that model is selected in the next step. When either temperature or relative humidity are removed from that model, the AIC goes up, indicating that the best model includes temperature and relative humidity as predictors and excludes rainfall. The new model still has the issue of multicollinearity between the two predictors (temperature and relative humidity). Further considering may be given for other methods to deal with this effect, including Ridge or Bayesian regression models (Neter et al. 1996). In the next step illustrated below, this model is analyzed.

severity.mult.lm2 <- lm(disease~temperature+rel.humidity,data=severity.mult)
summary(severity.mult.lm2)

Output

severity.mult.lm2 <- lm(disease~temperature+rel.humidity,data=severity.mult)
summary(severity.mult.lm2)

Call:
lm(formula = disease ~ temperature + rel.humidity, data = severity.mult)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.2463 -0.5882 -0.2936  0.2904  3.2503 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)   
(Intercept)  -78.2751    23.9119  -3.273  0.01361 * 
temperature    0.8308     0.2149   3.866  0.00616 **
rel.humidity   0.5438     0.1936   2.808  0.02621 * 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 1.486 on 7 degrees of freedom
Multiple R-Squared: 0.8328,     Adjusted R-squared: 0.7851 
F-statistic: 17.44 on 2 and 7 DF,  p-value: 0.00191 

Just as for simple linear regression (see Sparks et al., 2008), it’s important to check that the model assumptions are adequately met. The following code produces four plots to help in model evaluation. The plot function in R provides several types of diagnostic graphs that can be selected by changing the number called in which = _.

par(mfcol=c(2,2)) #Sets the graphic to display all four plots in one window
plot(severity.mult.lm2,which=1) #residuals versus fitted
plot(severity.mult.lm2,which=2) #QQ plot
plot(severity.mult.lm2,which=3) #plot of square-root standardized residuals #against fitted value
plot(severity.mult.lm2,which=5) #Residuals versus leverage plot

The four plots are of: (i) Residuals vs. fitted (or predicted) values, (ii) Scale-location (square-root of the standardized residual) vs. fitted values, (iii) Normal Quantile-Quantile (Q-Q) relation, and (iv) Residuals vs. leverage (influence of a single observation on model estimates).

Output

Click to enlarge.

Ideally the first two plots will show no relationship. A straight line in the Q-Q plot indicates a normal distribution for the residuals. Ideally the leverage of each observation would be equal. If Cook's distance is > 1 the observation may be considered an outlier. To obtain a summary of influential cases, using measures such as DFFITS, DFBETAS, and Cook's distance, there is a function called influence.measures that can be used with lm or glm objects (Neter et al. 1996).

The model can also be used to predict disease severity for new values of temperature and relative humidity. The following code produces the prediction for temperature = 95 and relative humidity = 30, along with confidence limits for the prediction.

new <- data.frame(temperature=95,rel.humidity=30) #new temperature and new rel humidity
predict(severity.mult.lm2,newdata=new,interval="confidence")

Output

          fit      lwr      upr
[1,] 16.96165 11.42027 22.50303

Application of maximum likelihood based methods for disease forecasting

With the development of more rapid computing technologies, it is possible to compute a large number of simulations to support new modeling methods for plant disease forecasting. The development of disease forecasting models is based on determining the value of unknown parameters for an appropriate model.

One method for estimation of parameters is maximum likelihood estimation (MLE), first developed by R.A. Fisher (Myung 2003). The basic idea behind MLE is that one can calculate the probability of observing a particular set of data given the model of interest and this is called the likelihood of a sample. The likelihood of observing the actual data for each of the potential parameter values is calculated to determine which set of parameter values maximizes the likelihood.

To illustrate the concept of MLE, let's examine a random sample from a binomial distribution with given parameter p. MLE is then used to determine the parameter estimate, as if there was no knowledge of the parameter. Finally, the parameter estimate based on MLE is compared to the true value.

Before getting into this example, some background on MLE is necessary. To begin, based on some prior knowledge, a model for the data may be developed, which may be defined as a probability distribution function, f. The likelihood, L, is defined as a function of the parameter(s) in the model given the data: L(θ|x). When the number of samples is > 1, the likelihood is a joint likelihood of observations. In practice, it is often easier computationally to work with the log-likelihood, defined as logL(θ|x). The logL is maximized by taking the first derivative with respect to the parameter in order to obtain the estimate of the parameter.

Consider a binomial sampling situation where a coin is tossed 10 times (in R: rbinom(10,1,0.5)). If using 1 to represent the heads and 0 to the tails of the coin, this is one possible outcome for those 10 events: 0 1 0 1 1 0 0 1 0 0.

Given this data, a natural question is: What is the probability that a head will be observed on the next toss? One approach may be to count the number of times that a head was observed in the initial probability experiment. Following this approach, given an observation of four heads out of 10 coin tosses, the probability of a head on the next toss may be 0.4. This probability may also be estimated using MLE as a simple example of the MLE method, assuming that the occurrence of heads follows a binomial distribution with a probability parameter p. Formally, Prob(x = 1) = p and Prob(x = 0) = 1 - p, where x is an observation. All 10 tosses of the coin are considered independent, such that the 10 observations x1, x2,...,x10 can be considered a random sample from a binomial distribution. The joint likelihood becomes:

L(x1, x2, ..., x3 | p)

= Prob(x1 = 0 | p) * Prob(x2 = 1 | p) *...* Prob(x10 = 0| p)

= (1 - p) * p * (1 - p) * p * p * (1 - p) * (1 - p) * p * (1 - p) * (1 - p)

= (1 - p)⁶ p⁴

and logL as

logL = 6log(1 - p) + 4 log( p)

Now, the goal is to maximize the logL function. This may be accomplished by taking the first derivative of logL with respect to p and setting it equal to zero. In this situation, an explicit solution for p is obtained as:

solving this equation, gives p = 0.4.

This implies that the probability of heads is 0.4. But how good is this estimate? For the example above, the observations were obtained from a random binomial distribution with probability p = 0.5. Does this imply that the MLE estimate is not useful? In this case the random sample was only of size 10, which is a fairly small sample. Sample size is an important consideration for use of MLE or other estimation methods for plant disease forecasting models. The smaller the number of observations, the bigger the difference between the estimate and the true parameters is likely to be.

      1 0 1 1 1 0 0 0 0 1 1 0 0 1 0 0 1 0

      1 1 1 0 1 0 0 0 0 1 1 1 0 1 0 0 1 0

      1 1 0 0 1 0 0 1 1 0 0 1 0 1

#The glm function is available in the Stats library.  
# For more information on glm, type: ?glm

#Example data where the event is presence or absence of disease (prevalence)
#(1=presence,0=absence), with associated air temperature over 3-months, a measure of maximum
#snowfall, and an indicator for whether the disease occurred the previous year.  We might assume
#that favorable winter weather increases the survival of a pathogen.  This
#information can be used to develop a model using logistic regression, which
#uses maximum likelihood to obtain the estimates of a parameter

prevalence=c(1,1,0,0,0,0,0,1,0,0,0,0,1,1,1,0,1,1,0,1,1,0,1,1,1,1,1,0,1,1,0,1,0,1,0,0,1,0,0,
1,1,0,0,0,1,1,1,0,0,0,0,1,1,1,0,1,0,1,1,0,1,0,1,1,0,1,1,0,1,0,0,0,1,1,1)

month1=c(-6.15,-4.78,-5.47,-4.79,-4.94,-4.68,-5.05,-4.58,
-4.53,-4.78,-4.85,-5.31,-4.95,-4.18,-3.95,-4.53,
-4.12,-3.91,-5.25,-5.78,-5.18,-5.22,-4.48,-5.87,
-4.63,-4.32,-4.40,-4.16,-4.22,-2.89,-4.60,-5.70,
-5.72,-4.27,-4.24,-5.04,-6.39,-4.73,-5.69,-3.85,
-4.86,-4.55,-3.43,-5.58,-5.09,-4.40,-4.15,-3.65,
-4.56,-5.43,-3.23,-5.86,-5.47,-5.55,-5.31,-3.86,
-5.48,-4.44,-6.27,-4.84,-4.33,-5.07,-6.51,-5.51,
-4.79,-5.93,-4.04,-5.42,-5.16,-5.12,-5.01,-3.79,
-5.93,-6.43,-4.76)

month2=c(-8.52,-7.14,-7.99,-8.53,-6.76,-6.71,-8.27,
-8.08,-9.29,-6.03,-5.43,-6.36,-8.71,-7.57,
-8.13,-6.99,-8.40,-6.08,-9.01,-5.78,-7.48,
-7.07,-6.88,-6.87,-7.17,-9.07,-4.66,-10.05,
-8.68,-7.73,-9.49,-6.85,-6.51,-7.47,-7.96,
-10.65,-7.62,-7.14,-8.38,-7.01,-7.02,-6.65,
-7.34,-7.50,-8.90,-8.80,-9.12,-8.13,-8.79,
-8.09,-5.58,-5.46,-7.72,-8.22,-6.96,-5.55,
-8.60,-8.88,-8.47,-8.55,-7.89,-9.25,-7.44,
-9.06,-7.50,-8.05,-7.52,-7.85,-9.33,-8.72,
-7.12,-9.40,-9.76,-8.41,-8.79)

month3=c(-5.09,-3.84,-5.24,-3.89,-5.22,-5.60,-5.14,-2.62,
 -2.44,-4.86,-5.46,-3.20,-5.35,-4.04,-3.78,-5.03,
-4.66,-5.46,-3.59,-3.80,-4.10,-3.17,-2.72,-4.32,
-5.69,-4.23,-3.83,-6.20,-6.12,-2.05,-4.41,-4.75,
-4.13,-4.95,-4.33,-5.83,-4.73,-4.66,-3.41,-2.53,
-5.41,-5.39,-5.87,-6.87,-4.99,-5.27,-5.88,-5.16,
-4.72,-3.80,-2.98,-4.23,-4.69,-2.31,-5.03,-3.40,
-2.79,-2.81,-3.98,-5.13,-4.77,-3.85,-3.97,-2.70,
-4.21,-2.82,-5.47,-3.57,-4.53,-4.39,-3.47,-4.35,
-3.56,-5.42,-3.38)

maximum.snow=c(163,201,213,183,169,197,136,180,137,
168,189,168,172,238,170,183,212,157,
184,245,247,148,151,166,153,212,197,
202,227,153,237,195,208,223,174,219,
217,128,119,177,111,175,136,186,200,
197,159,227,246,194,247,222,160,247,
172,217,194,222,168,243,157,229,149,
185,211,178,225,193,201,202,212,154,
196,272,152)

previous.prevalence=c(0,1,0,0,0,0,0,1,0,0,0,0,0,
1,1,0,0,1,0,1,0,0,1,0,0,0,
1,0,0,1,0,0,0,0,0,0,0,0,0,
1,0,0,0,0,0,0,0,0,0,0,0,1,
0,0,0,1,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0)

#Create a dataframe
example.forecasting=as.data.frame(cbind(prevalence,month1,month2,month3,maximum.snow,previous.prevalence))

#Begin by considering a model with the air temperature for month1 as the only predictor of prevalence

model1=glm(prevalence~month1,family=binomial,data=example.forecasting)
summary(model1)

Call:
glm(formula = prevalence ~ month1, family = binomial, data = example.forecasting)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.301  -1.230   1.010   1.135   1.275  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  -0.7501     1.5077  -0.498    0.619
month1       -0.1812     0.3057  -0.593    0.553

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 103.64  on 74  degrees of freedom
Residual deviance: 103.28  on 73  degrees of freedom
AIC: 107.28

Number of Fisher Scoring iterations: 3

#To determine if month1 is a useful predictor, we can perform a deviance test,
#which compares the residual deviance with the parameter added to the null deviance.
#Here, this is 103.64-103.28 with 1 df and is a chi-square test.

1-pchisq(103.64-103.28,1)

1-pchisq(103.64-103.28,1)
[1] 0.5485062

#We can pick a particular model, but also may use the step function to
#compare models and select the model with the best fit:

step(glm(prevalence~1,data=example.forecasting,family=binomial),scope=list(lower=~1,upper=~month1*month2*month3*maximum.snow*previous.prevalence),direction='forward')

Start:  AIC= 105.64 
 prevalence ~ 1 

                      Df Deviance     AIC
+ previous.prevalence  1   86.557  90.557
                    103.638 105.638
+ month3               1  102.206 106.206
+ month1               1  103.285 107.285
+ maximum.snow         1  103.413 107.413
+ month2               1  103.530 107.530

Step:  AIC= 90.56 
 prevalence ~ previous.prevalence 

               Df Deviance    AIC
+ month1        1   82.930 88.930
              86.557 90.557
+ month2        1   84.622 90.622
+ maximum.snow  1   86.394 92.394
+ month3        1   86.557 92.557

Step:  AIC= 88.93 
 prevalence ~ previous.prevalence + month1 

                             Df Deviance    AIC
+ month2                      1   80.784 88.784
                            82.930 88.930
+ maximum.snow                1   82.659 90.659
+ month3                      1   82.688 90.688
+ month1:previous.prevalence  1   82.930 90.930

Step:  AIC= 88.78 
 prevalence ~ previous.prevalence + month1 + month2 

                             Df Deviance    AIC
                            80.784 88.784
+ month3                      1   80.344 90.344
+ maximum.snow                1   80.726 90.726
+ month1:month2               1   80.783 90.783
+ month1:previous.prevalence  1   80.784 90.784
+ month2:previous.prevalence  1   80.784 90.784

Call:  glm(formula = prevalence ~ previous.prevalence + month1 + month2,      family = binomial, data = example.forecasting) 

Coefficients:
        (Intercept)  previous.prevalence               month1  
            -6.9184              19.8306              -0.7426  
             month2  
            -0.3709  

Degrees of Freedom: 74 Total (i.e. Null);  71 Residual
Null Deviance:      103.6 
Residual Deviance: 80.78        AIC: 88.78 

step(glm(prevalence~month1*month2*month3*maximum.snow*previous.prevalence,data=example.forecasting,family=binomial),scope=list(lower=~1,upper=~month1*month2*month3*maximum.snow*previous.prevalence),direction='backward')