Chapter 5

Regression Analysis, Frequentist

Chapter Outline
1 Introduction 161 A.2 Nelder-Mead Algorithm 188
2 Ordinary Least Squares (OLS) 164 A.3 Comparison of Algorithms 190
2.1 Estimation by OLS 165 A.4 Potential Problems 190
2.2 The Properties of OLS 167 B R Functions for Parameter
3 Maximum Likelihood 170 Estimation 191
4 What if the Assumptions B.1 The nls Function When the
of Equations (3)–(5) Model Consists of a Single
Are Violated? 171 Equation 191
4.1 Assumption 1: The Expectation B.2 nls for More Complex
of Model Error Is 0 for All X 171 Models 196
4.2 Assumption 2: Constant B.3 The optim Function for
Variance of Model Error 176 Function Minimization 199
4.3 Assumption 3. All Residual Exercises 202
Errors Are Independent 180 References 204
A Algorithms for OLS 186
A.1 Gauss-Newton Algorithm 186

1 INTRODUCTION
A regression equation is an equation that relates some dependent or response
variable y to some independent or explanatory variables X ¼ x1, …, xn, and
involves some unknown parameters θ ¼ θ1, …, θp. Regression analysis, which
involves estimating the unknown parameters and the quality of prediction of the
regression function, is a central topic in statistics and there is a large amount of
literature on the subject (for nonlinear regression see, e.g., Bates and Watts
(1988), Rhinehart (2016), or Seber and Wild (1989)). System models are in fact
regression models, so it is important to understand how such regression models
are treated in statistics.
In this chapter we discuss two frequentist approaches to parameter estima-
tion for regression models: least squares and maximum likelihood. We discuss
in detail the usual assumptions that the regression model is correctly specified
(defined rigorously below), and that all model errors have the same variance
and are independent. These assumptions are often not verified by system
models, so we discuss methods of dealing with situations where the

Working with Dynamic Crop Models. [Link]

© 2019 Elsevier B.V. All rights reserved. 161
162 SECTION A Background

assumptions are incorrect. Two appendices give more information on algo-

rithms for nonlinear least squares, and R functions for implementing these algo-
rithms. Chapter 7 goes into more detail about frequentist parameter estimation
specifically for system models, often referred to as “calibration.” Chapter 8 goes
into more detail about Bayesian parameter estimation.
The general form of a regression equation is
y ¼ f ðX; θÞ + ε (1)
where f(X; θ) is a function of X and θ and ε is a random variable, the model error.
We assume that we have data from a sample of n individuals from the popula-
tion of interest (the “target” population), which is described by Equation (1).
The response of each individual i in the sample is then
yi ¼ f ð X i ; θ Þ + ε i (2)
Equations (1), (2) simply say that there is some difference between the model
and the observed value, that we note ε. It contains no useful information until we
make some assumptions about the model error. The standard simple assump-
tions are shown in Table 5.1.
The assumption of Equation (3) says that for some parameter vector
θ ¼ θtrue, model error has expectation 0 for all X. Taking the expectation on both
sides of Equation (2), we see that this implies that

TABLE 5.1 Standard Assumptions About Model Error

For some θ ¼ θtrue

Assumption 1. Expected model error equals 0. E(εi) ¼ 0 (3)

Assumption 2. Constant variance of model error. var(εi) ¼ σ 2
(4)
Assumption 3. Independent model errors. All εi independent (5)

Eðyj XÞ ¼ f ðX; θtrue Þ (6)

where the notation y | X means y given that X so the expectation is just over y and
not over X. Equation (3) or equivalently Equation (6) define what we mean when
we say that a model is “correct.” We mean that there exists some parameter vector
θ ¼ θtrue such that, if we average the measured values over all individuals in the
target population for any specific value of X, the result is equal to f(X; θtrue).
This also defines what is meant by “true” parameter values. They are the param-
eter values such that Equation (6) holds, that is, such that, for every X, the
expectation of the measurement is equal to the value calculated using the model.
Note that a “correct” model is not necessarily a model that is a good predic-
tor. Model error has expectation zero, but for each individual it will not in
general be zero. The mean squared error of prediction for the correct model
is equal to the variance of model error, which can be small or large.
 Regression Analysis, Frequentist Chapter 5 163

The second assumption, in Equation (4), says that the variance of model
error is the same for all the Xi. Finally, Equation (5) says that all the model errors
are independent. That is, knowing the model error for one individual in the sam-
ple gives us no information about model error for another individual. If each
individual in the sample is chosen at random in the population, then clearly
the model error for individual i0 6¼ i will have no relation to the model error
for individual i. Thus, this assumption is automatically fulfilled if we have a
simple random sample.
To illustrate many of the notions discussed in this chapter, we will use the
data and model given in Example 1. Example 2 illustrates what a large sample
from the target population of Example 1 would look like, if the three assump-
tions of Equations (3)–(5) are satisfied.

EXAMPLE 1 Wheat Grain Weight Versus Degree Days After Anthesis

Darroch and Baker (1990) studied grain filling in three spring wheat genotypes in a
single field over three years. The data for cultivar Neepawa in 1988 are shown in
Figure 5.1. They found that a logistic function fit their data quite well, implying that
a good choice for f(X; θ) is:
W
f ðX; θÞ ¼ (7)
1 + e BC ∗DD
This model has a single input variable, degree days after anthesis, noted DD, and
three parameters, noted W, B, and C.

25
Grain weight (mg)

20

15

10

200 300 400 500 600

DD after anthesis
FIGURE 5.1 Data of grain weight versus degree days after anthesis.
164 SECTION A Background

EXAMPLE 2 A Large Sample Simulated on the Basis of the Three Assumptions

in Table 5.1
Consider a population such that
W true
yi ¼ + εi (8)
1 + exp B C ∗DDi
true true

with (Wtrue, Btrue, Ctrue) ¼ (28.2, 4.16, 0.0165) and


εi  N 0, σ 2 (9)
iid

with σ 2 ¼ (0.84)2 ¼ 0.706. This implies that the model in Equation (7) is a correct
model because there are parameter values such that the expectation of ε is 0 for
all values of the explanatory variable. The assumption (Equation 4) of constant var-
iance is also satisfied, and all the model errors are independent (assumption in
Equation 5).
Figure 5.2 shows what a large sample from this target population would look like.
It is assumed that the sampling dates (values of DD) are fixed. As can be seen, the
model with the true parameter values goes through the middle of the points for each
value of DD (assumption 1). The spread around the model is the same for all values of
DD (assumption 2). By construction the model errors are independent, though this
cannot be seen from the figure.

30

25
Grain weight (mg)

20

15

10

200 300 400 500 600

DD after anthesis
FIGURE 5.2 Large sample of measured values in the case where the three assumptions
about residual error in Equations (3)–(5) are satisfied. The solid line is the model of
Equation (7) with (W,B,C) ¼ (Wtrue,Btrue,Ctrue).

2 ORDINARY LEAST SQUARES (OLS)

According to this approach, the parameters are chosen to minimize the sum of
squared errors (SSE) between simulated and measured values:
 Regression Analysis, Frequentist Chapter 5 165

X
n
SSE ¼ ½yi  f ðXi ; θÞ2 (10)
i¼1

where n is the number of data points. We can say the same thing with the fol-
lowing equation:
( )
X
n
θ^OLS ¼ arg min ½ y i  f ðXi ; θ Þ 2
(11)
θ i¼1

which says that the OLS parameters are the arguments of the model that min-
imize the sum of squared errors. This is called “ordinary” least squares because
it is a simple sum of squared errors, not a weighted sum, that is to be minimized.
When we discuss how to deal with violations of the assumptions in
Equations (4), (5), we will present other versions of the least squares approach,
namely weighted least squares (WLS) and generalized least squares (GLS).
Suppose that the regression model is linear in the parameters, so that it can
be written
y ¼ θ0 + θ1 x1 + ⋯ + θp xp + e (12)
Let Y ¼ y1, …, yn be the column vector of observed values. The explanatory var-
iables of the sample can be put in matrix form
0 1
1 x11 … x1p
X ¼ @⋮ ⋮ ⋱ ⋮ A
1 xn1 … xnp
where xij is the value of the jth explanatory variable for individual i. The first
column with 1s multiplies the intercept term θ0 in the regression equation. The
regression equation in matrix notation is then
Y ¼ Xθ + ε
where θ is the column vector of parameters and ε ¼ ε1, …, εn is the column
vector of model errors.

2.1 Estimation by OLS

For a linear regression function, there is an analytical expression for the
ordinary least squares estimator of the parameter values, namely

1
^
θOLS ¼ XT X XT Y
An unbiased estimator of the variance of model error is
1
T

σ^2OLS ¼ θOLS Y  X^θOLS
Y  X^
np
The variance-covariance matrix of the parameter estimator

  1
var ^θOLS ¼ XT X σ^2OLS (13)
166 SECTION A Background

If the regression equation is a nonlinear function of the parameters, there is

no general analytical solution to Equation (11) and the solution must be calcu-
lated using an algorithm that approaches the solution in successive iterations.
R functions that do these calculations are presented in the appendix. Given
^θOLS , the estimator of σ 2 for nonlinear regression is
( )
Xn 
2
σ^ ¼
OLS
2
yi  f Xi ; θ^OLS =ðn  pÞ ¼ SSE=ðn  pÞ (14)
i¼1

where n is the sample size and p is the number of parameters.

The variance-covariance matrix of ^ θOLS in nonlinear regression can be
approximated based on Equation (13), replacing X by derivatives of the regres-
sion function with respect to the parameters, evaluated at ^θOLS . Thus


1
var ^θOLS  FT F σ 2 (15)
where
0   1
∂f ðx1 ; θÞ ∂f ðx1 ; θÞ
B ∂θ  true … ∂θp θtrue C
B 1 θ C
B
F¼B ⋮  ⋱ ⋮  C C
@ ∂f ðxn ; θÞ ∂f ðxn ; θÞ A
…
∂θ1 θtrue ∂θp θtrue
To obtain a practical estimator, use σ^2OLS in place of σ 2 and ^θOLS in place of θtrue.
The result of OLS parameter estimation for Example 1 is shown in Example 3.
Given the estimated parameters, we can plug them into the model in order to
make predictions for
any new value of the input variables. To predict for, say,
X∗, we would use f X∗ ; ^θOLS .

EXAMPLE 3 OLS Parameter Estimation for the Model and Data of Example 1
The R function nls was used to estimate the parameters of the logistic model of
Example 1, based on the data presented there. The OLS parameter values are
W^ OLS ¼ 28:17, B^OLS ¼ 4:160, and C^ OLS ¼ 0:01649. The estimated variance of model
error is σ^2 ¼ 0:8364. The variance-covariance matrix of θ^OLS is
OLS

W B C
W 0.1666 0.06338 3.0831e04
B 0.06338 0.1297 4.979e04
C 3.0831e-04 4.979e04 2.002e06

The fit of the model to the data and the model residuals (observed minus simu-
lated values) are shown in Figure 5.3.
 Regression Analysis, Frequentist Chapter 5 167

Seed weight 25

20

15

10

200 300 400 500 600

(A) Growing DD since anthesis

0.5
Residulas

0.0

−0.5

−1.0

200 300 400 500 600

(B) Growing DD since anthesis
FIGURE 5.3 (A) Data (dots) and model (line) for logistic model with parameters estimated
using OLS. (B) Residuals.

2.2 The Properties of OLS

If the assumptions of Equations (3)–(5) are satisfied, then in the case of the linear
model, ^ θOLS is the best possible linear unbiased estimator, meaning that among all
possible estimators that are unbiased and are linear combinations of the yi values
θOLS ), ^
(as is ^ θOLS has the smallest possible mean squared error as an estimator of
θ (Graybill, 1976). Also, σ^2OLS is an unbiased estimator of σ 2.
true
168 SECTION A Background

For a nonlinear model, OLS has similar properties but only asymptotically,
that is, as the size of the sample tends toward infinity (Seber and Wild, 1989,
p. 565). There are also additional constraints on the model, essentially that it
must be well behaved. The difference ^ θ  θtrue has a normal distribution asymp-
totically, with

1 
^
θ  θtrue  N 0, σ 2 F0 F (16)

One should neither be too impressed nor too unimpressed by the fact that these
properties only hold asymptotically for nonlinear models. One should not be too
impressed because samples are always finite in size, often quite small, and we
have promised no good qualities for OLS for nonlinear models for finite sample
size. On the other hand, asymptotic properties are far from useless. Imagine an
estimator that did not promise to converge to the true parameter values at least
asymptotically. We would have serious doubts about using such an estimator.
Example 4 illustrates some of the above properties using simulated data based
on the population described in Example 2. We drew repeated samples from the
population, and for each sample calculated the OLS estimates of the parameter
values and of model error variance. In this artificial example we know the true
values θtrue and σ 2, and so can compare the OLS estimates with the true values.
Because we know exactly the distribution of y in the population, we can also cal-
culate the mean squared error of prediction for any value of ^θOLS .

EXAMPLE 4 Effect of Sample Size on OLS Parameter Estimator

We drew 500 different samples of size 11 (one value at each value of DD), or of size
33 (three values at each DD) or of size 110 (10 values at each DD) from the popu-
lation described in Example 2. The estimates ^ θOLS and σ^2OLS were calculated for each
sample using the R function nls. Samples where nls failed to converge were
eliminated.
The first three lines of Table 5.2 show the mean (averaged over 500 samples) of the
parameter estimates. As sample size increases from 11 to 33 to 110, the means approach
closely θtrue, as expected because the OLS estimator is asymptotically unbiased. In fact,
they are already quite close to θtrue even for the smallest sample size. The next three lines
show that the standard deviations of the estimators diminish substantially with larger
sample size. This is an expression of the fact that the estimator is consistent; as sample
size increases, the estimated values cluster more and more closely around the true
value. The line “mean(^ σ 2OLS )” shows the mean over samples of estimated variance of
model error, which is quite close to the true value for all sample sizes.
The last line shows the average of mean squared error of prediction (MSEP). The
minimum possible value, attained for θ ¼ θtrue, is MSEP ¼ σ 2 ¼ 0.706. For small sam-
ple size, average MSEP is substantially larger because the estimated parameter
values are often quite far from the true values. As the sample size increases, there
is less spread of the estimated parameter values around the true value, and so aver-
age MSEP diminishes.
 Regression Analysis, Frequentist Chapter 5 169

TABLE 5.2 Effect of Increasing Sample Size on OLS Estimators

11 33 110
(1 Measurement (3 Measurements (10 Measurements
Sample Size at Each x) at Each x) at Each x)
Mean(W ^ OLS ) 28.28 28.20 28.20
Wtrue ¼ 28.2
Mean(B^OLS ) 4.158 4.148 4.157
Btrue ¼ 4.16

Mean(C^ OLS ) 0.01650 0.01646 0.01649

Ctrue ¼ 0.0165
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

ffi
^ OLS 0.394 0.232 0.131
var W
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

 0.369 0.210 0.111
var B^OLS
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
  0.00144 0.00082 0.00044
var C^ OLS


Mean σ^2OLS 0.711 0.706 0.706
σ ¼ 0.706
2

Mean(MSEP) 0.873 0.766 0.724

minimum
MSEP ¼ 0.706

Figure 5.4 shows histograms of B^OLS for the different sample sizes. This again
shows that as sample size increases, there is less and less variability in the estimated
value around Btrue.

Sample size = 11 Sample size = 33 Sample size = 110

100
100
80 150
80
Frequency

60
Frequency

100
Frequency

60
40
40
50
20 20

0 0 0
3.0 4.0 5.0 6.0 3.0 4.0 5.0 6.0 3.0 4.0 5.0 6.0
B B B

FIGURE 5.4 Histograms of the estimated values of parameter B in different samples for
increasing sample size.
170 SECTION A Background

3 MAXIMUM LIKELIHOOD
There is a second major way of obtaining estimates of parameter values other
than least squares, called maximum likelihood (ML). This is an extremely gen-
eral procedure that can be used whenever one can make reasonable assumptions
about the distribution of model error. (Note that for least squares, we did not
need to make assumptions about the distribution of model error, other than con-
stant variance.)
The likelihood that θtrue ¼ θ for some parameter vector θ, written L(θ j Y), is
defined as the probability that Y would equal the observed Y, if θtrue ¼ θ. The ML
estimator of θtrue is then the value of θ that maximizes L(θ j Y). That is,
the estimator is the parameter vector that maximizes the probability of obtaining
the values that were in fact measured.
Consider the simple case, where the model errors are all independent and
normally distributed with the same variance:
yi ¼ f ð
Xi ; θtrue
 Þ + εi
εi  N 0, σ 2 (17)
iid

The assumptions about model error here are the same as in Equations (3)–(5),
but with the additional assumption of a normal distribution.
According to Equation (17), if θtrue ¼ θ, then the model errors are εi ¼ yi 
f(Xi; θ). The probability of obtaining those error values is
Y
n
1
eðyi f ðXi , θÞÞ =ð2σ Þ
2 2
Probðε1 , …, εn Þ ¼ Lðθj Y Þ ¼
i¼1 ð2πσ 2 Þ1=2
X
1  ðyi f ðXi , θÞÞ2 =ð2σ 2 Þ
¼ e (18)
ð2πσ 2 Þn=2
where we have used the fact that for independent variables, the joint probability
is the product of the individual probabilities.
The ML estimator is the value of θ that maximizes the above expression. The
only part of the expression that involves θ is the exponent, and the expression is
maximized by minimizing
X
n
SSE ¼ ðyi  f ðXi , θÞÞ2
i¼1

This is exactly the same expression as for the OLS estimator. In this simple case,
the OLS estimator and the ML estimator are identical.
It is also possible to show that the value of σ 2 that minimizes Equation (18) is
X
n
σ^2ML ¼ ðyi  f ðXi , θÞÞ2 =n
i¼1
 Regression Analysis, Frequentist Chapter 5 171

It can be shown that this is a biased estimator, so it is usually replaced by the

unbiased estimator
X
n
σ^2ML, unbiased ¼ ðyi  f ðXi , θÞÞ2 =ðn  1Þ (19)
i¼1

This again is identical to the OLS estimator.

If the assumptions in Equation (17) are valid, the ML estimator is asymp-
totically consistent and asymptotically has minimum variance for an unbiased
estimator (Casella and Berger, 1990).

4 WHAT IF THE ASSUMPTIONS OF EQUATIONS (3)–(5)

ARE VIOLATED?
4.1 Assumption 1: The Expectation of Model Error Is 0 for All X
4.1.1 Checking Assumption 1
Assumption 1 is that there exist parameter values, θtrue, called the “true” param-
eter values, such that for each value of the explanatory variables X, the model is
equal to the expectation of the response. In essence, this says that the model has
the correct dependence on X. If this assumption is satisfied, we say that the
model is “correctly specified” or simply “correct.” If the assumption is not sat-
isfied, we say that the model is “misspecified.”
For fairly simple models
with few explanatory variables, a simple check is to
plot residual error ei ¼ yi  f Xi ; ^θOLS versus each explanatory variable to see if
the residuals seem to be centered at 0 for all X. As an example, consider two
different models fit to the data of Example 1. Example 3 showed the fit and
a plot of residuals for a logistic model. Example 5 shows the fit and correspond-
ing residuals for an exponential model. While the residuals of the logistic model
have no obvious structure and seem to be centered around 0, the residuals of the
exponential model do have structure. The residuals are systematically negative
for small and large values of DD and positive in between. In this case, it looks
like we could propose a model to relate residual error to DD. This implies that
the model f(X; θ) does not include all the dependence of y on DD because there is
still some structured dependence left in the residual error. It is this structure of
the residuals that indicates that the model does not correctly take into account
the effect of the input variable.

EXAMPLE 5 A Misspecified Model Fit to the Data of Example 1

The two parameters a and b of the exponential model f(x; θ) ¼ a(1  e b∗DD) were fit
to the data of Example 1 using OLS. The calculations were done using the R function
nls. The left graph of Figure 5.5 shows the fit of the model to the data. The right
graph shows the residuals.
There is some clear dependence on DD of the residuals. That is, the model does
not correctly describe the effect of DD on the response. The model is misspecified.

Continued
172 SECTION A Background

EXAMPLE 5 A Misspecified Model Fit to the Data of Example 1—cont’d

25
Seed weight

20

15

10

200 300 400 500 600

(A) Growing DD since anthesis

2
Residulas

−2

−4

200 300 400 500 600

(B) Growing DD since anthesis
FIGURE 5.5 (A) Misspecified model. Data (dots) and model (line). (B) Misspecified model.
Plot of residuals.

4.1.2 Consequences if Assumption 1 Is Not Satisfied

White (1981) studied the consequences of using OLS to estimate model parameters
for misspecified models. The major assumptions he made are that the input variables
are random variables and model errors are all independent and identically distrib-
uted. The results concern the asymptotic properties of OLS estimation.
A first conclusion concerns the values of the estimated parameters. We
define here “optimal” parameter values, which are the parameter values that
 Regression Analysis, Frequentist Chapter 5 173

minimize mean squared error of prediction (MSEP ) for the misspecified model,
asymptotically. White (1981) showed that in general the optimal parameter
values are not equal to the true parameter values. In fact, because the misspe-
cified model need not have the same functional form as the correct model, it
may have parameters that are defined completely differently than the correct
model. However, even if the misspecified model has parameters that are anal-
ogous to those of the correct model, the estimated values do not tend toward the
true parameter values. This is because the optimal parameter values are com-
pensating for errors in model formulation.
A second conclusion concerns MSEP for a misspecified model. White (1981)
showed that for a misspecified model, the OLS estimators tend toward the optimal
parameters that minimize MSEP as the amount of data tends toward infinity. On
the one hand, this is perhaps not so surprising. OLS minimizes the sum of squared
errors. If we do this for larger and larger samples, it is not surprising that we
approach the parameters that minimize mean squared error for the entire popula-
tion. On the other hand, this is a very powerful result. It says that even if the model
is misspecified, the OLS parameters will make the model into the best it can be.
Example 6 illustrates the above two effects of misspecification with an artificial
example. Two different misspecified models are shown. In the first, the functional
dependence on the explanatory variable x is wrong. The asymptotic parameter
values in this case are not comparable to the true parameter values; in the true
response, the parameters describe the exponential increase of y while in the model
the parameters are the coefficients of a linear and quadratic term. In the second
model, the dependence on x has the correct form, but one of the model parameters
is fixed at an incorrect value. This second model is also misspecified because there
is no value of the single estimated parameter that will make model error 0 for all x.
The estimator of the single parameter does not tend asymptotically to the true
value; it tends toward a value that compensates for the error in the fixed parameter.

EXAMPLE 6 Simple Examples of Misspecified Models

Suppose that the true response is

yi ¼ 0:09

exp ð0:4x
 Þ + εi (20)
ei  N 0, 0:12
The distribution of x values in the target population is assumed to be x U(5, 15) (the
x values are drawn at random from a uniform distribution with lower and upper
limits of 5 and 15, respectively). All the errors and explanatory variables are mutually
independent. The correctly specified model is then

hðx; θÞ ¼ θ1 exp ðθ2 x Þ

with true parameter values


1 , θ2
θtrue true
¼ ð0:09, 0:4Þ
We consider two different misspecified models:

h1 ðx; θÞ ¼ θ1 + θ2 x 2
h2 ðx; θÞ ¼ θ1 exp ð0:3x Þ

Continued
174 SECTION A Background

EXAMPLE 6 Simple Examples of Misspecified Models—cont’d

The first model has the wrong functional form (quadratic rather than exponential). The
second has the correct functional form, but one of the parameters is set to an incorrect
value. The parameters for these two models were estimated using a sample of size
30,000, which should be a good approximation to the asymptotic situation. The
results, shown in Table 5.3, show that the estimated parameter values are not equal
to the true parameter values. The fitted models are shown in Figure 5.6. OLS parameter
estimation has made the models as close as possible to the true response.

TABLE 5.3 OLS Parameters Estimated for Misspecified Models

Model ^
θOLS,1(θtrue50.09) ^
θOLS,2(θtrue50.4)
h1(x; θ) (wrong function) 7.298 0.1498
h2(x; θ) (wrong coefficient) 0.3464 Fixed at 0.3

Wrong function Wrong coefficient

True True
35 Model 35 Model

30 30

25 25

20 20
y

15 15

10 10

5 5

0 0

6 8 10 12 14 6 8 10 12 14
x x
FIGURE 5.6 Two misspecified models showing the true expectation of the response (solid line)
and simulated values using parameter values estimated by OLS on a sample of size 30,000 (dashed
line).
 Regression Analysis, Frequentist Chapter 5 175

A third important consequence of violation of assumption 1 concerns cali-

bration for different target populations. Suppose that we have two different tar-
get populations with the same correct model and the same true parameter
values. This might be, for example, corn fields in two different regions. For
the correct model, the OLS estimator will tend asymptotically toward the true
parameter vector for both populations. However, for a misspecified model, the
OLS estimator will not tend toward the same value for the two populations
because the estimator compensates for errors in the model and the errors depend
on the input variables. Extrapolation to different populations is always error-
prone, but for misspecified models it is even more so.
Example 7 illustrates this. The example shows that the OLS parameter esti-
mate for a misspecified model is different for two different target populations,
even though the true response is the same for both populations.

EXAMPLE 7 The Effect of the Target Population on the Estimated Parameters

in a Misspecified Model
The true response is assumed to be
yi ¼ 0:09

exp ð0:4x
 Þ + εi (21)
ei  N 0, 0:12
We consider the misspecified model
hðXi ; θÞ ¼ θ1 exp ð0:3Xi Þ
We consider two different target populations, namely x U(5, 15) and x  U(10, 20),
where the notation x U(a, b) means that x has a uniform distribution with lower
limit a and upper limit b. In each case we draw a sample of size 30,000 from the
target population, and estimate θ1 using OLS. The estimated parameter value in each
case is shown in Table 5.4.

TABLE 5.4 OLS Parameters for Two Different Target Populations

Target Population ^
θOLS,1
xU(5, 15) 0.3464

xU(10, 20) 0.5706

The results show that even though the true response is the same for all values of x,
the estimated parameters are different depending on which range of x is sampled.
This is a result of misspecification.
176 SECTION A Background

4.1.3 What to Do When Assumption 1 Is Violated?

If there is evidence that assumption 1 is violated, one should change to a better
fitting model. For fairly simple models with few input variables, the model can
be specifically chosen to satisfy assumption 1. For such data-driven models, we
can in general avoid major violations of assumption 1.

4.2 Assumption 2: Constant Variance of Model Error

Assumption 2 is that model variance is the same for all values of the input vari-
ables. If this is not the case, then assumption 2 is violated. One common type of
violation is where the variance of model error increases as the size of the
response increases. For example, consider biomass over the growing season.
Often model error will be small early in the season when biomass is small,
and increase with larger biomass values.

4.2.1 Checking Assumption 2

Violation of assumption

2 can often be detected from residual plots where resid-
uals ei ¼ yi  f Xi ; ^
θOLS are plotted against the predicted response. The spread
of the residual errors should be comparable for all values of the predicted
response.
Example 8 illustrates a case where this assumption seems to be violated. The
data are wheat grain weights from an experiment where grain weights were
measured at numerous dates in a range of environments. The model is a logistic
model. It is quite clear that the variability of the residuals increases as the pre-
dicted value increases. That is, the variance of residual error does not seem to be
constant, so assumption 2 is not satisfied.

EXAMPLE 8 A Case With Nonconstant Variance of Model Error

Grain weight was measured at numerous dates in various fields at Clermont-Ferrand
in France (Robert et al., 1999). The results concern 16 different wheat cultivars, a
range of levels of nitrogen fertilization, and two different years.
A logistic function was fit to the data using OLS. Figure 5.7 shows standardized
residuals for grain weight
as a function
 of degree days
after anthesis
 (DD) (A stan-
dardized residual is ^εsi ¼ ^εi  ^ε =σ^ε where ^εi ¼ yi  f xi ; ^ θOLS , ^ε is the mean of
the ^εi and σ ε is the standard deviation of the ^εi ). The solid line is a smoothed approx-
imation to average residual error. The average of the standardized residuals is close
to 0 throughout, so assumption 1 can be assumed to hold. However, the spread of the
residuals increases with increasing predicted value, indicating that assumption 2 is
violated.
 Regression Analysis, Frequentist Chapter 5 177

A Logistic L3
Standardized residual

–2

10 20
30 40 50
Predicted value
FIGURE 5.7 Plot of standardized residuals for data and model of Example 8. (From
Robert, N., Huet, S., Hennequet, C., Bouvier, A., 1999. Methodology for choosing a model
for wheat kernel growth. Agronomie, 19, 405–417.)

4.2.2 Consequences if Assumption 2 Is Not Satisfied

What if one ignores unequal variances and estimates the parameters using OLS
anyway? If assumption 1 is satisfied but assumption 2 is not, then the OLS esti-
mator is still asymptotically unbiased and consistent. That is, the expectation of
the parameter estimators will still be equal to the true parameter values, and the
estimators will tend to be closer and closer to the true values as the sample size
increases. However, for a given sample size, there will be more variability in the
estimators around the true values than if we transformed to constant variance of
model error, as described below. If there is relatively little data so that the
parameter estimators have at best large variance, then having even larger var-
iance can be a major disadvantage.

4.2.3 What to Do When Assumption 2 Is Violated?

[Link] Weighted Least Squares
Suppose that we know or can estimate the variance of model error at each xi,
noted σ i2. For example, this might be based on measurements on multiple
individuals at the same xi. If assumption 1 is satisfied, then the variability
between individuals with the same x is a measure of variance of model error
for that xi. If we divide each yi and each f(xi; θ) by σ i, then model error after this
transformation is
ε∗i ¼ yi =σ i  f ðxi ; θÞ=σ i ¼ εi =σ i
The variance of this new model error is



var ε∗i ¼ 1=σ 2i varðεi Þ ¼ 1
178 SECTION A Background

Thus, after the transformation, the residual variance is identical (and equal to 1)
for all xi. Now assumption 2 is satisfied and we can use OLS.
The sum of squared errors to be minimized is now
X
n X
n

SSE ¼ ½yi =σ i  f ðxi ; θÞ=σ i 2 ¼ 1=σ 2i ½yi  f ðxi ; θÞ2 (22)
i¼1 i¼1

The approach that consists of minimizing the above SSE is called weighted least
squares (WLS) because it is equivalent to weighting each term in OLS by 1/σ i2.
The R function nls allows one to specify a weight for each measurement, and
thus can be used for WLS.
According to Equation (22), we should give less weight to xi values with
large residual variance. This is logical. We should attach less importance to
xi values where there is more variability around f(xi; θ) because then the mea-
sured value contains less information about the true parameter values. If we do
not use different weights, then we give as much importance to points with little
information as to points with high information.

[Link] Power Transformation

Suppose that the residuals after parameter estimation by OLS indicate that
residual variance is approximately proportional to some power of the expected
value. In this case we can use a power transformation of both measured and sim-
ulated values, as proposed by (Box and Cox, 1964). The measured and simu-
lated values are transformed as

λ 
∗ yi  1 =λ λ 6¼ 0
yi ¼
log ðyi Þ λ¼0
( 
λ
f ð x ; θ Þ  1 =λ λ 6¼ 0
f ∗ ðx i ; θ Þ ¼ i
log ðf ðxi ; θÞÞ λ¼0
where lambda is an additional parameter whose value needs to be chosen by the
user or estimated from the data. The quantity to be minimized by OLS then
becomes
X
N  2
SSE ¼ y∗i  f ∗ ðxi ; θÞ (23)
i¼1

The variance of the transformed variable is approximately (Seber and Wild,

1989, p. 69)

varðyi ∗ Þ  varðyi Þ½Eðyi j xi Þ2λ2 (24)

If var(yi)  [E(yi j xi)]κ and we choose λ ¼ 1  κ/2, then varðyi ∗ Þ  1.
Two cases are of particular interest. If the variance of yi is proportional to
E(yi j xi)2, then according to Equation (24) choosing the transformation with
 Regression Analysis, Frequentist Chapter 5 179

λ ¼ 0 (the log transformation) will make varðyi ∗ Þ approximately constant. If the

variance of yi is proportional to E(yi j xi), then choosing the transformation with
λ ¼ 0.5 will make varðyi ∗ Þ approximately constant.
In the general case, one needs to estimate the additional parameter λ. One
possibility is to try various values of λ and choose the transformation that gives
residual variances that seem more or less constant for all Xi. Alternatively, in the
maximum likelihood framework, one can estimate simultaneously the
parameter vector θ and the variance parameter λ using maximum likelihood.
This involves solving for the vector (θ, σ 2, λ) that maximizes the likelihood.
This is illustrated in Example 9

EXAMPLE 9 Estimation of a Variance Parameter Using ML

Using the data presented in Example 8, Robert et al. (1999) tested the logistic equa-
tion with two alternative assumptions concerning the variance of model error. The
variance was either considered to be constant or to be of the form var(εi) ¼ σ 2f(xi; θ)λ.
The parameters of the logistic equation, and also σ 2 and λ, were estimated using ML.
The estimated value of λ was 1.29.
Figure 5.8 shows the standardized residuals for the raw data (left) and after the
power transformation (right). The residuals on the right have a more uniform distri-
bution, so it seems that the power transformation has allowed us to more closely sat-
isfy the assumption of constant variance.

Logistic L3 Logistic L3V

Standardized residual

Standardized residual

2
2
1
0
0
–1
–2 –2
–3
10 20 30 40 50 10 20 30 40 50
Predicted value Predicted value

FIGURE 5.8 Residual errors without power transformation (left) and with power transforma-
tion with λ ¼ 1.29 (right). (From Robert, N., Huet, S., Hennequet, C., Bouvier, A., 1999. Method-
ology for choosing a model for wheat kernel growth. Agronomie, 19, 405–417.)

[Link] Concentrated Likelihood

Another method of handling unequal variances is the method of concentrated
likelihood, a maximum likelihood method (Seber and Wild, 1989). This method
is adapted to the case where the model calculates multiple responses, each with
a different variance of model error.
Suppose that all model errors are independent and normally distributed, and
that the variances are all equal for a single response but different for different
responses (e.g., for yield and LAI). The likelihood is then
180 SECTION A Background

Y 1 ½yij fj ðxi ; θÞ =ð2σ2j Þ

2
LðYj θ, ΣÞ ¼  1=2 e (25)
i, j 2πσ 2j

where the product is over the different responses (subscript j) and over measure-
ments of that response (subscript i).
The concentrated likelihood method proceeds in two stages. First, treating θ
as fixed, one calculates the variances of model error that maximize the likeli-
hood. This gives

X
nj
 2
σ^2j ðθÞ ¼ 1=nj Yij  fj ðXi ; θÞ (26)
i¼1

where nj is the number of measurements of response j. Then Equation (26) is

substituted into Equation (25) to give a “concentrated” likelihood that involves
only θ. Maximizing that likelihood is equivalent to minimizing the expression
below:
( )nj =2
Y
X nj
 2
^
θ ¼ arg min 1=nj Yij  fj ðXi ; θÞ (27)
θ j i¼1

Equation (27) gives the ML estimator of θ. The minimization can be done, for
example, using the R function optim. The estimated variances of model error
can be calculated from Equation (26), by putting ^θ in place of θ.
The criterion to be minimized in Equation (27) is a product of squared
differences and not a sum, like OLS. Note that this criterion is invariant
to changes in units. Changing the units of any of the outputs would simply
multiply the criterion by a constant, which would not change the estimated
parameter values. This is what we want. We do not want the estimated
parameters to change because we express biomass in g/m2 instead of tons/
ha for example.
An example of the use of concentrated likelihood is the study by Wallach
et al. (2011). The data consisted of measurements of both LAI and aboveground
biomass in multiple fields. It was assumed that model variance was the same for
all observations of each response, but different for the two responses.

4.3 Assumption 3. All Residual Errors Are Independent

Assumption 3 is that all model errors are independent. That is, knowing residual
error for one measurement gives no information about model errors for other
measurements. Violations of this assumption often result from the sampling
procedure, in particular when the same individual is sampled several times.
For example, the data might be biomass measured several times during the
 Regression Analysis, Frequentist Chapter 5 181

course of the season in each of several fields. The fields might be chosen at ran-
dom, but once the field is chosen, then multiple measurements are done in the
same field. The particularities of that field could then influence all the measure-
ments, which would lead to some relation between model errors within the
same field.
The variance-covariance matrix of model errors, noted Σ, is central to con-
sideration of nonindependence of model errors. This is a square matrix of size n
by n, where n is the number of measurements. This matrix has the variances of
the model errors on the diagonal and the covariances in the off-diagonal posi-
tions. If all model errors have the same variance and are uncorrelated (assump-
tions 2 and 3 are satisfied), then Σ has the form
0 1
σ2 0 … 0 0
B0 σ2 … 0 0C
B C
Σ¼B
B⋮ ⋮ ⋱ ⋮ ⋮C C (28)
@0 0 … σ2 0A
0 0 … 0 σ2

There are no off-diagonal terms because all model errors are independent and
therefore all the covariances are zero. All the elements on the diagonal are iden-
tical because we are assuming equal variances. If assumption 3 is violated, then
some of the off-diagonal elements of Σ are nonzero.

4.3.1 Checking Assumption 3

Nonindependent model error shows up on residual plots in that repeated mea-
surements on the same individual will have related residuals. Example 10 illus-
trates this. Seed weights measured in three different years are presented. All the
residuals except the first for year 1986 are positive, and all those for 1988 are
negative. This type of similarity between different residuals for the same field is
a typical symptom of nonindependence of residual error.

EXAMPLE 10 A Case With Nonindependent Model Errors

The data here are grain weights of the spring wheat cultivar Neepawa in three dif-
ferent years (Darroch and Baker, 1990). The parameters of the logistic model
W
seed:weight ¼
1 + e BC ∗DD
were estimated by OLS using all the data. The left panel of Figure 5.9 shows the fit of
the model to the data while the right panel shows the residuals.
It can be seen that residuals for the same year tend to have similarities; for
example, all are negative for 1988. This is an indication that assumption 3 is
violated.

Continued
182 SECTION A Background

EXAMPLE 10 A Case With Nonindependent Model Errors—cont’d

6
1986
4 1987
Seed weight (mg)

30 1988

Residual
2
20
0
1986 −2
10 1987
1988 −4
5
100 200 300 400 500 600 700 100 200 300 400 500 600 700
DD after anthesis DD after anthesis

FIGURE 5.9 Seed weight data for three years and logistic model fitted using OLS (left),
residuals (right).

4.3.2 Consequences if Assumption 3 Is Not Satisfied

What if one ignores correlations and estimates the parameters using OLS? If
assumption 1 is satisfied but assumption 3 is not, then the OLS estimator is still
asymptotically unbiased and consistent, as for violations of assumption 2. That
is, asymptotically the expectation of the parameter estimators will still be equal
to the true parameter values, and the estimators will tend to be closer and closer
to the true values as the sample size increases. However, there will be more var-
iability (between possible samples) in the parameter estimators than if we had
independent model errors. In addition, we are then ignoring the fact that the
actual amount of information in the data is less than with independent model
errors. That is, the effective sample size is smaller than the number of data
points (see Chapter 2). As a result, estimates of the standard error of the param-
eter estimators that ignore nonindependence will be unrealistically small.

4.3.3 What to Do if Assumption 3 Is Violated?

[Link] Generalized Least Squares (GLS)
Suppose that one knows or can estimate the variance-covariance matrix Σ. If we
have multiple measurements at each value of X, those values could be used to
estimate the variances and covariances. For example, if in all three fields in
Example 10 the measurements were made at the same DD values, then we could
estimate the variance-covariance matrix using those values. Another possibility
is to begin with OLS, ignoring nonindependence of model errors. The model
residuals are a first estimate of model errors. The variance-covariance matrix
of those model residuals is then taken as an estimate of the variance-covariance
matrix of model errors.
Given an estimate of Σ one can transform the measured and simulated
values so that the transformed variables are independent. The first step is to
express the variance-covariance matrix as Σ ¼ RTR where R is an upper trian-
gular matrix (This is known as the Cholesky decomposition, and is done by
 Regression Analysis, Frequentist Chapter 5 183

the R function chol). Let Y be the vector of all measurements, so Y ¼ (y1, …, yn),
let F be the corresponding vector of simulated values, so F ¼ (f(x1; θ), …,
f(xn; θ)), and Ε the vector of residuals, so Ε ¼ (ε1, …, εn). At the second step,
the measured and simulated values are transformed as follows:

1
1
Y ∗ ¼ RT Y and F∗ ðθÞ ¼ RT FðθÞ (29)
The transformed residuals are Ε∗ ¼ Y∗  F∗(θ) ¼ (RT)1Ε. The variance-
covariance matrix of the transformed residuals is

 

1
1 1
varðΕ∗ Þ ¼ var RT Ε ¼ RT ΣR1 ¼ RT RT RR1 ¼ I

where I is the identity matrix (the diagonal elements are all 1, the off-diagonal
elements all 0). Thus after transformation, the variance covariance matrix of the
residuals has the form of Equation (28) (with σ 2 ¼ 1). This means that we can
now use OLS on the transformed data.
The SSE criterion to be minimized is

SSE ¼ ½Y ∗  f ∗ ðx; θÞT ½Y ∗  f ∗ ðx; θÞ

1
¼ ½Y  f ðx; θÞT R1 RT ½Y  f ðx; θÞ (30)
T 1
¼ ½Y  f ðx; θÞ Σ ½Y  f ðx; θÞ

The first line of Equation (30) says that we can use OLS on the transformed
measurements and simulated values. The last line of Equation (30) says that
equivalently, we can work with the original measurements and simulated
values, but then the SSE to be minimized includes the inverse of the
variance-covariance matrix. Minimizing this last expression is called general-
ized least squares.

[Link] Determinant Criterion

Suppose that the blocks in the variance-covariance matrix are all identical. For
example, suppose that we have measurements in multiple fields, and all mea-
surements are done at the same times. This could, for example, be measure-
ments of LAI at certain times, biomass at other times, and yield at harvest.
Then all fields would have the same variance-covariance matrix. We also
assume that model errors have a normal distribution. (This is a maximum like-
lihood technique.) In this case we can estimate parameters using the determi-
nant criterion (Bates and Watts, 1988).
We write the measured values and the model values as a matrix. Each
column
0 corresponds1to a 0 single individual (e.g., a 1single field). Thus
y1,1 … yN, 1 f ðx1,1 ; θÞ … f ðxN,1 ; θÞ
Y ¼@ ⋮ ⋱ ⋮ A F¼@ ⋮ ⋱ ⋮ A Z ¼ Y  F where
y1, M … yN, M f ðx1,M ; θÞ … f ðxN, M ; θÞ
184 SECTION A Background

M is the number of measurements per individual and N is the number of

individuals. Z is the matrix of residuals.
In this case the likelihood can be shown to be
1
etr½ðZ Σ1 Z =2
T
LðYj θ, ΣÞ∝ N=2
(31)
jΣj
where Σ is the variance covariance matrix for each individual (the same for all
individuals), the vertical bars indicate a determinant, and the trace (tr) is the sum
of the diagonal elements of a matrix. This likelihood can be maximized in two
steps. First, one solves for the value of Σ that maximizes the likelihood for
given θ. The result is
T
^ ðθ Þ ¼ Z Z
Σ
N
Then one substitutes that into Equation (31). The resulting expression is

  
L Y; θ, Σ^ ðθÞ ¼ k0  ðN=2Þ ln Z T Z (32)
where k0 is a constant. Thus, the maximum likelihood estimator of θ is the value
of θ that minimizes the determinant jZTZ j.
This result is remarkably general. It applies whatever the structure of the
variance-covariance matrix, as long as it is the same for all individuals. The other
requirement is that the number of individuals must be greater than M (the number
of measurements per individual) and greater than p (the number of parameters).

[Link] Random Parameter Models

Here we assume that the model parameters are random variables. They have
fixed values for a given individual but can vary between individuals. Given
the specific parameters for an individual, the model errors for that individual
are assumed independent. Covariances between model errors for different mea-
surements on the same individual arise because measurements on the same indi-
vidual share the same parameter values. One can use the R function nlme to
estimate the parameters for random parameter models.
As a simple example, we return to the data in Example 10. We assume that in
the logistic model, the parameter W is a random variable and that the values of
W for different individuals are independent. Two measurements on the same
individual are then correlated because they share the same value of W:
!

 Wi Wi
cov yij , yij0 ¼ cov + εij ,   + εij0
1 + eBC∗DDij 1 + eBC∗DDij0
varðW Þ
¼  
½1 + e ∗ ij  1 + eBC∗DDij0
BC DD

where the subscript ij refers to the jth measurement on individual i. The random
variables here are W, εij, and εij0 . The model errors have been assumed mutually
 Regression Analysis, Frequentist Chapter 5 185

independent and independent of W, so all the covariances involving the model

errors are zero. However, there remains the term cov(W,W) ¼ var(W), so
cov(yij, yij0 ) is not zero. Thus, the fact of having a random parameter introduces
covariances between different measurements for the same individual.
Random parameter models are often an appealing way to treat nonindepen-
dence of model errors. First, the assumption that differences between individ-
uals arise from differences in their parameter values can often seem quite
reasonable. In addition, random parameter models can be a parsimonious
way of describing nonindependence, if the model has relatively few parameters.
In Example 10, even if we treat all three parameters as random, there are
only six variances and covariances to estimate, namely var(W), var(B), var(C),
cov(W,B), cov(W,C), and cov(B,C).
An example of a random parameter model is presented in Example 11. The
data are wheat yield measurements for several N application rates in each of
multiple fields. Model errors may be correlated for different N rates in the same
field, due to specificities of that field not accounted for by the model. The ran-
dom parameter model assumes that each field has its own parameters, drawn
from a normal distribution. The problem then is to estimate the mean of each
parameter and the variance-covariance matrix of the parameters.

EXAMPLE 11 Nitrogen Response Model With Random Parameters

In Makowski and Wallach (2002), the following equations were used to model (a) N
uptake (Nup) as a function of applied N (Napp) and end of winter soil mineral N
(Nwinter), (b) yield (Yld) as a function of N uptake and applied N, (c) grain protein
content (Prot) as a function of N uptake and yield, and, (d) residual soil nitrogen
at harvest (Res) as a function of applied N.
8
>
> CN app + a1 + a2 Nwinter Napp < Nmax
<


Nup ¼ C Napp  Nmax
>
>
 Napp  Nmax
: max
BY +
1 + T Napp  Nmax
8

< Ymax + A Nup  BY max Napp < Nmax
Yld ¼
:Y
max Napp  Nmax

Nup
Prot ¼ P1 + P2
Yld
(
Rmin Napp < Nmax
Res ¼

Rmin + R Napp  Nmax Napp  Nmax

The model has 11 parameters. All model errors were assumed independent.
Parameter estimation was done in two ways. In the first approach, all parameters were
treated as fixed. Estimation then concerned the 11 fixed parameters and four variances
of model error (for N uptake, yield, grain protein content, and residual soil N). In the
second approach the nine parameters a1, C, Ymax, A, B, P1, P2, Rmin, and R were treated

Continued
186 SECTION A Background

EXAMPLE 11 Nitrogen Response Model With Random Parameters—cont’d

as random parameters. Estimation then concerned the expectations and variances of
these nine random parameters, seven covariances between them (the remaining covari-
ances were assumed to be 0), the values of the two fixed parameters, and the four
variances of model error. Parameter estimation for the random parameter model used
the S-PLUS function nlme, which is analogous to the R function of the same name.
The parameter vector obtained in each approach was used to calculate optimal
fertilizer rates for each field in the dataset. Observed yields were used to estimate the
gross profit (yield * price  amount of N * cost) that would have been earned for each
field, if the calculated optimal N dose had been applied. It was found that the ran-
dom parameter model led to substantially higher gross margins than the model with
fixed parameter values shows that correcting for violations of assumption 3 can have
important practical benefits.

A ALGORITHMS FOR OLS

The advantage of having software with built-in functions for parameter estima-
tion, like R, is that one does not need to program an algorithm that will actually
do the parameter estimation. As long as the algorithms work well, there is no
need to know how the algorithms operate. But when something goes wrong,
it is useful to have some basic understanding of how the algorithms function.
We consider here two commonly used algorithms: the Gauss Newton algo-
rithm (which is the default algorithm for the nls function in R) and the Nelder-
Mead or simplex algorithm (which is the default algorithm for the optim
function in R). These two algorithms operate very differently. The Gauss-
Newton algorithm is specifically built to find the parameter values that mini-
mize a sum of squared errors. The Nelder-Mead algorithm on the other hand
can be used to find the parameters that minimize any function.

A.1 Gauss-Newton Algorithm

We first recall the OLS solution for a linear model (i.e., a model where the
parameters enter linearly). A linear model can be written Y ¼ Xθ in matrix nota-
tion, where Y is the vector of n measured values, θ is the vector of p parameters,
and X is the matrix of input variables. In terms of individual elements
0 1 0 10 1
y1 x1,1 … x1, p θ1
@ ⋮ A ¼ @ ⋮ ⋱ ⋮ A@ ⋮ A (33)
yn xn,1 … xn, p θp
where xi, j is the value of input variable j for individual i. For the linear model
there is an analytical matrix expression for the least squares parameters:

1
^
θOLS ¼ XT X XT Y (34)
 Regression Analysis, Frequentist Chapter 5 187

The basic Gauss Newton algorithm approximates a model that is nonlinear in

the parameters by a linear model and then uses Equation (34) to find new param-
eter values. Because the model is not really linear, the resulting parameters are
not the exact solution. Therefore, one repeats the procedure starting at the new
parameter values. The procedure is iterated until the changes in the parameter
values become very small. In more detail, the algorithm works as follows:
ðiÞ
a. Define an initial estimation for θ, noted ^
θ with i ¼ 0.
b. Use a Taylor series expansion to derive a linear approximation to the non-
linear model. The approximation is:
 ðiÞ  X 
∂f ðx; θÞ  
p
ðiÞ
f ðx; θÞ  f x; ^
θ +  θj  ^θj
j¼1
∂θj ^θðiÞ

∂f ðx; θÞ
where is the derivative of f(x; θ) with respect to the jth parameter,
∂θj ^θðiÞ
ðiÞ
evaluated at ^θ . We now have an equation like Equation (33), with the
replacements
0  1
ðiÞ 0 1
y1  f x1 ; ^
θ y1
B C
B ⋮ C ! @⋮A
@  ðiÞ
 A replaces
yn  f xn ; ^
θ yn
0   1
∂f ðx1 ; θÞ ∂f ðx1 ; θÞ 0 1
B ∂θ  ðiÞ … ∂θp ^θðiÞ C x1, 1 … x1, p
B 1 ^θ C
B ⋮  C @ ⋮ ⋱ ⋮ A
B ⋮  ⋱ C !
@ ∂f ðx1 ; θÞ ∂f ðx1 ; θÞ A replaces xn, 1 … xn, p
…
∂θ1 ^θðiÞ ∂θp ^θðiÞ
0 ð iÞ
1 0 1
θ1  ^
θ1 θ1
B C
@ ⋮ A ! @⋮A
ð iÞ replaces
θp  ^
θp
θp

ðiÞ
c. Calculate the ordinary least squares estimate of θ  ^θ as in Equation (34),
but with the above replacements. This gives the next approximation,
ði + 1Þ ðiÞ
 ð iÞ

^
θ ¼^
θ + θ  ^θ .
d. Repeat steps b–c until convergence is achieved (i.e., until there is very little
change in the parameters or in the sum of squared errors).
The behavior of this algorithm is illustrated in Example 12, using the logistic
model and data of Example 1. Note that the first iteration leads to values that
are much worse (higher MSE) than the initial values. Both parameters after the
first iteration are negative. In this particular case, the model can still be evalu-
ated with negative parameter values, and the algorithm goes on to find the OLS
188 SECTION A Background

solution. With other models, negative parameter values might cause an error. To
avoid this problem, it is important to start with good initial estimates of the
parameter values. Often, one can roughly estimate some or all of the parameters
directly from the data. It is worthwhile doing so. If the algorithm does not con-
verge, one can try different starting values.

EXAMPLE 12 Path of Algorithm Used by nls

Consider the data and model of Example 1. In order to work in only two dimensions
(much easier to visualize), here we fix the parameter W ¼ 28.171259, which is its
OLS value, and estimate the two remaining parameters of the logistic curve,
B and C, using the Gauss-Newton algorithm of the R function nls.
Figure 5.10 shows a contour plot of mean squared error, and the trajectory of the
nls iterations. The starting point is indicated. The first iteration leads to a value
(B ¼  1.698, C ¼  0.014) that is, in fact, much worse (higher SSE) than the starting
value. The algorithm then reduces the step size, then reduces it again to find a smal-
ler SSE value than the initial value. The algorithm then changes direction and goes
quite directly to the OLS solution ( ). Overall, seven iterations were required.

nls
0.03

0.02
start
67
0.01 5
4
10
C

20
500
30
0.00 50
40
100
3

300

−0.01

2
−0.02

−2 −1 0 1 2 3 4 5
B
FIGURE 5.10 Path of Gauss-Newton algorithm for estimation of two parameters of
logistic model.

A.2 Nelder-Mead Algorithm

The Nelder-Mead algorithm (also called the simplex algorithm) does not use
derivatives, and is considered more robust than the Gauss-Newton algorithm.
It may converge when the Gauss-Newton algorithm does not. This is the default
algorithm of the R function optim.
 Regression Analysis, Frequentist Chapter 5 189

As with the Gauss-Newton algorithm, the user must give initial estimates of
the parameter values. This algorithm then begins by creating a simplex, which is
a geometric figure with p + 1 vertices in p-dimensional space where p is the
number of parameters. For two parameters this is just a triangle. For three
parameters it is a tetrahedron.
The algorithm evaluates the SSE at each vertex of the simplex. Then the ver-
tex with the highest SSE (the worst point) is reflected through the centroid of the
other points (reflection). If this new point is better than the current best point,
the algorithm elongates the distance of reflection to go further (expansion). If
the new point is no better than the second worst point, then the distance of
reflection is reduced (contraction). If no new point better than the second worst
point is found, the original simplex is shrunk toward the best point (shrinkage).

EXAMPLE 13 Steps in the Execution of the Nelder-Mead Algorithm

The first steps of the Nelder-Mead algorithm applied to the same model and data as
Example 12 are shown in Figure 5.11. Once again we estimate two parameters, B
and C, of the logistic function Equation (7). The initial point given by the user is indi-
cated by “start.” First (upper left panel), the simplex (here a triangle) is created by the
algorithm. Then (upper right panel) the worst point (highest SSE) is reflected through
the center of the line connecting the other two points (reflection). This new point is
worse (higher SSE) than the other two points, so it is not accepted. The distance of
reflection is reduced until an SSE value lower than the second worst point is found
(lower right panel). The new point is used to create a new triangle, where the pre-
viously second worst point is now the worst (bottom right panel).

Nelder-Mead, initial triangle Reflection

0.4 0.4
worst

0.0 0.0
C

start

−0.4 −0.4
2.8 3.0 3.2 3.4 2.8 3.0 3.2 3.4
B B

Reflection New triangle

0.4 0.4

0.0 0.0
C

worst

−0.4 −0.4
2.8 3.0 3.2 3.4 2.8 3.0 3.2 3.4
B B
FIGURE 5.11 Steps in the Nelder-Mead algorithm.
190 SECTION A Background

A.3 Comparison of Algorithms

Example 14 shows a comparison of the Gauss-Newton and Nelder-Mead algo-
rithms as implemented in R. Here, all three parameters of the logistic model
were estimated, using the data of Example 1. The final parameter values calcu-
lated by the two algorithms are very similar. Nelder-Mead requires many more
model evaluations than Gauss-Newton, which is fairly typical behavior.

EXAMPLE 14 Comparison of the Gauss-Newton and Nelder-Mead Algorithms

as Implemented in R
The three parameters of the logistic model were estimated using the data of
Example 1. The estimation was done using the R function nls with its default algo-
rithm (Gauss-Newton), and with the R function optim using its default algorithm
(Nelder-Mead). The nls function announced “Number of iterations to convergence:
8,” but this does not count the function evaluations needed to estimate the deriva-
tives. Each iteration requires one evaluation to evaluate SSE plus three evaluations to
estimate the derivatives with respect to the three parameters W, B, and C.
(In addition, there are a few extra function evaluations.) Both algorithms converged
to essentially the same point, but Nelder-Mead required almost eight times as many
evaluations of the model.

Number of
Evaluations
Method W B C of the Model
Gauss-Newton 28.1713 4.16022 0.0164920 37
(nls)
Nelder-Mean 28.1712 4.15961 0.0164895 285
(optim)

There are many variants of the Gauss-Newton and simplex algorithms, and
also many other types of algorithms, some of which are available with nls or
optim as options. For example, an option with optim is to use a simulated
annealing algorithm, which is a global optimization algorithm (Goffe et al.,
1994). Global algorithms are less sensitive to initial values than Gauss-Newton
type algorithms and are more likely to converge to the optimal solution. An
important drawback of these algorithms is that they often require a long calcu-
lation time. An application of simulated annealing to a crop model is described
by Mavromatis et al. (2001).

A.4 Potential Problems

There are two different types of problems that arise with these algorithms. The
first is that the algorithm may fail to converge. A commonly encountered error
message for nls is of the form
 Regression Analysis, Frequentist Chapter 5 191

Error in nls(y  W/(1 + exp(B  C * x)), obs, list(W ¼ 15, B ¼ 1, C ¼ 0.001)):

singular gradient

The “singular gradient” message means that one of the derivatives is zero, and
so the algorithm cannot continue. This may be because the parameters are uni-
dentifiable or because the starting values are poor and lead the algorithm to
evaluate the model at parameter values that create errors. Nelder-Mead may
also fail to converge if the model is evaluated at a point that causes an error.
The second problem is that the algorithm may seem to converge, but the result
is not the true minimum value of the function but rather a local minimum. This is
more pernicious because one might not even suspect that there is a problem.
A protection here is to start the algorithm at several different initial parameter
values.
It is important to keep the ratio of the number of estimated parameters to the
number of observations relatively low. Many numeric problems result from
trying to estimate numerous parameters simultaneously from a small number
of data. This is in addition to the danger of obtaining very poor parameter
estimates.

B R FUNCTIONS FOR PARAMETER ESTIMATION

We discuss here two R functions that can be used for parameter estimation. The
first is nls (for nonlinear least squares), which is specifically designed to find
parameter values that minimize a sum of squared errors between observed
values and values simulated using a nonlinear model. The second is optim,
which is an R function that can be used to find the parameter values that min-
imize any function.

B.1 The nls Function When the Model Consists of a Single

Equation
The simplest general form of the nls function is the following:
nls(formula, data, start)

where formula is a formula expression that relates the measurements to the

model, data gives the data (measured values and explanatory variables of
the model), and start gives the initial guesses for the parameter values. The
following example calculates the OLS parameters for the grain weight data
and logistic model of Example 1.
obs<-c(8,13,17,21,26,27,28,27,27,29,28)
DD <-c(195,240,280,330,380,410,480,530,580,620,680)
data<-[Link](y =obs,x=DD)
start <-list(W =29,B=4.2,C= 0.017)
logistic<-nls(yW/(1+exp(B-C*x)),data,start)

The first three lines store the measured grain weight values in the vector obs, the
values of the explanatory variable DD in the vector DD, and then create a data
192 SECTION A Background

frame with two columns; the values of obs are in a column named y and the
values DD in a column named DD. The names of the columns are important.
They must be the same names that appear in formula. The fourth line creates
a list called start, with initial guesses of the value of each parameter. Here
again the names must be the same as the names that appear in formula.
The last line is the call to the nls function. The first argument, the formula,
says that “y is modeled as W/(1 + exp(B  C*x).” The tilde indicates that this is a
formula. The R objects that appear in this expression are y, W, B, C, and x. Every
one of these objects must appear by name and be assigned a value in either the
second or third argument. If the name appears in the third argument (the object
start in our example) it is treated as a parameter to be estimated. In our case the
parameters are W, B, and C. The remaining names in formula are y and DD; they
are defined in the second argument (the object data in our example).
In fact, you do not really have to define every variable that is not a parameter
in data. If objects used in the formula are not found in data, then R will search
for these objects among the objects that have been defined previously. How-
ever, it is good programming practice to include all the variables in data to
insure that R is using the values you want.
In the above program, the result of executing nls is stored in the object
logistic. The basic results of nls can be printed by simply putting logistic
on a line by itself. The result is as follows, where the values given for W, B, and
C are their estimated values:
logistic #print the results of nls
Nonlinear regression model
model: y  W/(1 + exp(B – C * x))
data: data
WBC
28.17126 4.16022 0.01649
residual sum-of-squares: 5.596
Number of iterations to convergence: 5
Achieved convergence tolerance: 3.181e-06

One can obtain more detailed information by using the instruction summary
(logistic):

summary(logistic,correlation=T)
Formula: y  W/(1 + exp(B – C * x))
Parameters:
Estimate Std. Error t value Pr(>jtj)
W 28.171259 0.408285 69.00 2.17e-12 ***
B 4.160219 0.360084 11.55 2.86e-06 ***
C 0.016492 0.001415 11.65 2.68e-06 ***
Signif. codes : 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.8364 on 8 degrees of freedom
 Regression Analysis, Frequentist Chapter 5 193

Correlation of Parameter Estimates:

WB
B -0.43
C -0.53 0.98
Number of iterations to convergence: 5
Achieved convergence tolerance: 3.181e-06

Here we have as additional information the standard deviations of the parameter

estimators and the statistical significance of the parameters (are they signifi-
cantly different than 0, usually not a very useful question for nonlinear models).
We have asked for the correlations between parameters (“correlation ¼ T” in
the summary function) and so we also have the correlations. The correlation
between the estimators of B and C is quite high (0.98).

B.1.1 Recovering and Analyzing the nls Results

The result of executing the nls function is an object (an nls object, logically
enough) that contains the results of the estimation. There are specialized
instructions for recovering or analyzing the results of nls.
The function coef operates on nls objects and returns the estimated values
of the parameters.
# get coefficients
coef(logistic)
WBC
28.17125873 4.16021931 0.01649193

The function fitted operates on nls objects and returns the values calculated
using the calibrated model. Below we use this function to obtain a graph of
observed versus simulated values, to which we add the 1:1 line (Figure 5.12).
fitted(logistic)
[1] 7.888914 12.666716 17.252816 22.052749 25.116203 26.226212
27.527684 27.885441 28.045235 28.105961 28.146948
# graph observed versus predicted values and add 1:1 line
plot(DD,fitted(logistic))
abline(0,1)

The function residuals returns the differences between the observed and sim-
ulated values. Below we plot residuals versus simulated values, and add the line
at 0. We also calculate mean squared error (MSE) (Figure 5.13).
# graph residuals
[Link]()
plot(DD,residuals(logistic), xlab="growing DD since anthesis",
ylab="residulas")
abline(0,0)
194 SECTION A Background

Fitted ([Link]) 25

20

15

10

10 15 20 25
y
FIGURE 5.12 Observed seed weight values versus simulated values, using the fitted function.

0.5
Residulas

0.0

−0.5

−1.0

200 300 400 500 600

Growing DD since anthesis
FIGURE 5.13 Residuals of logistic model fitted to seed weight data.
 Regression Analysis, Frequentist Chapter 5 195

MSE<-mean(residuals([Link])^2)
MSE
[1] 0.5087239

The function predict uses the regression formula, with the estimated parame-
ters, to evaluate the model for new values of the explanatory variables. The first
argument is the name of the nls object. The second argument is a data frame or
list with the values of the new explanatory variables. In the following code, we
use the calibrated model to predict for each integral value of DD from 195 to
680. Then we plot the observed values and the predicted values (Figure 5.14).
# plot observed values and predicted response for each integer value
of DD
[Link]<-195:680
[Link]<-predict(logistic,list(x=[Link]))
plot(DD,obs, xlab="growing DD since anthesis",ylab="seed weight")
lines([Link],[Link])

The function AIC operates on an nls object and calculates the value of AIC (a
measure of model quality) for the model with the estimated parameters.
# AIC for logistic model with p= 4 parameters (W, B, C and residual
variance)
AIC(logistic)
[1] 31.7823

25
Seed weight

20

15

10

200 300 400 500 600

Growing DD since anthesis
FIGURE 5.14 Measured seed weight values and curve of values simulated by logistic model.
196 SECTION A Background

B.2 nls for More Complex Models

Often one wants to estimate parameters for more complex models that do not
have the form of a single equation. This is obviously the case for complex sys-
tem models. The structure of the nls instruction is then the same as for a one-
equation model, but in place of the explicit formula one has a call to the function
that evaluates the model. That function must calculate and return a vector of
simulated values that corresponds to the vector of measured values.
To illustrate, we use the same example as above, which only has a single
equation, the logistic equation, but we put it into a function. The program below
has first the function that evaluates the model and then the main program with
the nls instruction.
# function to calculate grain weight using a logistic equation
# x is a vector, so pred is a vector of simulated values
[Link]<-function(x,W,B,C) {
pred<-W/(1+exp(B-C*x)) # the vector of predicted values
return(pred)
}
# main program
obs<-c(8,13,17,21,26,27,28,27,27,29,28)
DD<-c(195,240,280,330,380,410,480,530,580,620,680)
data<-[Link](y=obs,x= DD)
start<- list(W=29,B=4.2,C =0.017)
logistic<-nls(y[Link](x,W,B,C),data,start)

In the first argument of the nls function, the formula says that the model is eval-
uated in the function [Link], which has four arguments. Therefore, the
definition of [Link] must also show four arguments. The arguments are
passed by position (first argument in the call from nls is used as first argument
in [Link], etc.). The names do not need to be identical in the call from
nls and in the function definition, though in this case they are. As in the simple
case, all the object names that appear in the formula (in our case y, x, W, B, and
C) must be assigned values either in the second or third arguments to nls.

B.2.1 nls With a Wrapper Function

Suppose that the model has already been programmed as a function in R. Sup-
pose further that the model function simulates for just a single context (soil,
year, management), and for that context it outputs daily values of all the state
variables. This is typical of crop models. The measured data, on the other hand,
typically come from several different contexts, and for each context there are
only data for selected days and certain output variables. Thus, the model func-
tion is not adapted to nls, which requires that the model function calculate a
vector that is directly comparable to the vector of measurements.
 Regression Analysis, Frequentist Chapter 5 197

One possibility would be to modify the model function. But in general, we

want to keep the model function unchanged because it will also be used for other
purposes than parameter estimation. The solution here is to write a second func-
tion that is called by the nls function and in turns calls the model function. This
new function has three main responsibilities. First, at each iteration (each new
parameter vector sent by nls), it must define the values of the model parameters,
those that are kept fixed and those that are being varied. Second, it must call the
model function several times, once for each context for which there are data.
Third, it must separate from the results of the model function just the simulated
values that correspond to the measured values, and return those values. This
intermediary function is a “wrapper” function in the sense that it is covering
the model function so that the nls instruction does not directly see it.
In the following example, we estimate the parameters for the simple maize
model. Often for system models, we only want to estimate a subset of all the
parameters. Here we will estimate just two of the parameters, namely LAImax
and alpha.
The data are LAI values at 60 days after emergence for 10 different fields in
2005. They are stored in the vector [Link]. The input variables are daily
weather (there are 15 different weather variables). Weather data are stored in
a three-dimensional array (days by climate variables by site) called weather.
The model parameter values are stored in the vector param. These are the fixed
values for those parameters that are not being estimated, and the starting guesses
for the parameters that are being estimated.
[Link]
[1] 6.087119 5.191289 3.618078 3.536665 4.709734 7.004336 3.648694
6.958858
[9] 6.219985 8.472625
weather[1:5,1]
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
[,14] [,15]
[1,] 10 43.4 -0.55 2005 1 6.5 12.6 4.3 0 2.2 5.1 7.3 85.6 145 0.67
[2,] 10 43.4 -0.55 2005 2 3.2 12.3 5.1 0 3.3 5.1 8.0 81.6 145 0.98
[3,] 10 43.4 -0.55 2005 3 7.5 9.2 1.4 0 1.8 1.0 4.0 80.8 145 0.56
[4,] 10 43.4 -0.55 2005 4 8.2 7.4 0.1 0 2.4 -0.5 2.7 79.3 145 0.61[5,]
10 43.4 -0.55 2005 5 7.6 6.6 0.6 0 2.1 -0.4 2.5 80.6 145 0.5
param
Tbase RUE K alpha LAImax TTM TTL
7.00e +00 1.85e+00 7.00e-01 2.43e-03 7.00e+00 1.20e+03 7.00e+02

The wrapper function that calls the model to predict values of LAI60 is shown
below. It first sets the values of the parameters to the new values sent from nls,
for those parameters that are being estimated. Then there is a loop over sites. For
each value of i, the wrapper prepares the weather data for the ith site, then calls
the maize model and extracts from the results the calculated LAI on day 60 after
198 SECTION A Background

sowing. It returns a vector of predicted values in the same order as the observed
values.
# function called by nls to calculate LAI on day 60
[Link]<-function(nsites,param,weather,sdate,ldate,
LAImax,alpha)
{
# replace the default values of LAImax and alpha with the current values
from nls
param["LAImax"]<-LAImax
param["alpha"]<-alpha
for (i in 1:nsites)
{
# extract weather for current site, and give names to climate variables
[Link]<-[Link](weather[,c(1,4:8),i])
names([Link])<-c("site","year","day","I","Tmax","Tmin")
# Call the maize model for this site-year, and extract LAI on day 60
[Link]<-[Link](param,[Link],sdate,ldate)
LAI60[i]<-[Link]$LAI[[Link]$day==(sdate+60)]
}
# return the vector of predicted values
return(LAI60)
}

The following shows the nls instruction that calls the wrapper function, and the
results.
# nls calls the [Link] wrapper function
# lines that start with + are continuations of previous instruction
nls.LAI60 <-nls([Link][Link](nsites,param,weather,
sdate,ldate,LAImax,alpha),
+ list(nsites =nsites,param =param,weather=weather,sdate= sdate,
ldate=ldate),
+ list(LAImax =7,alpha=0.0024))
summary(nls.LAI60)
Formula: [Link]  [Link](nsites, param, weather, sdate, ldate,
LAImax, alpha)
Parameters:
Estimate Std. Error t value Pr(>jtj)
LAImax 7.9410018 0.7813556 10.163 7.52e-06 ***
alpha 0.0021693 0.0002872 7.554 6.58e-05 ***
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.8102 on 8 degrees of freedom
Number of iterations to convergence: 7
Achieved convergence tolerance: 2.034e-06
 Regression Analysis, Frequentist Chapter 5 199

B.3 The optim Function for Function Minimization

There are cases where the nls function is unable to converge to the OLS param-
eter values. This may be a problem of the starting values or because the function
is ill-behaved. In this case it is worthwhile to attempt to use the optim function
with a simplex algorithm, which is generally more robust than the Gauss-
Newton algorithm of nls. optim is a general purpose function minimization
algorithm. It is not specifically designed for OLS problems.
The general form of the optim function is
optim([Link],[Link],[Link],
method ¼”Nelder-Mead”)
The first argument is a vector of initial values of the parameters to be estimated,
the second argument is the name of the function that evaluates the model and
returns the sum of squared errors for the current parameter values, the third
argument is not needed here (hence the two commas), and the next argument
or arguments is passed to [Link]. The method argument specifies
the minimization algorithm to be used.
We illustrate using the grain weight data and the logistic model of
Example 1. Below is the function ([Link]) that calculates the sum
of squared errors in this case. The first argument of the function must be the
vector of values of the parameters that are to be estimated. At each iteration,
optim calls the function with the current values of those parameters. The func-
tion can have any other argument. They must be in the same order as in the
optim instruction. In the example, the additional arguments are x (growing
degree days, needed to evaluate the logistic model) and y (observed grain
weights, needed to calculate the sum of squared errors).
[Link] <-function(params,x,y) {
W <-params[1]
B <-params[2]
C <-params[3]
pred<-W/(1+exp(B-C*x)) # the vector of predicted values
SS <-sum((y-pred)^2) # the sum of squared errors
return(SS)
}

Then the use of optim is as follows:

# estimate parameters using optim
params<-c(29,4.2,0.017) # initial guesses of parameter values
logistic<-optim(params,[Link],DD,obs,method="Nelder-
Mead")
logistic
logistic$par
[1] 28.17115372 4.16033032 0.01649234
200 SECTION A Background

$value
[1] 5.595963
$counts
function gradient
204 NA
$convergence
[1] 0
$message
NULL

The value of $value is the final best value of the quantity that is minimized, in
our case the sum of squared errors. $counts gives the number of times the model
function was called (204), and $convergence= 0 indicates successful conver-
gence. Note that the parameter values are in the object logistic$par.
The function that evaluates the model and that is called by optim cannot be
the same as the function that is called by nls because the quantities returned are
different in the two cases. nls requires that the model function return a vector of
predicted values while optim requires the sum of squared differences. That
means, unfortunately, that one cannot simply try nls first and then just replace
nls by optim if nls does not converge. One must replace the function or the
wrapper used with nls by one adapted to optim.
The optim function is not as convenient for OLS parameter estimation as
nls. It does not produce information about the accuracy of the parameter esti-
mates as nls does. Furthermore, the summary function, the predict function,
and other useful functions that take nls objects as arguments do not accept
optim objects as arguments.
There is thus a real advantage to using nls if possible in preference to optim.
If, however, nls does not converge and one uses optim, one can still obtain the
additional information that nls provides, though with additional effort. One
possibility is to use the parameter values estimated by optim as initial guesses
for nls. Perhaps with these presumably very good initial guesses, nls will con-
verge and one can take advantage of the information it provides. If not, one can
write R code to obtain the additional information. We show below how to cal-
culate several of the quantities of interest.

B.3.1 Additional Information When optim Is Used

Because the optim function does not provide a simple way to obtain all the results
of interest after OLS parameter estimation, these must be programmed by the
user. Here we show how to calculate the estimator of residual variance, the vector
of residuals, and the standard errors and correlations of the parameter estimators.
The example here again uses the grain weight data and logistic model of
Example 1. We suppose that we have run optim and that the estimated values
of the three parameters of the logistic model are in the object logistic$par.
 Regression Analysis, Frequentist Chapter 5 201

For these calculations, we need to have a function that evaluates the model
for all the contexts for which we have data and returns the simulated values. For
our example, that means that the function must evaluate the model for each
value of DD for which we have data.
# function that evaluates the model
[Link]<-function(params,x)
{
W <-params[1]
B <-params[2]
C <-params[3]
pred<-W/(1+exp(B-C*x))
return(pred)
}

We can calculate the vector of residuals and the estimated residual variance and
standard deviation as follows:
# residuals. logistic is the optim object with results of OLS
# obs is observed grain weights. DD is degree days
residualsOptim <-[Link](logistic$par,DD)
# Calculate estimated residual variance and residual standard deviation
sig2<-sum(residualsOptim^2)/(length(DD)-3)
sig2[1]
0.6994954
sig<-sqrt(sig2)
sig
[1] 0.8363584

The estimated standard deviations of the parameter estimators, and their corre-
lations, are given by Equation (16). This involves numerical derivatives for
which we use the R package numDeriv. This package has a function called
jacobian, which calculates the matrix of derivatives we need. This function
has three arguments of interest to us. The first is the name of the function that cal-
culates the model predictions, in our case [Link]. The second is the
vector of parameter values. The derivatives will be calculated with respect to each
parameter at the point represented by the values of the vector. In our case, this sec-
ond argument is logistic$par, which contains the OLS estimates provided by
optim. The third argument (three commas later because we are skipping two argu-
ments) has the other arguments required by the function [Link]. In
our case there is only one, namely DD. If there were other additional arguments,
they would appear in order, separated by commas. The R program is then
library(numDeriv)
# X is a matrix of derivatives. Size n (data points) x p (parameters)
X <-jacobian([Link],logistic$par,DD)
202 SECTION A Background

# Calculate var-covar matrix Sigma

inverseXtX<-solve(t(X)%*%X)
Sigma<-inverseXtX*(sig2)
# Standard deviations of estimators are square roots of diagonal
elements
stdevs<-sqrt(diag(Sigma))
stdevs
[1] 0.408279270 0.360092422 0.001415028
# Correlations
correlations<-((1/stdevs)%*%t(1/stdevs))*Sigma
correlations
[,1] [,2] [,3]
[1,] 1.0000000 -0.4311114 -0.5336516
[2,] -0.4311114 1.0000000 0.9771435
[3,] -0.5336516 0.9771435 1.0000000

The results here are identical to those furnished by the summary function nls.

EXERCISES
Easy
1. In the simple maize model of Chapter 1, the equation for the daily change in
biomass up to flowering is ΔB(d) ¼ B(d + 1)  B(d) ¼ RUE[1  eKLAI(d ))]
I(d). The two explanatory variables are leaf area index on day d,
LAI(d), and solar radiation on day d, I(d). What exactly does it mean to
say that this model for ΔB(d) is correctly specified for some particular target
population?
2. Consider the graph of Figure 5.15 for residuals of some model as a function
of the values simulated by the model.
a. Which of the assumptions underlying OLS (Table 5.1) seems to be
satisfied?
b. Which assumptions do not seem to be satisfied?
c. What action would you suggest? Explain.
3. Consider the graph of Figure 5.16 for residuals of some model as a function
of the values simulated by the model.
a. Which of the assumptions underlying OLS (Table 5.1) seems to be
satisfied?
b. Which assumptions do not seem to be satisfied?
c. What action would you suggest? Explain.
Moderate
4. Consider the graph of Figure 5.16 for residuals of some model for four dif-
ferent individuals with 10 measurements per field.
a. What would be the form of the variance-covariance matrix here?
b. Can you use the determinant criterion here? If not, why not?
 Regression Analysis, Frequentist Chapter 5 203

1.5

1.0

0.5
eps

0.0

−0.5

−1.0

6 8 10 12 14 16 18
yhat
FIGURE 5.15 Graph of residuals of a model versus simulated values.

4 Individual 1
Individual 2
Individual 3
Individual 4

2
eps

−2

−4

6 8 10 12 14 16 18
yhat
FIGURE 5.16 Residuals for multiple measurements for each of four individuals.

5. Suppose that in the sugar cane model described below, the observed yields
are noted yi , i ¼ 1, …, 6 and the corresponding model predictions are noted
f(xi; α, κ, x0).
a. What is the quantity to be minimized for OLS parameter estimation?
204 SECTION A Background

b. Assuming that the assumptions of Table 5.1 are satisfied and that the
residuals have a normal distribution, what quantity does the maximum
likelihood estimator minimize?
6. The following model describes the relation between sugar cane yield y and
applied Nitrogen x using a Mitscherlich equation:
y ¼ αf1  exp ½κ ðx  x0Þg
where α is the asymptotic limit of yield for large x and x0 is applied N for
zero yield. The table below shows data for yield for various N doses (from
Schabenberger and Pierce, 2002, p. 238).
a. Use the nls function to estimate the parameters α, κ, and x0. (Reasonable
starting values are around α ¼ 200, κ ¼ 0.01, x0 ¼ 20.) Use at least two
different starting vectors. What are the OLS parameter estimates?
b. What is the estimated value of residual variance?
c. Create a graph with simulated values every day and observed values as
points.
d. Create a graph of observed versus simulated values. Add the diagonal
line with intercept 0 and slope of 1.
e. Create a graph of residuals. Add the horizontal line for residuals ¼ 0.

Nitrogen (kg/ha) Sugar Cane Yield (Mg/ha)

0 72
25 106
50 133
100 157
150 184
200 191

REFERENCES
Bates, D.M., Watts, D.G., 1988. Non Linear Regression Analysis and Its Applications. Wiley, New
York.
Box, G.E.P., Cox, D.R., 1964. The analysis of transformations. J. R. Stat. Soc. Ser. B (Stat
Methodol.) 26, 211–252.
Casella, G., Berger, R.L., 1990. Statistical Inference. Wadsworth and Brooks, Pacific Grove, CA.
Darroch, B.A., Baker, R.J., 1990. Grain filling in three spring wheat genotypes: statistical analysis.
Crop Sci. 30 (3), 525–529.
Goffe, W.L., Ferrier, G.D., Rogers, J., 1994. Global optimization of statistical functions with sim-
ulated annealing. J. Econ. 60, 65–99.
Graybill, F.A., 1976. Theory and Application of the Linear Model. Wadsworth and Brooks,
Belmont.
Makowski, D., Wallach, D., 2002. It pays to base parameter estimation on a realistic description of
model errors. Agronomie 22, 179–189.
 Regression Analysis, Frequentist Chapter 5 205

Mavromatis, T., Boote, K.J., Jones, J.W., Irmak, A., Shinde, D., Hoogenboom, G., 2001. Developing
genetic coefficients for crop simulation models with data from crop performance trials. Crop
Sci. 41, 40–51.
Rhinehart, R.R., 2016. Nonlinear Regression Modeling for Engineering Applications. John Wiley
and Sons, Ltd., Chichester, UK. [Link]
Robert, N., Huet, S., Hennequet, C., Bouvier, A., 1999. Methodology for choosing a model for wheat
kernel growth. Agronomie 19, 405–417.
Schabenberger, O., Pierce, F.J., 2002. Contemporary Statistical Models for the Plant and Soil
Sciences. CRC Press, Boca Raton.
Seber, G.A.F., Wild, C.J., 1989. Nonlinear Regression. Wiley, New York.
Wallach, D., Buis, S., Lecharpentier, P., Bourges, J., Clastre, P., Launay, M., et al., 2011. A package
of parameter estimation methods and implementation for the STICS crop-soil model. Environ.
Model. Softw. 26, 386–394.
White, H., 1981. Consequences and detection of misspecified nonlinear regression models. J. Am.
Stat. Assoc. (374), 419–433.