Regression

Lecture notes
Spring 2016

by Prof. Nicolai Meinshausen

Original version by Prof. Hansruedi Künsch

Seminar for Statistics

ETH Zurich

February 2016
2
Contents

1 Linear Regression 5
1.1 Introduction: underlying question . . . . . . . . . . . . . . . . . . 5
1.1.1 Examples and historical background . . . . . . . . . . . . 5
1.1.2 Linear model and examples . . . . . . . . . . . . . . . . . 7
1.2 Prerequisites for the linear model . . . . . . . . . . . . . . . . . . 12
1.3 The least squares estimator . . . . . . . . . . . . . . . . . . . . . 15
1.3.1 Normal equations . . . . . . . . . . . . . . . . . . . . . . . 15
1.3.2 Geometric interpretation . . . . . . . . . . . . . . . . . . . 17
1.3.3 Link to the MLE for the normal distribution . . . . . . . 19
1.3.4 Why not regress on each variable individually ? . . . . . . 19
1.4 Properties of LS estimation . . . . . . . . . . . . . . . . . . . . . 20
1.4.1 Moments without any normality assumptions . . . . . . . 21
1.4.2 Distribution under the assumption of normality . . . . . . 22
1.4.3 Asymptotic normality . . . . . . . . . . . . . . . . . . . . 23
1.5 Tests and confidence intervals . . . . . . . . . . . . . . . . . . . . 24
1.5.1 Basic test statistics . . . . . . . . . . . . . . . . . . . . . . 24
1.5.2 Confidence band for the entire hyperplane . . . . . . . . . 26
1.5.3 Comparison of nested models, analysis of variance . . . . 27
1.5.4 Coefficient of determination . . . . . . . . . . . . . . . . . 30
1.6 Simple linear regression . . . . . . . . . . . . . . . . . . . . . . . 31
1.6.1 Results for the special case of simple linear regression . . 31
1.6.2 Regression and correlation . . . . . . . . . . . . . . . . . . 32
1.6.3 Switching X and Y; regression to the mean . . . . . . . . 34
1.7 Residual analysis, verification of models . . . . . . . . . . . . . . 36
1.7.1 Normal plot . . . . . . . . . . . . . . . . . . . . . . . . . . 36
1.7.2 Tukey-Anscombe plot . . . . . . . . . . . . . . . . . . . . 39
1.7.3 Time series plot, Durbin-Watson test . . . . . . . . . . . . 40
1.7.4 Interior analysis . . . . . . . . . . . . . . . . . . . . . . . 41
1.7.5 Generalized least squares, weighted regression . . . . . . . 43
1.8 Model selection . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
1.8.1 Model selection using “stepwise regression” . . . . . . . . 45
1.8.2 Model selection criteria . . . . . . . . . . . . . . . . . . . 46
1.9 The Gauss-Markov theorem . . . . . . . . . . . . . . . . . . . . . 50

2 Nonlinear and nonparametric methods 53

2.1 Robust methods . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

3
4 CONTENTS

2.1.1 Influence of individual observations on the LSE . . . . . . 53

2.1.2 Huber and L1 regression . . . . . . . . . . . . . . . . . . . 55
2.1.3 Regression estimators with restrictions on influence . . . . 57
2.1.4 Regression estimators with high breaking point . . . . . . 58
2.2 Nonlinear least squares . . . . . . . . . . . . . . . . . . . . . . . . 59
2.2.1 Asymptotic confidence intervals and tests . . . . . . . . . 61
2.2.2 More precise tests and confidence intervals . . . . . . . . 62
2.3 Generalized linear models . . . . . . . . . . . . . . . . . . . . . . 64
2.3.1 Logistic regression . . . . . . . . . . . . . . . . . . . . . . 64
2.3.2 General case . . . . . . . . . . . . . . . . . . . . . . . . . 65
2.4 Cox regression . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
2.5 Nonparametric regression . . . . . . . . . . . . . . . . . . . . . . 68
2.5.1 Some procedures for the one-dimensional case . . . . . . . 68
2.5.2 Bias-variance tradeoff . . . . . . . . . . . . . . . . . . . . 70
2.5.3 Curse of dimensionality . . . . . . . . . . . . . . . . . . . 71

A Results from probability theory 1

A.1 Computation of moments . . . . . . . . . . . . . . . . . . . . . . 1
A.2 The normal distribution . . . . . . . . . . . . . . . . . . . . . . . 2
A.2.1 Univariate normal distribution . . . . . . . . . . . . . . . 2
A.2.2 Multivariate normal distribution . . . . . . . . . . . . . . 5
A.2.3 Chi-squared, t and F distributions . . . . . . . . . . . . . 7

B Literature 9
Chapter 1

Linear Regression

1.1 Introduction: underlying question

1.1.1 Examples and historical background

Before modern times a widely-held but short-sighted assumption was that vari-
ations in repeated measurement must be due to one measurement being entirely
correct and all others entirely wrong. If for example the position of a star
is measured five times and this yields five different results, it was thought that
the measurement had only been performed correctly once, and incorrectly on
the other four occasions.
This way of thinking only changed about 400 years ago, when the concept of a
“random error” was introduced. This new concept arose from the idea that
all measured data contain a small error that makes them deviate from the
truth, but that they all are approximately true. Thus stochastics was born.
This opened up the possibility of investigating approximate, stochastic relation-
ships between variables – which is the subject of this course.
The method of least squares (or LS for short) was published by Legendre in
1805. Gauss also discussed this method in a book published in 1809 – and there
mentioned that he had been using this method since 1795. This statement
cannot be proven, and so it is not clear who has the honour of having first
discovered least squares.
This method was originally used to solve problems of celestial mechanics, where
data were fit to orbits determined by theory.
Astronomical example (Ceres): On the basis of the Titius-Bode law (dis-
covered by Titius in 1766), which provides an empirical description of the pat-
tern of planetary distances from the Sun, it was thought that the space between
Mars and Jupiter must be inhabited by a further planet. On Jan 1st, 1801,
Giuseppe Piazzi found the missing body and named it Ceres. It is the largest
asteroid.

5
6 Linear Regression

Ceres moves quite rapidly and was soon lost once more. It subsequently was
Gauss who used the method of least squares to compute a sufficiently exact
orbit from the few available data, so that Ceres could be rediscovered.
Later on, this method was used in a much more general way, including in the
social sciences. Yule, for example, carried out an investigation in 1899 into
whether poor people were best served by being put into poorhouses or by being
supported in their usual surroundings. For this purpose he used the regression
equation

∆P aup = a + b · ∆Out + c · ∆Old + d · ∆P op + error.

Here ∆ denotes the changes between two successive censuses, “Paup” is the
number of people receiving poor relief, “Out” is the ratio of people getting
poor relief outside of poorhouses to the number of people in poorhouses, “Old”
is the proportion of over-65s in the general population and “Pop” is the total
population. This regression equation was fitted using data from two censuses
for each of a number of administrative districts (so-called “unions”). The dis-
tricts had largely autonomous social policies. Yule formed 4 categories of such
“unions” (rural, mixed, urban and metropolitan), and for each category the
coefficients a, b, c and d were estimated separately. Thus the censuses of 1871
und 1881 yielded a estimated coefficient of b = 0.755 for metropolitan districts.
In other words, an increase in the variable Out goes hand in hand with an in-
crease in the number of poor people – even when other influential factors such
as “Old” are accounted for. This led Yule to the conclusion that keeping people
in their usual surroundings and supporting them there actually leads to even
more poverty.
This is not a conclusive argument, however. It has been shown that those
districts with more efficient administration also built more poorhouses at that
time. At the same time, efficient administration leads to a reduction in poverty.
i.e. the effects of efficient administration and of the establishment of poorhouses
cannot be separated. We call these two variables “confounded”. Other economic
factors are also potential confounders. Generally speaking, a lot of care must
be taken when attempting to conclude causality from an observed association.
In particular, we cannot draw conclusions about interventions (changes in the
variables) on the basis of the regression equation. For further discussion of this
point, I refer to D. Freedman’s article “From association to causation: Some
remarks on the history of statistics”, Statistical Science 14 (1999), 243-258,
from which the example above is taken.
Another (fictitious) example: is there a link between the number of storks and
the rate of human births? Some corresponding data are given in Figure 1.1.
The statistics show a highly significant connection between the number of storks
and the birth rate. From this one might wrongly conclude that babies are
brought by storks. (“Cause and effect, causal link”).
In this case the confounding variable is time. This happens quite often (meaningless
correlation of time series). In this example it is quite obvious that this asso-
1.1 Introduction: underlying question 7

Figure 1.1: Connection between birth rates and stork numbers

ciation does not imply a causal relationship. However, if we look at the amount
of brutality in TV shows and the rate of violent crime, we are quite sure to find
a statistical link (both increase with time), and a causal relationship is possible
a priori – though such a relationship is not easy to prove!

1.1.2 Linear model and examples

Multiple regression:
Given: a single dependent variable (target variable) which up to
measurement errors (or random fluctuations) depends on several
“independent” or “explanatory” variables (or experimental conditions).
Wanted: the parameter values that describe this linear dependence,
and the error variance.

The same model expressed in formulae:

Yi = θ1 xi1 + . . . + θp xip + εi (i = 1, . . . , n)

Terminology:
• The numbers (Yi ; i = 1, . . . , n) form the vector y of realizations of the
dependent variable (also known as the “target variable” or “response”).
• The numbers (xij ; i = 1, . . . , n) form the vector x(j) of realizations of the
j-th independent (explanatory) variable (experimental condition)
(j = 1, . . . , p).
• The (xij ; j = 1, . . . , p) form the vector xi of explanatory variables (exper-
imental conditions) of the i-th observation (i = 1, . . . , n).
• The numbers (θj ; j = 1, . . . , p) form the unknown parameter vector θ.
• The numbers (εi ; i = 1, . . . , n) form the vector ε of (unknown) errors,
which we shall assume to be random.
8 Linear Regression

• n is the number of observations (the sample size), while p is the number

of explanatory variables.

While the θj and εi are unknown quantities, all the xij and yi are known.

Vectorial notation of model:

Yi = xTi θ + εi (i = 1, . . . , n)

Matrix notation of model:

Y = X × θ + ε
n×1 n×p p×1 n×1

where X is an (n × p) matrix with rows xTi and columns x(j) .

The first explanatory variable is usually a constant, i.e. xi1 = 1 for all i. The
model thus contains an intercept. To give an interpretation to the parameter θ1
in such a case, we assume the errors εi to have mean zero. In other situations
such an assumption is also generally made.

We furthermore assume that there are more observations than covariates (p <
n) and that the matrix X has the maximum rank possible, p, i.e. that the p
columns of X are linearly independent. If this were not the case, the parameters
would not be identifiable (different choices of parameter may yield the same
model). Sometimes models with linear dependence in the columns are used all
the same, and the identifiability forced by auxiliary conditions.

A word about the notation: we endeavour to consistently write vectors in bold

type. On the other hand we are less consistent in distinguishing between random
variables (in upper case) and their realizations (in lower case) as we often switch
between these two interpretations, and as fixed matrices are also written in
upper case.

Examples:

(1) The location model:

 
1
 1 
p = 1, X= , θ1 = µ.
 
..
 . 
1
1.1 Introduction: underlying question 9

(2) The 2–sample model:

 
1 0

 1 0 
.. ..

 

 . . 
  
 1 0  µ1
p = 2, X= , θ= .

 0 1 
 µ2

 0 1 

 .. .. 
 . . 
0 1

Here the most common questions are “Is µ1 = µ2 plausible?”, and “How
large is their difference?”, respectively. From the introductory course,
we can already handle these using the 2–sample t-test or the 2–sample
Wilcoxon test.
Statistically speaking, the 2-sample model is the simplest one imaginable.
(In practice it is mostly simpler than the 1-sample (location) model, as
systematic and semi-systematic errors often cancel out.)
(3) One-way (simple) analysis of variance with k levels (groups)
This is a generalization of the previous example from 2 to k groups. In
this case, p = k, the parameters are the group means µj (1 ≤ j ≤ k)
and the matrix X looks similar to above. Another parametrization is also
commonly used, namely
µj = µ + α j ,
where µ is the overall mean and αj is the j-th group effect. In this setup,
the k + 1 parameters cannot all be determined and the columns of X are
linearly dependent. The identifiability
P of these parameters is forced by
introducing a condition such as αj = 0.
This type of model is treated at greater depth in the course Analysis of
Variance (“Angewandte Varianzanalyse und Versuchsplanung”).
(4) Regression through the origin: Yi = βxi + εi (i = 1, . . . n).
 
x1
 x2 
p = 1, X =  . , θ1 = β.
 
 .. 
xn

(5) Simple linear regression: Yi = α + βxi + εi (i = 1, . . . n).

 
1 x1
 1 x2   
α
p = 2, X =  . . , θ= .
 
 .. ..  β
1 xn
10 Linear Regression

(6) Quadratic regression: Yi = α + βxi + γx2i + εi (i = 1, . . . n).

1 x1 x21
 
 
 1 x2 x2  α
2 
p = 3, X= . . ..  , θ =  β .

.
 . . . .  γ
1 xn x2n

Conventional interpretation of quadratic regression: We fit a

parabola to the two-dimensional point cloud (x1 , y1 ), . . . , (xn , yn ), cf. Fig.
1.2, left side.
Alternative explanation: A fixed parabola is given by the pairs (x1 , x21 ), . . . , (xn , x2n ).
In the third dimension (that is, y), we now look for a suitable plane con-
taining it, cf. Fig. 1.2, right side.
Bottom line: The function we are fitting is quadratic in the known
covariates, but linear in the unknown coefficients (and thus it is a linear
model).

Figure 1.2: Quadratic regression (left), interpreted as multiple linear regression

(right).

(7) Power or exponential dependence:

A dependence of the form Yi = αxβi + εi , or Yi = α exp(βxi ) + εi , where
α and β are unknown parameters, does not fit our model. However, the
deterministic part does not change when we take logarithms on both sides.
This leads us to the related model
log(Yi ) = log(α) + β log(xi ) + εi .
This is an example of the general linear model with target variables log(Yi )
and  
1 log(x1 )  
 1 log(x2 )  log(α)
p=2 X=   θ= .
...  β
1 log(xn )
Inverting this transformation, we obtain
Yi = αxβi · ηi ,
1.1 Introduction: underlying question 11

where ηi = exp(εi ). In other words, the errors on the original scale are
multiplicative instead of additive. When the dependence is of power
or exponential type, multiplicative errors are usually more plausible, as
then the error size is proportional to the scale of the target variables.

Figure 1.3: Dependence of tremor size on distance for various loadings

One concrete example of this is a dataset whose target variable is the

size of the tremor caused by a controlled explosion, and the explanatory
variables are the explosive load and the distance of the explosion from
the location of the measurement. The corresponding data are plotted in
Figure 1.3. It is clearly visible that the effect of distance is non-linear.

The physical expectation is that the tremor size should be inversely pro-
portional to the squared distance. In this case, we have a power model
with a known parameter β. Figure 1.4 plots the logarithmic tremor sizes
against the logarithmic distances for a fixed load size. They have an ap-
proximately linear relationship with an roughly constant variation. One
immediate question is whether or not the slope of the linear fit really is
−2, the value postulated by theoretical physical considerations. This will
be one of the questions we shall look into during this course.

These examples have shown us two important principles:

• A model is called linear if it is linear in its parameters. The original

explanatory variables can be subjected to arbitrary transformations.

• We can often obtain a linear model simply by transforming both sides of

a deterministic relationship. However, we must then think about whether
or not the errors can plausibly be additive on the transformed scale.
12 Linear Regression

Figure 1.4: Distance and tremor size for explosions with explosive loading
3.12. Logarithmic axes are used

Our goals in regression analysis:

• Producing a best “fit”. Fitting a (hyper)plane over the covariates

through the response points and minimizing deviations from these re-
sponses. The standard procedure for this is the method of least squares,
but the deviations can also be sized up in other ways.

• Good parameter estimates. These answer the question: How does the
response change when an explanatory variable changes?

• Good predictions. These answer the question: What response can we

expect to get under a new set of experimental conditions?

• An indication of the uncertainty underlying the three previous

problems using tests and confidence intervals.

• Developing a simple and functioning model. This is usually the

result of an iterative process.

1.2 Prerequisites for the linear model

For a linear model fit using least squares to be meaningful, we have to make
certain assumptions. These are also needed for the validity of the statistical
tests and confidence intervals we shall derive. Before we list these conditions
1.2 Prerequisites for the linear model 13

in descending order of importance, let us note that the model places no prereq-
uisites on the explanatory variables. These can be continuous or discrete, and
they may be transformed and combined in an arbitrary manner. Furthermore
there is no difference in principle between the deterministic fixing of the values
of explanatory variables by the experimenter, and their being realizations of
random variables themselves. The theory derived in the following sections will
always regard the explanatory variables as deterministic. That is to say: if
these variables are random, our statements are to be understood as conditional
ones given the values of the explanatory variables.
As we shall see later on, these assumptions can partly be checked statistically.

The data are useful for gaining the information sought

1. (“representative”, “meaningful”). They are a random sample
from the population under investigation.

• If this assumption is incorrect, the whole analysis is worthless from

the outset. The data might just as well be discarded.
• To judge the usefulness of the data we require some insight into the
problem. The real underlying situation is the deciding factor – and
that cannot be decided by using statistical methods.
• Our actual target variables might not be measurable in any precise
way. How does one measure intelligence, for example? We attempt
this using tests and their evaluation, and we then define intelligence
by the result of intelligence tests. The connection of these two quan-
tities remains an unanswered question, though. A further example
of this: the state measures its citizens’ wealth by the assets declared
on their tax forms.
Transferring results from measurable quantities to the ac-
tual underlying quantities of interest is a whole problem
unto itself.
2. The regression equation is correct. That is:

E [εi ] = 0 ∀i

Specifically: no significant explanatory variables should be missing and

the relationship between the target variable and the explanatory variables
should be a linear one (after suitable transformations).

3. The errors are uncorrelated. This means (under Assumption 2):

E [εi εj ] = 0 ∀i, j (i 6= j)

If the errors are correlated, the least squares fit is still of some use, but
its precision is not what we think it is. The levels we obtain for tests,
and the confidence intervals we compute, are wrong. We will discuss this
in more detail later.
14 Linear Regression

4. All covariates xi are exact.

That is, the covariates xi are all known without any errors. If, however,
they do contain measurement or observational errors, the method of least
squares makes systematic errors. There are ways of correcting for this
if at least the ratio of error variances is known for each variable. In
the literature, this is discussed under the heading of “errors in variables
models”.
At first this may seem to contradict the previous statement that the xi
may be random. However, the assumption here is that we know exactly
which value of xi led to the corresponding yi . This is a different question
to how xi came about.
5. The error variance is constant (“homoscedasticity”). That is:

E ε2i = σ 2
 
∀i

All measurements should have the same precision. (In particular, there
should be no “bad errors” with a much higher variance.) A constant error
variance can often be reached by a simple transformation of the target
variable. If the homoscedasticity assumption does not hold, the method
of least squares quickly becomes imprecise (compared to other methods).
This will also be discussed at greater depth later.

6. The errors (εi ; i = 1, . . . n) follow a joint normal distribution.

(The same then also holds for the Yi .) Such a normal distribution of
the errors is often a consequence of the general properties of the normal
distribution (cf. Central Limit Theorem in the Appendix). It should
however not be assumed without question.
Assumptions 2, 3, 5 and 6 can partly be checked using statistical tools. We
shall discuss suitable methods for this later in these notes.
The conditions listed above are usually only satisfied in an approximate way.
The art of statistics is to develop a feel for which deviations from the assump-
tions are significant, and which statements and methods are still meaningful
when the model is false.
If for instance (X1 , X2 , . . . Xp , Y ) is a (p + 1)-dimensional random vector that
follows an arbitrary distribution, and we fit a linear model using n indepen-
dent realizations of that random vector, then we are actually estimating the
coefficients of the best linear predictor, defined as
 
Xp
arg min E (Y − θ0 − θj Xj )2  .
θ0 ,...θp
j=1

Thus least squares can nearly always be used if prediction is our only goal. Some
care must be taken in interpreting the parameters and specifying the precision.
1.3 The least squares estimator 15

Finally, here is an example where Assumption 3 is not true (and neither are
Assumptions 1 and 2). The dependent variable here is the number of live
births in Switzerland since 1930, and the explanatory variable is time (as well
as certain functions of time, if e.g. quadratic trends are also of interest).

Figure 1.5: Effect of the contraceptive pill

As we see in Figure 1.5, the data on live births in Switzerland following World
War II have an approximately linear trend up to 1964. A closer look, though,
reveals that the data are not symmetrically distributed around the regression
line, but that they form “groups” on either side of it; successive maxima and
minima are about twenty years apart (one generation).
Furthermore, the years before 1964 are not “representative” for the following
years – that is, the model is no longer valid. It is generally quite dangerous
to extrapolate a fitted linear model to an area where no observations of the
explanatory variables are available.

1.3 The least squares estimator

Consider the following model: Y = Xθ + ε

We would like the “best possible” estimate of θ. The least squares estimate θ
b
is defined as the quantity that minimizes the L2 –norm of the error:

ky − X θk
b = min ky − Xθk.
θ

Thus we minimize the Euclidean distance of the error y − Xθ from the vector
zero.

1.3.1 Normal equations

To compute the least squares estimate, we calculate the partial derivatives of

ky − Xθk2 by θ (which form a vector) and require them to be zero, thus
16 Linear Regression

obtaining the equation

(−2) X T (y − X θ)
b =0 ((p × 1) − dimensional zero vector),

which is the same as: XT Xθ

b = XT y .

These are the normal equations. We have p linear equations for p unknowns
(note that X T X is a p × p matrix). The entries of X T X are the scalar products
of the columns of X. Thus solving the normal equations is especially simple
when the columns x(j) of P X are orthogonal. Another interpretation can be
found by writing X X as ni=1 xi xTi : this way, we see that X T X is n times
T

the matrix of empirical second moments of the covariates (xi ).

If we now assume that the matrix X has full rank (which means rank p), then
X T X is invertible. In such a case the least squares solution is unique and can
be written as
b = (X T X)−1 X T y.
θ
This is a useful formula for theoretical considerations, but its susceptibility to
rounding errors make it unsuitable for numerical computation. One numerically
stable algorithm uses QR decomposition and Givens rotations.
Special case: simple linear regression yi = α + βxi + εi . Then
   
1 x1 y1
 1 x2   
α  y2 
X= . .  θ= y =  . ,
   
 .. ..  β  .. 
1 xn yn

and plugging this into the normal equations gives us

Pn
= Pni=1 yi
P
·β
Pn+ ( i=1 xi )P
nα
.
( i=1 xi )α + ( ni=1 x2i ) · β = n
i=1 xi yi

To solve this system of equations in a simple manner we employ the “orthogo-

nalization” technique, i.e. we introduce the new variable
n
1X
x −→ x
e := x − x (where x := xi is the “arithmetic mean”)
n
i=1

e = α + βx, we then get yi = α

Writing α e + βe
xi . As
n
X n
X n
X
x
ei = (xi − x) = xi − nx = nx − nx = 0,
i=1 i=1 i=1

we immediately obtain:
Pn Pn
i=1 yi x
ei yi
α
e=
b =y ; β = Pi=1
b
n .
n e2i
i=1 x
1.3 The least squares estimator 17

A simple transformation back then finally gives us the desired quantities:

Pn Pn
(x i − x)yi (x − x)(yi − y)
i=1
βb = Pn 2
Pn i
= i=1 2
; b = y − βx
α b
i=1 (xi − x) i=1 (xi − x)

In a multiple linear regression setup that includes an intercept term (thus mak-
ing the corresponding column in X contain only the number 1), we can orthog-
onalize in a similar way:
X p
yi = α + θj x
eij .
j=2

(for j > 1) and α = θ1 + pj=2 θj xj .

Pn P
eij = xij − xj , with xj =
Here x i=1 xij /n
From this we can conclude that y = θb1 + pj=2 θbj xj , i.e. the point (x2 , . . . xp , y)
P
lies on the fitted plane.

1.3.2 Geometric interpretation

We can interpret this estimation procedure by looking at rows or at columns.

In the former case, if the model contains an intercept (i.e. xi1 ≡ 1), then we
have n points (yi , xi2 , . . . xip ) randomly spread around a (p − 1)-dimensional
hyperplane in a p-dimensional space. (If our model lacks an intercept, we have
n points spread around a hyperplane through the origin in (p + 1)-dimensional
space.) However, this random spread only occurs parallel to the y-axis. The
least squares estimator therefore estimates the parameters of the hyperplane
so as to minimize the sum of squared distances of the points from the plane
parallel to the y-axis.
More mileage is to be had by interpreting the column vectors in the model.
Here the vector y of observations is a single point in the n-dimensional space
Rn . If we vary the value of the parameter θ, the product Xθ describes a
p-dimensional subspace of Rn , i.e. a p-dimensional hyperplane through the
origin. Then an obvious way to estimate θ is to make Xθ minimize y on
this hyperplane. Choosing the L2 norm to yield our metric in Rn amounts
to choosing the Euclidean distance, which geometrically implies that we are
performing an orthogonal projection of y onto this hyperplane. In particular,
the least squares estimator is characterized by the property that r = y − X θb
(the vector of residuals) is orthogonal to all the columns of X:

b T X = 0.
(y − X θ)

This amounts to a geometric interpretation of the normal equations (cf. Figure

1.6).
Now the orthogonal projection X θ b is the estimate of E [y] given by our model.
Its components are called the fitted values, and they are generally denoted by
y
b.
18 Linear Regression

Figure 1.6: The vector r of residuals is orthogonal to the hyperplane spanned

by the vectors 1 and x.

To compute y
b , we can use the formula
y
b = Xθb = X(X T X)−1 X T y = P y and thus: y
b = Py
| {z }
=:P

It is easy to check that the matrix P has the following properties:

P T = P, P2 = P
and
X
Pii = tr(P ) = tr(X(X T X)−1 X T ) = tr((X T X)−1 X T X) = tr(Ip×p ) = p.
i
These are necessary and sufficient conditions for P to be an orthogonal projec-
tion from Rn to Rp .
It is quite evident that the matrix P only depends on the explanatory variables
(the experimental conditions), but not on the target variables (the responses).
It is also known as the hat matrix, as it “puts the hat on” y. Furthermore, the
diagonal entries Pii tell us how much influence the observation yi (at the point
xi ) has over the fitted value ybi .
The residuals r, also denoted by b
ε, can be written in a manner similar to the
fitted values. That is:
ε=y−y
r=b b = (I − P ) y = Qy and thus: r = Qy
| {z }
=:Q
1.3 The least squares estimator 19

Q is an orthogonal projection yet again; it is orthogonal to P and has values in

the residual space:

QT = Q2 = Q, P Q = QP = 0 (0 as an n × n matrix), tr(Q) = n − p.

1.3.3 Link to the MLE for the normal distribution

Assuming that the explanatory variables are given, the assumptions of our
model (independence and normality of the errors) imply that the conditional
density of y1 , ..., yn is
n p
2
Y 1 X
Ly,X (θ, σ ) = ϕ((yi − θj xij )/σ).
σ
i=1 j=1

When performing maximum likelihood estimation, we consider the density to

be a function of the parameters σ and θ (fixing yi and xij ). This function
is the so-called likelihood function, and we estimate the parameters θ and σ
so that they maximize this likelihood function (or equivalently, we maximize
its logarithm). The result of this process is called the Maximum Likelihood
Estimator (MLE).
It is immediately obvious that maximizing with respect to θ does not depend
on the value of σ, and that it amounts to minimizing ky − Xθk2 . If all the
errors εi are i.i.d. and they all have the distribution N (0, σ 2 ), the MLE of θ is
exactly the same as the least squares estimator.
The MLE for σ 2 , on the other hand, is
Pn
2 (yi − ybi )2
σbM L = i=1 .
n

This is not the estimator generally used, however. Instead, it is scaled to become
unbiased. We shall see that the correct scaling factor here is n/(n − p). In other
words, the estimator of error variance we shall use is
Pn
2 (yi − ybi )2
b = i=1
σ .
n−p

1.3.4 Why not regress on each variable individually ?

The following (artificial) example shows why multiple regression cannot simply
be replaced by several simple regression procedures.
Let there be 2 covariates x1 , x2 , and assume that we have the following obser-
vations:
x1 0 1 2 3 0 1 2 3
x2 −1 0 1 2 1 2 3 4
y 1 2 3 4 −1 0 1 2
20 Linear Regression

Figure 1.7: Multiple regression versus simple regression

Left side of Fig. 1.7: We plot the values of y against the corresponding
values of the covariates x1 und x2 . Subsequently, we find a plane that fits these
8 points exactly (in three-dimensional space):

y = 2x1 − x2 σ 2 = 0)
(b

The coefficients (2 and -1, respectively) tell us how y changes if we change

exactly one of x1 or x2 and fix the other one.
We conclude that y decreases as x2 increases (x2 larger ⇒ y smaller).
Right side of Fig. 1.7: We simply regress y on x2 and forget about the values
of x1 (they are not kept fixed). The ensuing regression line is:

1 4
y = x2 + σ 2 = 1.72)
(b
9 3

This line tells us how changes x2 affect y if x1 is allowed to vary.

We conclude that y increases as x2 increases (x2 larger ⇒ y larger).
The reason for this difference in the behaviour of y depending on x2 is that the
covariates x1 and x2 are strongly correlated. That is: when x2 increases, so
does x1 .
In summary:

Combining several simple regressions (each using the method of least

squares) generally only gives us the same result as a multiple regression
if the explanatory variables are orthogonal.

1.4 Properties of LS estimation

Let us first take an intuitive look at the precision of the regression plane. We
assume a known linear model and simulate a random point cloud from this
model. This cloud of points we can now use to fit a regression plane using the
1.4 Properties of LS estimation 21

method of least squares. If we now take a second (or third, . . .) point cloud,
we will generally get a different estimated regression plane – even though the
underlying model, and thus the theoretical plane, are the same (cf. Fig. 1.8).
In other words, the parameter estimates and the fitted regression plane are
random ! Because of this, we do need some idea of their precision.

Figure 1.8: Three different estimated regression lines for the same underlying
model.

1.4.1 Moments without any normality assumptions

The next results do not require the errors εi to be normally distributed. The
assumptions in this section are the following:

E [ε] = 0 
Usual model: Y = Xθ + ε with
Cov[ε] = E εεT = σ 2 In×n


Results:
h i
(i) E θb = θ, as
h i
b = E (X T X)−1 X T y = E (X T X)−1 X T (Xθ + ε) = θ + 0 = θ.
   
E θ

(ii) E [b
ε] = 0, E [b
y] = E [y] = Xθ.
(iii) Cov[ 2 T −1
h θ] = σ (X X)i , as 
b
b − θ)T = E (X T X)−1 X T εεT X(X T X)−1 = σ 2 (X T X)−1 .

b − θ)(θ
E (θ

y] = Cov[P y] = σ 2 P P T = σ 2 P
(iv) Cov[b (because P is a projection matrix).
ε] = σ 2 Q
(v) Cov[b (similarly).
(vi) Cov[b
ε, y
b ] = 0, (as QP = 0).
The covariance matrices in (iv) and (v) are only positive definite. We can see
from (v) that the residuals ri = εbi are correlated – unlike the true errors – and
22 Linear Regression

that their variance

εi ] = σ 2 (1 − Pii )
Var[b
is not constant. We can furthermore conclude that
" n # n
X X
2 2
E ri = σ (1 − Pii )
i=1 i=1
= σ (n − tr(P )) = σ 2 (n − p).
2

Therefore Pn 2 b 2
2 i=1 ri ky − X θk
σ
b = =
n−p n−p
is an unbiased estimate of σ2, as we already claimed.
b2 . For this
Note that we cannot make any statements about the variance of σ
we would need to know the fourth moment of the errors εi .

1.4.2 Distribution under the assumption of normality

Assumptions in this section:

Usual model: Y = Xθ + ε now assuming ε ∼ N n (0, σ 2 In×n )

Results:
b ∼ N p (θ, σ 2 (X T X)−1 ) (as θ
(i) θ b is a linear combination of normally
distributed quantities and therefore itself follows a normal distribution).
b ∼ N n (Xθ, σ 2 P ),
(ii) y ε ∼ N n (0, σ 2 Q)
b (for the same reason as above).
iii) y
b and bε are independent (as they are uncorrelated and both normally
distributed).
(iv) Pn 2
i=1 ri
∼ χ2n−p .
σ2
(see below).
b2 is independent of θ
(v) σ b = (X T X)−1 X T y (This is a consequence of iii),
b = (X T X)−1 X T y
as θ b ).
Proof of iv): Regard a coordinate system with an orthogonal basis, such
that the first p vectors in the basis span the column space of X. Denote the
corresponding transformation matrix by A; in other words, the columns of A
contain the coordinates of the new basis vectors in the old coordinate system.
Then A is orthogonal, and if we use stars to denote the new coordinate system,
we have y∗ = AT y, ε∗ = AT ε etc.. By construction,

b ∗ = (y1∗ , y2∗ , . . . yp∗ , 0, . . . 0)T ,

y
ε∗ = (0, . . . , 0, ε∗p+1 , . . . ε∗n )T
b
1.4 Properties of LS estimation 23

(which can easily be checked, noting that the last n − p rows of AT X are all
identically equal to zero). Because of the orthogonality of A, we find that
n
X n
X n
X
εb2i = εb∗2
i = ε∗2
i .
i=1 i=1 i=p+1

From this the claim follows, as the orthogonality of A means that ε∗ also follow
a normal distribution N n (0, σ 2 In×n ).

1.4.3 Asymptotic normality

The above results on the distribution of estimators are the key to subsequent
statements about uncertainty, i.e. confidence intervals or tests. One might
therefore ask how decisive the assumption of normal errors is. It has been seen
that the results are still approximately correct if the errors are not normally dis-
tributed. This is investigated mathematically by looking at the limit behaviour
of the distribution when the number of observations goes to infinity.
Consider the following situation: We have n data (y1 , x1 ), . . . , (yn , xn ) which
satisfy the linear model. Here each xi is a p-dimensional column vector (i.e. xTi
is the i-th row of X). We assume that the errors εi are i.i.d. but not necessarily
normally distributed, and look at the limit case n −→ ∞.
For the asymptotic approximation to hold, we need some weak conditions on
the explanatory variables xi :
• The smallest eigenvalue of X T X = ni=1 xi xTi , namely λmin,n , converges
P
to ∞.
• maxj Pjj = maxj xTj ( ni=1 xi xTi )−1 xj converges to zero.
P

The first condition states that increasing n always yields more information,
while the second condition prohibits any xj from dominating the others.
Theorem 1.4.1. If the errors εi are i.i.d. with mean 0 and variance σ 2 , and if
(xi ) satisfies the conditions just given, then the LS estimators θ
b are consistent
(for θ), and the distribution of

(X T X)1/2 (θ
b − θ)

converges weakly to N p (0, σ 2 I).

Proof: The i-th component θbi is unbiased and has variance σ 2 ((X T X)−1 )ii ,
which converges to zero by the first assumption. Then consistency follows from
Chebyshev’s inequality.
To show the weak convergence of a random vector in Rp , it suffices to show the
weak convergence of the distributions of all its linear combinations (Theorem
of Cramér and Wold, for which see the literature). So we consider
n
X
b − θ) = cT ((X T X)−1/2 )X T ε = an T ε =
cT ((X T X)1/2 )(θ ani εi ,
i=1
24 Linear Regression

where
an = X((X T X)−1/2 )T c.
We have a sum of n independent, but not identically distributed terms ani εi .
Furthermore, the distribution of these summands changes with n. This is pre-
cisely the situation for which Lindeberg’s Theorem was established (see the
introductory class for details). We have
n
X n
X
Var ani εi = σ 2 a2ni = σ 2 an T an = σ 2 cT c.
i=1 i=1

Without loss of generality, we can assume this variance to be 1. We now need

only check the condition
n
X
a2ni E ε2i 1[|εi |>η/ani ] → 0
 

i=1

for all η > 0. As all εi have the same distribution and a finite second moment,
it follows that  d→∞
E ε2i 1[|εi |>d] = E ε21 1[|ε1 |>d] −→ 0.
  

As i a2ni = 1 also holds, it suffices to show that maxi |ani | converges to zero.
P
From the Schwarz inequality, we obtain

a2ni ≤ kck2 xTi (X T X)−1 xi .

Thus the claim follows straight from the second condition. 2

We can also show that σ b2 is consistent. However, to show the asymptotic
normality of σ 2
b we require the existence of a fourth moment for the εi . The
asymptotic variance would then depend quite strongly on the exact value of
this variance.
Conclusions:
The tests and confidence intervals for θ and for the means E [y], which we shall
derive using the assumption of normal errors, still have roughly the correct level
if this normality assumption does not hold. However, we know nothing about
the efficiency of these methods in such a case, and the confidence intervals for σ
found in the literature usually have a very wrong level in non-normal situations.

1.5 Tests and confidence intervals

1.5.1 Basic test statistics

Assume the linear model with ε ∼ N n (0, σ 2 In×n ). As seen in the last sec-
b exactly follows an N p (θ, σ 2 (X T X)−1 )–
tion, these assumptions imply that θ
distribution, and σb2 is independent of θ.b Some of the consequences of this
are:
1.5 Tests and confidence intervals 25

(a) For each individual parameter θi , we have

θb − θi
p i ∼ tn−p .
b ((X T X)−1 )ii
σ

(b) The entire parameter vector θ gives rise to an F -distribution as follows:

b − θ)T (X T X)(θ
(θ b − θ)
2
∼ Fp,n−p
pb
σ

(c) Similarly for each linear transformation ϑ = Bθ by a (q × p) matrix B:

b − ϑ)T V −1 (ϑ
(ϑ b − ϑ)
∼ Fq,n−p ,
qbσ2

where V = B(X T X)−1 B T .

(d) For the expectation of the i-th observation (that is, the true location of
the hyperplane for the i–th experimental condition), we have

ybi − E [yi ]
√ ∼ tn−p where pii := (P )ii
σ
b pii

(e) For the expectation of a new observation under an arbitrary experimental

condition x0 (the true location of the hyperplane for the new experimental
condition x0 ), we have

yb − E [y0 ]
q 0 ∼ tn−p
b xT0 (X T X)−1 x0
σ

(f) A random new observation y0 = y0 (x0 ) under the experimental condition

x0 satisfies
y0 − yb0
q ∼ tn−p .
b 1 + xT0 (X T X)−1 x0
σ

These facts allow us to carry out statistical tests in the usual way (using the
quantities on the left side as test statistics) and to compute confidence intervals
for individual parameters, for linear combinations of them or for the unknown
true location of the hyperplane at some x0 . Furthermore, statement (f) enables
the construction of prediction intervals for future observations.
We illustrate this by the tremor size example that we introduced earlier. Ta-
ble 1.1 shows the computer output when we take the logarithm of tremor size as
our response variable and the logarithms of distance and load size as explana-
tory variables.
The output
p contains the estimated coefficients as well as the residual standard
b ((X T X)−1 )ii and the results of tests of the null hypotheses θi = 0.
error σ
26 Linear Regression

Coefficients:
Value Std. Error t value Pr(> |t|)
(Intercept) 2.8323 0.2229 12.71 0.000
log10(dist) -1.5107 0.1111 -13.59 0.000
log10(ladung) 0.8083 0.3042 2.66 0.011
Residual standard error: 0.1529 on 45 degrees of freedom
Multiple R-Squared: 0.8048
F-statistic: 92.79 on 2 and 45 degrees of freedom
p-value 1.11e-16

Table 1.1: Computer output for the tremor size example

The other numbers here will be explained in the subsequent sections. Now the
t distribution with 45 degrees of freedom is very close to the standard normal
distribution. Thus the numbers in Table 1.1 show us that the true coefficient
of the logarithm of distance is less than 2 (in absolute terms), while the true
coefficient of the logarithm of load size may well be 1.

1.5.2 Confidence band for the entire hyperplane

We can also compute a 95% confidence set, for instance, in which the true
hyperplane lies. Before we do this, we will first mention an obvious strategy
which fails. For any point x0 (at which the fitted hyperplane has the expected
response value yb(x0 )), we can construct an interval around yb(x0 ) as above, inside
which we expect the value on the true hyperplane to lie with 95% confidence.
If we do this for each possible value of x0 , do we then get a 95% confidence set
for the true hyperplane?

The answer is “no”, of course. At each individual x0 , the probability of the

true hyperplane passing through this confidence set is exactly 95%, but for two
such intervals, the probability of the hyperplane passing through both of them
is at least 90% and at most 95%.
Similarly in extreme cases: 10 points ⇒ 50 - 95 %
20 points ⇒ 0 - 95 %
This is not a good way to get a confidence set for the true hyperplane!

There is a better way, and it is as follows: The Schwarz inequality for the scalar
product ha, bi = aT (X T X)−1 b implies that

|b b − θ)| = |xT (X T X)−1 (X T X)(θ

y0 − E [y0 ] | = |xT0 (θ b − θ)|
0
≤ (xT (X T X)−1 x0 )1/2 ((θ
0
b − θ)T (X T X)(θ
b − θ))1/2 .

Using statement (b) above tells us that with probability 1 − α, we have

b2 (xT0 (X T X)−1 x0 ) p Fp,n−p (1 − α)

y0 − E [y0 ])2 ≤ σ
(b
1.5 Tests and confidence intervals 27

simultaneously for all x0 . This gives us the simultaneous confidence set we are
looking for. Its shape is that of a hyperboloid, and it is the envelope of all
hyperplanes whose parameters are compatible with the data according to b).

1.5.3 Comparison of nested models, analysis of variance

Prerequisites:
“Basic hypothesis” H : y = Xθ + ε
(X : n × p, rank(X) = p, ε ∼ N n (0, σ 2 I)).
“Special null hypothesis” H0 : the above, and additionally Bθ = b
(where the dimensions of B are (p − q) × p,
and rank(B) = p − q < p).
Example:
 
1 0 0 ... 0 0 ... 0

 0 1 0 ... 0 0 ... 0 

B=
 0 0 1 ... 0 0 ... 0 ,
 b = 0.
 ... ... ... ... ... ... ... ... 
0 0 0 ... 1 0 ... 0
This null hypothesis can be written in words as “The first p − q coefficients θi
are all zero.” We are testing whether the first p − q variables are superfluous
to the model.
By statement (c) of Section 1.5.1,
b − b)T (B(X T X)−1 B T )−1 (B θ
(B θ b − b)
σ2
(p − q)b
is a suitable test statistic for this null hypothesis. Under this hypothesis, its
distribution is Fp−q,n−p . However, we can use the following geometric argument
to obtain a different shape and interpretation of this test statistic.
We assume that b = 0 (this is not a significant restriction, as we can replace
the original observations by new ones y − Xθ, taking some θ that satisfies
Bθ = b). Then we can project y, first into the p–dimensional space spanned
by the columns of X, and from there into the q–dimensional subspace defined
by the additional condition Bθ = 0.
Let the corresponding sums of squares of the residuals (under H and H0 ) be
SSE and SSE0 .
We know that:
• SSE/(n − p) always is an unbiased estimator of σ 2 , assuming the basic
hypothesis H and the null hypothesis H0 .
• SSE and SSE0 − SSE are sums of squares in orthogonal subspaces, and
under the null hypothesis H0 , we know that (SSE0 − SSe )/(p − q) is
an unbiased estimate of σ 2 . If only the basic hypothesis H is true, the
expectation of this difference is greater than σ 2 .
28 Linear Regression

 1 y

  e
 
e
 e SSe (n-p)-dim.
  e

r
  e

  -e ŷ
Q  %
Q  %
Q
Q
Q SSe0  %
 %
Q
Q (n-q)-dim.  %
Q
Q  %SSe0 − SSe

 r % (p-q)-dim.
Q  %
Q 
Q
Q r  %
QQ% ŷ(0)
s

Figure 1.9: Comparison of models

• Because of the orthogonality of the subspaces, we have

2 2
b (0)
y−y b k2 + y
= ky − y b (0)
b−y

.
Thus given H0 :
2
(SSE0 − SSE)/(p − q) b−y
y b (0) /(p − q)
= ∼ Fp−q,n−p
SSE/(n − p) b k2 /(n − p)
ky − y

and we can use the expression on the left side as a test statistic for H0 .
At first glance, the two test statistics we have just derived are different: having
the same distribution does not make them identical. However, the following
lemma shows that both these expressions are in fact identical.
Lemma 1.5.1. The least squares estimator θ b(0) under the supplementary con-
dition Bθ = b is

θ b − (X T X)−1 B T (B(X T X)−1 B T )−1 (B θ

b(0) = θ b − b).

Furthermore,
b − b)T (B(X T X)−1 B T )−1 (B θ
SSE0 = SSE + (B θ b − b).

Proof: We introduce a vector λ containing the p − q Lagrange multipliers for

the p − q supplementary conditions. Our task is now to minimize

(y − Xθ)T (y − Xθ) + (Bθ − b)T λ

over θ and λ. This gives us the equations

X T (y − X θ
b(0) ) + B T λ = 0, Bθ
b(0) = b.
1.5 Tests and confidence intervals 29

It is easy to check that θ

b(0) satisfies these conditions. Moreover, Pythagoras’
Theorem implies that

b(0) )T (y − X θ
(y − X θ b(0) ))T (X(θ
b−θ
b(0) ) = SSE + (X(θ b−θ
b(0) )).

Plugging in, we obtain the second claim. 2

Just as multiple regression cannot merely be replaced by simple regressions

on individual variables, the test of the null hypothesis β1 = β2 = 0 may yield
completely different results than the two tests of the null hypotheses β1 = 0 and
β2 = 0. For instance, it may happen that the latter null hypotheses are both
accepted without any problem, but the combined null hypothesis β1 = β2 = 0
is thrown out quite clearly. This means that either of the explanatory variables
can be left out – but not both. The solution to this apparent paradox lies in
the heavy correlation of the two variables. Thus either of them can replace the
other quite simply.

It frequently happens that an explanatory variable is categorical (place of origin,

type, colour, sex, . . . ). In the tremor size example, explosions were carried out
at six different locations; this may impact the result by means of variation
in ground consistency. Such a variable is also called a factor. The simplest
model for this merely postulates a different intercept for each category and
assumes that all other coefficients are the same for all categories. This model
is written by introducing indicator variables for each category and using these
as additional explanatory variables. For the matrix X to still have full rank,
we must omit either the first column xij ≡ 1 or the indicator variable for the
first category. For such a categorical variable, one meaningful null hypothesis is
that the coefficients of all indicators are zero; this can be testes using an F -test.
The results for the tremor size example are given in Table 1.2. The third row
compares the full model to the model without the explanatory variables given
by the location factor “St”. It shows the extreme significance of the location.

Df Sum of Sq RSS F Value Pr(F)

log10(dist) 1 2.79 5.07 108 0

log10(loading) 1 0.59 2.86 23 7.62e-06
St 5 2.10 4.38 16 0
Table 1.2: Tests of the effects of individual terms in the tremor size example

The decomposition

2 2
b (0)
y−y b k2 + y
= ky − y b (0)
b−y

is also known as analysis of variance, and thus Table 1.2 is known as an Analysis
of Variance (ANOVA) table.
30 Linear Regression

1.5.4 Coefficient of determination

One particularly special case of the preceding results is the following: testing
to see whether the response actually depends on the covariates x.
   
1 x12 . . . x1p 0 1 0 ... 0
 1 ... ... ...   0 0 1 ... 0 
X=  ... ... ... ... 
 B=  . . . . . . . . . . . . . . .  (p−1)×p


1 . . . . . . xnp 0 0 0 ... 1

   
y y
 0   y 
Under the null hypothesis H0 : θ
b(0) =  b (0) = 
y  =: y (n×1)
  
 ..  ..
 .   . 
0 y

b (0) k2 = ky − yk2
SSE0 = ky − y

b k2 + kb
= ky − y y − yk2

Analysis of variance (ANOVA) table:

Sum of squares Degrees of freedom Mean square E [Mean square]

2
Regression y − yk2
kb p−1 y − yk2 /(p − 1)
kb σ 2 + ||E[y]−E[y]k
p−1
Error ky − yb k2 n−p b k2 /(n − p)
ky − y σ2
Total around
overall mean ky − yk2 n−1 − −
We test the significance of the dependence on covariates by means of the test
statistic
y − yk2 /(p − 1)
kb
F = .
b k2 /(n − p)
ky − y
Under H0 , this follows an Fp−1,n−p distribution. In the class Analysis of Vari-
ance (“Angewandte Varianzanalyse und Versuchsplanung”), similar (but more
complex) decompositions of the sum of squares ky − yk2 are analyzed in special
linear models with a range of F -tests.
We could also examine the distribution of the test statistic F under the alter-
native hypothesis H ∩ (¬H0 ), and thus study the power of the F -test. The
distribution we obtain is the so-called non-central F distribution Fp−1,n−p,δ2
with non-centrality parameter δ 2 = kE [y] − E [y] k2 (although the liter-
ature defines this non-centrality parameter in more than one way). See the
literature for more details.
One important quantity is the quotient
y − yk2
kb
R2 := .
ky − yk2
1.6 Simple linear regression 31

This is the coefficient of determination, the proportion of variance explained

by the model. It measures the goodness of fit of the model with explanatory
variables x(j) . It is not difficult to see that R2 is also the maximum squared
correlation of y with an arbitrary linear combination of the columns x(j) . The
coefficient of determination is also equal to the square of the multiple corre-
lation coefficient between y and the x(j) . The linear combination maximizing
the correlation with y is the least squares estimate y b itself.
Remark: R2 and F are at first the most important numbers in the computer
output.

1.6 Simple linear regression

1.6.1 Results for the special case of simple linear regression

We have already derived the least squares estimators. Here we merely give
explicit formulæ for the most important test statistics and confidence intervals:
Test of the null hypothesis β = β0 at level γ: Reject the null if
|βb − β0 |
√ > tn−2;1−γ/2 ,
σ
b/ SSX
where
n
X
SSX = (xi − x)2 .
i=1
Correspondingly, the confidence interval for β is
σ
βb ± tn−2;1−γ/2 · √
b
.
SSX
The confidence interval for the expectation of a new observation at x0 (i.e. the
value of the regression line at x0 ) is:
s
b 0 ± tn−2;1−γ/2 · σ 1 (x0 − x)2
α
b + βx b + .
n SSX
The confidence interval for the entire regression line (simultaneously for all x)
is s
p
b ± 2F2,n−2;1−γ · σ 1 (x − x)2
α
b + βx b + .
n SSX
Naturally, the simultaneous confidence interval is better than the individual
one. Finally, the prediction interval for a new observation at x0 :
s
b 0 ± tn−2;1−γ/2 · σ 1 (x0 − x)2
αb + βx b 1+ + .
n SSX
The prediction interval is wider than the confidence interval. All three intervals
have boundaries that form hyperbolae.
32 Linear Regression

1.6.2 Regression and correlation

The concept of correlation used to be applied more often than regression was.
Let Y and X be random variables, i.e. the data x1 , . . . , xn are no longer
considered to be fixed.
Definition 1.6.1. The correlation (“Pearson product moment correlation
coefficient’’) is defined as:
Cov(X, Y )
ρ = ρ(X, Y ) = p (if Var(X) 6= 0, Var(Y ) 6= 0)
Var(X) · Var(Y )

Properties of correlation:
(i) −1 ≤ ρ ≤ +1 (Schwarz inequality)
(ii) |ρ| = 1 ⇔ The joint distribution of X and Y is concentrated on a line
(and the sign of ρ matches the sign of this lines gradient).
(iii) If ρ = 0, X and Y are said to be uncorrelated.
(iv) ρ can be estimated by
Pn
i=1 (xi− x)(yi − y)
r = ρb = pPn Pn ,
2 2
i=1 (xi − x) · i=1 (yi − y)
and for this estimate ρb, we have:
• −1 ≤ ρb ≤ 1
• |b
ρ| = 1 ⇔ all the points lie on a single line
• sign (b
ρ) = sign (β)
b

In Figure 1.10, we can see some typical scatterplots using a variety of correlation
coefficients.

Figure 1.10: Scatterplots with various correlation coefficients.

The z transformation (“variance-stabilizing transformation for the

correlation coefficient”) (Fisher). Let (X, Y ) be jointly normally distributed.
Define  
1 1 + ρb
z := tanh−1 (b
ρ) = log .
2 1 − ρb
Then for an arbitrary value of ρ, we have the following good approximation (for
about n > 10)  
−1 1
z ∼ N tanh (ρ), .
n−3
Graphical interpretation of the z transformation:
1.6 Simple linear regression 33

• If the true value of ρ is near 0, the variance of ρb is high.

• If the true value of ρ lies near ±1, the variance of ρb is small.
The z transformation rescales so as to make the variance constant (i.e. it
“compresses in the middle” und “stretches at the edges”).
If we want to test ρ = 0 against ρ 6= 0, we have 3 tests to choose from:
1. Table or diagram (see Figure 1.11)
2. t- or F -test of β = 0
3. tanh−1 transformation.
The first and third methods also allow testing for any other fixed value of ρ
(and thus the construction of confidence intervals).

Figure 1.11: Confidence limits of the correlation coefficient. On the horizontal

axis: the correlation coefficient r of the sample; on the vertical axis: the true
correlation coefficient ρ. The labels on the curves denote the sample size (from
F.N. DAVID: Tables of the Ordinates and Probability Integral of the Distri-
bution of the Correlation Coefficient in Small Samples, The Biometrika Office,
London 1938)

Rank correlation: Since Pearson’s correlation is not robust towards outliers

(see Figure 1.12), some sort of rank correlation is often used. There are two
types: that of Spearman and that of Kendall:
“Spearman’s rank correlation” is simply the Pearson correlation of the ranks
34 Linear Regression

Figure 1.12: Various types of point clouds with correlation r = 0.7

of the Xi with those of the Yi . As the sum of ranks (or of their squares) is fixed,
the formulæ can be simplified thus:
6 ni=1 Di2
P
rS = 1 − Di := Rg(Xi ) − Rg(Yi ).
n(n2 − 1)
Kendall’s rank correlation is defined as
Tk − Td
rK = 2 · ,
n(n − 1)

Tk = # concordances = # pairs with (xi − xj )(yi − yj ) > 0

where
Td = # discordances = # pairs with (xi − xj )(yi − yj ) < 0
Addendum: partial correlations
Let X, Y and Z be random variables. Then the partial correlation between X
and Y given Z is defined as:
ρXY − ρXZ ρY Z rXY − rXZ rY Z
ρXY.Z := q , or estimated as: rXY.Z := q
(1 − ρ2XZ )(1 − ρ2Y Z ) 2 )(1 − r 2 )
(1 − rXZ YZ

This measures the strength and direction of the linear dependence between X
and Y after accounting for the linear dependence of X and Y on Z.

1.6.3 Switching X and Y; regression to the mean

If both X and Y are considered to be random, we can write the least squares
estimation line as follows:
σ
bY
y − y = ρb (x − x).
σ
bX
1.6 Simple linear regression 35

Since ρb is always less than 1 in absolute terms, the prediction for Y is always
closer to the mean than the corresponding value of X, if we measure distances in
terms of standard deviations. If for instance ρb is positive and our observation
of X lies e.g. 1 standard deviation above the mean, the prediction of the
corresponding value of Y will be less than one standard deviation above the
mean. In other words, we always predict a regression (return) to the mean; this
has given rise to the name “regression”.
This phenomenon is continually being rediscovered, and is often interpreted at
length in the framework of cultural pessimism. However, as our formulae show,
this phenomenon is a very general one that occurs continually and requires
no special interpretation. It lies in the very nature of prognoses that as more
observations are available in the middle, any prognosis tends towards the mean.

Figure 1.13: Regression lines “y versus x” und “x versus y”

Thus such a regression to the mean is nothing special. This becomes even
clearer if we switch the roles of X and Y . Due to symmetry, the regression line
for X by Y looks like
σ
bX
x − x = ρb (y − y).
σ
bY
In other words, if we look back and ask what value of X led to Y being one
standard deviation above the mean, the answer will be “less than one standard
deviation”. Now we might be tempted to see this as a sign of progress rather
than regression!
Drawing both regression lines in the same plot, we can see how the gradients

σ
bY 1σ bY
ρb and
σ
bX ρb σ
bX

differ. These two regression lines are quite evidently not the same; instead, we
have the scissor-like gap shown in Figure 1.13).
36 Linear Regression

If (X, Y ) follows a two-dimensional normal distribution, the contours of the joint

density are ellipses. The regression lines intersect a contour at points in which
it has vertical or horizontal tangents, respectively, as these points are maxima
of the conditional density of X given Y = y or Y given X = x, respectively.

1.7 Analyzing residuals, verifying the model and deal-

ing with breaches of assumptions

Residual analysis is the process of graphically (and some times also numerically)
analyzing the residuals, i.e. the error estimates

ri := εbi = yi − ybi ,

in order to verify the assumptions on the model after fitting it, and to develop
a better model.

1.7.1 Normal plot

Assumptions on the distribution can quite generally be checked with a quan-

tile/quantile plot (QQ plot). If we are specifically checking against the normal
distribution, we call it a normal plot.
First we introduce the normal plot for i.i.d. random variables X1 , . . . , Xn . The
“empirical cumulative distribution function” is defined as

1
u = Fn (x) = #{Xi ≤ x}.
n
This is a step function which approaches the true distribution function when n
becomes large (Lemma of Glivenko and Cantelli). In particular, we have

x−µ
Fn (x) −→ Φ( )
σ
if the Xi follow a normal distribution. Thus if we set

z := Φ−1 (Fn (x)),

then we obtain
x−µ
z≈
σ
for sufficiently large n. In a normal plot, we plot x against z at selected points. If
Xi really does follow a normal distribution, the normal plot will roughly exhibit
a straight line whose intercept and slope are µ and σ, respectively. However,
the random fluctuations of the data do lead to some deviation away from an
exact line. We can get an idea of the size of such fluctuations by performing
simulations – see Figure 1.14.
1.7 Residual analysis, verification of models 37

Figure 1.14: QQ plots for a normally distributed random variable X, from

which a) 20, b) 100 and c) 1000 samples are taken.

If the assumption of normality is incorrect, the normal plot shows systematic

deviations from a straight line. Some typical cases are given in Figure 1.15.
However, the interpretation of this plot is not always clear, as the boundaries
between the different types of situations are continuous. For instance, the case
illustrated in Figure 1.16 can be interpreted either as a mixture of two groups
or as a short-tailed distribution.

Figure 1.15: a) Heavy-tailed distribution b) Skew distribution c) Outliers

There is also a formal test of normality based on the normal plot, the “Shapiro-
Wilks test”. It essentially measures the correlation of the point cloud seen in
the normal plot.
38 Linear Regression

Figure 1.16: Two groups or short-tailed distribution?

Making these plots

• For small n (ca. n ≤ 100):

Plot all observations individually along the vertical axis. This automat-
ically gives us the order statistics (ordered observations) X(1) ≤ X(2) ≤
. . . ≤ X(n) . By definition, we have Fn (X(i) ) = i/n, but the effects of
the jumps of Fn are generally toned down by plotting a slightly modi-
fied quantity on the horizontal axis, such as Φ−1 ( i−1/2 −1 i
n ) (or Φ ( n+1 ),
i−(3/8) i−(1/3)
Φ−1 ( n+(1/4) ) or Φ−1 ( n+(1/3) ).

• For large n (ca. n ≥ 100):

Choose some equidistant values of x from the range of the sample (hori-
zontal axis).

• The horizontal axis is frequently labelled u = Φ(z). This non-linear scale

ensures that the plot of (u, x) is a line. Thus the points Fn−1 (u) on the
vertical axis can easily be plotted without any calculations. In the time
before computers (when special “probability graph paper” was used, this
was especially important.

• The axes are often switched around, i.e. z is plotted against x.

We have so far not paid any attention to the errors εi while performing regres-
sion. However, √ we can use the normal plot of the residuals εbi or the standardized
residuals εbi / 1 − Pii . Recall that the residuals εbi neither have constant vari-
ance, nor are they uncorrelated. Standardizing makes their variance constant
again, though this effect is usually negligible.

In summary: the normal plot checks the normality of the residuals against
possible skewness, heavy- (or short-)tailedness, outliers and other peculiarities.
1.7 Residual analysis, verification of models 39

1.7.2 Tukey-Anscombe plot

The Tukey-Anscombe P plot is a plot of the residuals against the fitted values ybi .
We always have ri ybi = 0, i.e. the Tukey-Anscombe plot always has sample
correlation zero. If this plot exhibits a non-linear structure, this is an indication
that the model assumptions are broken. If the residuals were plotted against
the yi , their correlation would make any interpretation more difficult.

For a simple linear regression, this is (essentially) equivalent to plotting the ri

against the xi (unlike multiple regression, for which the plot of ri against ybi is
more informative than the component-wise plots of ri against xij ). Likewise,
this plot is quite similar to the original scatterplot of yi against xi for a simple
linear regression. It is, however, easier to recognize deviations from a horizontal
line than deviations from a sloping line.

The ideal look of a Tukey-Anscombe plot is shown in Figure 1.17. One frequent
deviation from the assumption of constant error variance is a variance that
increases as the target variable does. The effect of this on the Tukey-Anscombe
plot can be seen in Figure 1.18 a)–c). If the Tukey-Anscombe plot contains
some structure in the form of a “trend”, this is an indication that the regression
function has not been specified correctly (i.e. the mean error is not zero.) Figure
1.18d is a typical example of a case in which a quadratic term is most likely
missing.

Figure 1.17: The Tukey-Anscombe Plot for an example in which the model
assumptions hold true.

If the Tukey-Anscombe plot shows some sort of connection between the error
variance and ybi or the variables in x, the target variables should be transformed
or a “weighted regression” performed (see Section 1.7.5). If the spread of the
errors increases linearly with the fitted values, a logarithmic transformation will
stabilize the variance; and if the error spread is proportional to the square root
of the fitted values, taking the square root of the target variable stabilizes the
variance. (This can be shown using Taylor expansions.)
40 Linear Regression

Figure 1.18: a) Linear increase in standard deviation, b) non-linear increase

in standard deviation, c) 2 groups with differing variance, d) missing quadratic
term.

1.7.3 Time series plot, Durbin-Watson test

If the errors are dependent, the levels of the tests and confidence intervals are
no longer correct. This can be seen quite easily: If ε ∼ N n (0, Σ), a simple
calculation shows that
b ∼ N p (θ, (X T X)−1 (X T ΣX)(X T X)−1 ).
θ

The size of the correlation effects between the errors also depends on the co-
variates X and on the shape of the covariance matrix Σ. In many cases, this
effect is considerable.
We cannot get far without making some assumptions about the dependence
structure: even if we knew the errors εi , we cannot estimate the n(n + 1)/2
entries of the covariance matrix from n data without some extra information.
If the observations form a time series, their covariance often is a (usually mono-
tonely decreasing) function of the time between observations. Dependence of
this type can be detected by plotting the residuals ri against the observation
times ti . Where these are unknown, the position k(i) of the observation in the
series may be used instead.
If the points vary randomly around the horizontal axis in the time series plot,
everything is fine. However, if adjacent ri are similar, this indicates that the
errors may be serially correlated. Sometimes we even observe a jump in the
level of the residuals. In such a case, the model has evidently changed suddenly
at a particular point in time.
It is possible to test independence against an alternative of serial correlation.
Two such tests are:
(i) The run test, which counts the number of continuous sub-sequences
(“runs”) in which the residuals have identical signs. When independence
is assumed, there should not be too many or too few runs.
(ii) The Durbin-Watson test, which uses the test statistic
Pn−1
(ri+1 − ri )2
T = i=1Pn 2 .
i=1 ri
1.7 Residual analysis, verification of models 41

Expanding this, we find that

Pn−1 !
r r
i=1 i i+1
T ≈2 1− P n 2 .
i=1 ri

(The deviation stemming from this approximation is marginal.) The quo-

tient in this formula is an estimate of the correlation of εi and εi+1 (as-
suming that all the εi have the same variance). If the εi are independent,
T is approximately 2; lower values of T indicate positive dependence.
Finding the critical values for this test is made all the more difficult by
the dependence of the distribution of ri , and thus the distribution of T , on
the experimental design, i.e. on the chosen xi . The test only looks at the
extremes over all designs, which leads to two different tabulated values
(see e.g. Sen und Srivastava (1990), p. 326). If T is smaller than the
lower tabulated value, the null hypothesis of “independence” is rejected.
If T is larger than the upper tabulated value, the null hypothesis is kept,
and between these two critical values the situation depends on the xi , i.e.
the test “abstains”.
The drawback of the Durbin-Watson test is its exclusive focus on the
correlation between observations that immediately follow each other.
Recall that the consequences of serial correlation also depend on the experi-
mental conditions (the values of the covariates). If we can freely choose the
order of our observations, then – in a univariate case – it would be tempting to
make our observations in ascending order of x, as this often is the easiest to do.
In such an experimental setup, however, the effect of positive error correlation
on the slope estimate is particularly large – and thus it is not a good choice. In
other words, error trends in time are mixed with the effect of the covariate.
It is much better to choose values of x that are as “orthogonal to time” as
possible, i.e. such that the xi and the ti (or k(i)) are uncorrelated. Then error
trends will cancel each other out during slope estimation (linear trends will
cancel out exactly, non-linear ones approximately) and positive error correlation
will thus only have a small effect on the variance of the slope estimate. The
same is not true for the intercept, though: there the variance of the arithmetic
mean is large when the observations have high correlation. The easiest way to
obtain near-orthogonality of covariates to time is by randomizing the order of
the values of x.

1.7.4 Interior analysis

“Interior analysis” is the local estimation of error variance using replicates or

“near replicates”. “Replicates” are repeated measurements at the same xi , and
“near replicates” are measurements made close to xi (cf. Figure 1.19). “Interior
analysis” checks for a “lack of fit”, i.e. a systematic error in the class of fitted
models. In other words, it checks the assumption E [εi ] = 0.
42 Linear Regression

y6 
c
c y32 c
s
c y34
c y3. sc y
33
c c

c c y31 s Replicates
s 
c c
c c Data
c  point yij
c Near replicates s Mean yi.
-
x2 x1 x6 x3 x4 x5 x
Figure 1.19: Replicates

The model with ni replicates at xi can be formulated thus:

Yij = xTi θ + εij(i = 1, . . . , k; j = 1, . . . , ni ),

where εij i.i.d. ∼ N (0, σ 2 ). The sample size is n = ki=1 ni . We can compare
P
this model to the larger model where Y and X are connected by an arbitrary
function f . Introducing the quantities E [Yij ] = f (xi ) = µi as independent
parameters, we have:
Yij = µi + εij (i = 1, . . . , k; j = 1, . . . , ni ).
The least squares estimate of µi is simply the (arithmetic) mean of the obser-
vations at xi :
Xni
µ
bi = yi. = yij /ni .
j=1
Thus we have two nested models and can perform the usual F -test. The or-
thogonal decomposition for the ANOVA table is:
X ni
k X X ni
k X k
X
2 2
(yij − ybi ) = (yij − yi· ) + ni (yi· − ybi )2
i=1 j=1 i=1 j=1
| {z } |i=1 {z }
(1) (2)

• (1) measures the random error, and has ki=1 (ni − 1) = n − k degrees of
P
freedom
• (2) measures the random error and the “lack of fit”, and has k − p degrees
of freedom
If we have two near replicates, we can correct them by moving each of them
parallel to the regression line so that they both have the average of their co-
variates x (cf. Figure 1.20). This then gives us two replicates. When p is large,
however, it can be hard to find near replicates.
1.7 Residual analysis, verification of models 43


 
 
 
g
 
  w9 
 


: w
g 


Figure 1.20: Correcting near replicates

1.7.5 Generalized least squares, weighted regression

Model: Y = Xθ + ε as previously, but now with ε ∼ N (0, σ 2 Σ)

.
We assume Σ to be known and σ 2 unknown – i.e. the error covariance matrix
is given up to a multiplicative constant. We furthermore assume that Σ is
positive definite. Then there exists a regular matrix A for which AAT = Σ,
i.e. A is a square root of Σ (cf. Appendix A).
Reduction to the standard model: We compute the transform

ỹ := A−1 y = A−1 (Xθ + ε) = A −1

| {zX} θ + A
−1
| {z ε} = X̃θ + ε̃
X̃ ε̃

Then we have

E [ε̃] = E A−1 ε = A−1 E [ε] = 0

 

Cov[ε̃] = Cov[A−1 ε] = A−1 Cov[ε](A−1 )T

= A−1 σ 2 (AAT )(A−1 )T = σ 2 I.

In other words, the “tilde model” obtained by the linear transformation using
A−1 satisfies the conditions for the standard multiple regression model we al-
ready know. The key point here is that A is invertible (which is ensured by the
positive definiteness of Σ).
Applying the known theory to the “tilde model”:
For the “tilde model”, we estimate θ using least squares, that is, we minimize

kỹ − X̃θk2 = (y − Xθ)T A−T A−1 (y − Xθ) = (y − Xθ)T Σ−1 (y − Xθ).

This is the same as performing a least squares estimation for the original data
(y, X) using a different scalar product. The estimate obtained by this is
b = (X̃ T X̃)−1 X̃ T ỹ = (X T Σ−1 X)−1 X T Σ−1 y,
θ

which is the so-called generalized least squares estimate of θ. Its distribu-

tion is
b ∼ N p (θ, σ 2 (X T Σ−1 X)−1 ).
θ
44 Linear Regression

Tests, confidence intervals, etc. can be constructed in a similar way to the

standard case. The results of Section 1.9 will show that if Σ 6= I, the generalized
least squares method has a smaller variance than the standard least squares
method.
One important special case occurs when Σ is a diagonal matrix, i.e. when the
errors are uncorrelated, but have varying degrees of precision:
 
v1 0
 v2 
Σ= (vi > 0 ∀i)
 
.. 
 . 
0 vn

In such a case we introduce weights wi which are proportional to v1i , i.e. we

minimize i wi ri2 . The more precise an observation, the greater a weight is
P
given to it by the generalized least squares procedure.
It is, however, a fairly rare occurrence that we would actually know the error
covariance matrix up to a constant factor. If we do not have such informa-
tion, we often first use standard least squares, estimate a covariance matrix Σb
from the residuals, and then use this estimated covariance matrix to perform
generalized least squares. This procedure is frequently used when the errors
exhibit correlation in time, and it is usually referred to as the Cochrane-Orcutt
procedure.

1.8 Model selection

We assume that our observations have been generated by the model

yi = f (xi ) + εi (i = 1, . . . , n),
where εi are i.i.d. with E [εi ] = 0, Var(εi ) = σ 2 . Moreover, we model the
regression function as
p
X
f (xi ) ≈ θj xij .
j=1

Here we always assume an intercept to be present, i.e. xi1 = 1 for all i.

Question: Which variables should we include in the model?
Naive answer: “More is always better!”.
At a general level, this answer is wrong, as it is quite possible that some of
these p variables are “superfluous” or make only a very marginal contribution
to the explanatory power of the model. At the same time, estimating more
coefficients raises the random error of our parameter estimates and predictions.
Thus this answer is not correct.
Our choice of variables may be determined by theoretical considerations in
the our field of work (such as in physics). This is not model selection in a
1.8 Model selection 45

statistical sense. Furthermore, we do not consider transformations of y and x,

unlike before. The only question is whether or not a variable (that may have
been transformed earlier on) is included or omitted.
Our search for the “best model” depends on the question we are trying to
answer:
(i) Using regression to find an explanatory model
(ii) Using regression to make predictions
We shall first discuss stepwise methods, and then we shall look at methods which
evaluate all of the 2p−1 possible models according to some suitable criterion and
choose the “best” one. Stepwise methods of course take less effort.

1.8.1 Model selection using “stepwise regression”

(1) Forward stepwise regression:

Start with the model containing only the intercept, i.e. with only the
constant xi1 :
yi = θ1 + εi .
(Of course the estimate of θ1 is then merely the average of all observa-
tions). Now the variables are included one by one (and “stepwise”); at
each step the model with the most significant F –value over the previous
model is included.
 y  
  ee
  
 
e
 e
 
r
 
e
 e
ŷ(q+1)
Q 
Q 
Q
Q 
Q 
Q
Q 
Q 
Q
Q 
Q 
Q 
Q r
Q  ŷ
Q (q)
Q
Q

Figure 1.21: Previous model (q-dimensional) and new enclosing model ((q +1)-
dimensional, i.e. with one more variable).

Stopping condition: Repeat until none of the F –statistics are signifi-

cant (for a given significance level). As we are carrying out repeated tests,
we should take care not to misinterpret this level.
(2) Backwards stepwise regression:
46 Linear Regression

Start with the full model:

yi = θ1 + θ2 xi2 + . . . + θp xip + εi

Now the variables are taken out of the model one by one (and “stepwise”).
At each step, we remove the variable whose F –value in the comparative
test is smallest, and continue until all such F are significant.
Discussion of stepwise regression:
• “Backward” selection requires more effort, and under some circumstances
it may even incur more numerical problems (if e.g. p ≥ n, we cannot even
carry out this procedure). However, it gives us the greatest certainty of
finding a good model.
• “Forward” and “backward” selection are often combined (using two dif-
ferent significance levels, so as to avoid entering a loop of removing and
including the same x(j) ).
• The stopping rule does not necessarily gives us a “best” model according
to the criteria we shall discuss in the next section.
• The sequence of included or removed variables should not be regarded as
some kind of rank of the importance of the variables.
• “Forward” and “backward” selection may yield entirely different solutions.
Example of this last point: We choose three explanatory variables such that
• X1 and X2 each have only very little correlation with Y , but Y is (almost)
a linear combination of X1 and X2 .
• X3 correlates strongly with Y .
Forward selection will first select X3 and then stop (or perhaps choose {X1 , X3 }
or {X2 , X3 }), while backward selection will choose {X1 , X2 } and then stop.

1.8.2 Model selection criteria

Mallows’ Cp statistic

The Cp statistic is an estimate of the mean quadratic prediction error of a fitted

model which is averaged over the observed experimental conditions xi (i =
1, . . . , n). It also takes into account the bias of a badly-fitting model, but it
requires a good (bias-free and sufficiently exact) estimate of σ 2 (e.g. one taken
from a “full” model that might contain far too many variables for its fit to be
used, one estimated using replicates or “near replicates”, or one known from
experience). It automatically “penalizes” superfluous variables and can be used
as a quality measure (estimated) of a model.
By the assumptions made at the beginning of this section, the variables yi are
independent with mean E [yi ] = f (xi ) = µi and variance Var(yi ) = σ 2 . Each
model is described by the subset M ⊂ {1, 2, . . . , p} of variables included in it.
Here we assume that xi,1 ≡ 1 for each index i, and that each M contains the
1.8 Model selection 47

index 1. We denote the corresponding matrix of covariates by X M (a submatrix

of X), i.e.
X M = (xij ; 1 ≤ i ≤ n, j ∈ M ).
Next we use least squares to estimate the parameter of the model M by

bM = ((X M )T X M )−1 (X M )T y
θ

and the mean vector µi = E [yi ] by

M
bM = X M θ
y b .

So we a fit a linear model M with |M | variables – which may or may not be

correct. For instance, the fit may have the shape of a straight line, even though
the expectations E [yi ] lie on a parabola; cf. Figure 1.22.

Figure 1.22: Theoretical best linear fit (dashed line) for a model that is actually
quadratic (solid line).

If the model is correct, we have

= X M ((X M )T X M )−1 (X M )T µ,
 M
E y b

which amounts to the best approximation of µ by the variables in the model. Of

course, this approximation gets better and better as we include more variables
in our model (it usually even improves quite clearly, as in practice the influence
of a variable is rarely zero, even if it is small). Furthermore, we have

yM ) = σ 2 X M ((X M )T X M )−1 (X M )T ,
Cov(b

i.e. the random fluctuations are as if the model were true. In particular, we
have
Xn
yiM ) = σ 2 tr(X M ((X M )T X M )−1 (X M )T ) = |M |σ 2 .
Var(b
i=1
The more variables we include in our model, the greater the sum of variances
of ybiM .
48 Linear Regression

Our quality measure for a model is the sum of mean squared errors ybiM from
the true values µi :
n
 M  X
y − µ||2 =
 M
yi − µi )2 .

SM SE = SM SE(M ) = E ||b E (b
i=1

Now for any random variable Z and each constant c , we know that

E (Z − c)2 = E ((Z − E [Z]) + (E [Z] − c))2 = Var(Z) + (E [Z] − c)2 + 2 · 0,

   

and thus
n
X n
X n
X
yiM ) ybiM 2 2
(E ybiM − µi )2 .
   
SM SE = Var(b + (E − µi ) = |M |σ +
i=1 i=1 i=1

The first of these summands is small for models containing few variables, while
the second is small for models containing many variables. SM SE is often scaled
by σ 2 , so that
SM SE(M )
Γp (M ) =
σ2
becomes the term of interest. The inequality Γp (M ) ≥ |M | always holds, with
equality being attained exactly when the model M is bias-free (but potentially
contains superfluous terms).
We can also regard ybiM as a prediction for a new observation Yn+i = µi + εn+i .
In this case the sum of prediction square errors is
n
X n n
 X  X
E (Yn+i − ybiM )2 = E (Yn+i − µi )2 +
 M
yi − µi )2 = nσ 2 +SM SE.
  
SP SE = E (b
i=1 i=1 i=1

(To be more precise, we should actually refer to this as “sum of mean squared
prediction errors”.)
Thus minimizing SM SE or Γp or SP SE always leads to the same model.
However, we cannot compute any of these quantities without knowledge of σ
and µ. One naive estimate of SP SE is the sum of squared errors
n
X
b M ||2 =
SSE(M ) = ||y − y (yi − ybiM )2 ,
i=1

i.e. the criterion used in least squares estimation. However, this quantity
always decreases when additional variables are included in the model, and it
underestimates SP SE:
n
X n
X
b M ||2 = Var(yi − ybiM ) + (E [yi ] − E ybiM )2
   
E ||y − y
i=1 i=1
n
X
= (n − |M |)σ 2 + (E [yi ] − µi )2 = SP SE(M ) − 2|M |σ 2 .
i=1
1.8 Model selection 49

It is therefore better to estimate SP SE by SSE(M )+2|M |b σ 2 , where σ

b2 denotes
2
an estimate of σ (e.g. from the full model M = {1, 2, . . . p}), and to then choose
the model minimizing this estimate of SP SE. Similarly, we can use

SSE(M )
Cp (M ) := − n + 2|M |,
b2
σ
as an estimate of Γp . Random fluctuations mean that Cp can become smaller
than |M | (or even negative), unlike Γp .

Figure 1.23: Cp plot: a) best prediction estimate, b) best bias-free prediction

estimate, c) Point at which the random fluctuation of Cp causes it to be below
the model size.

Akaike’s information criterion AIC

For an arbitrary model (not necessarily a regression model) with k parameters,

we define
AIC(α) = −2`bk + α k
where `bk denotes the maximum log–likelihood in the model (i.e. the log-
likelihood at the MLE). Thus the first term measures the goodness of fit, and
the second penalizes the complexity of the model. The multiplicative constant
in the penalty term is usually taken to be 2. When selecting between different
candidate models, we can now choose the one which minimizes the AIC.
In particular, if we have a linear model with normally distributed errors and a
selection M ⊂ {1, 2, . . . p} of explanatory variables, the log-likelihood is
n 1
log fM (y, θ) = − (log(2π) + log(σ 2 )) − 2 (y − X M θ M )T (y − X M θ M ).
2 2σ
Assuming σ 2 to be given (e.g. by experience or by prior investigation), we then
see that
1 bM )T (y − X M θ
bM ),
−2`bM = 2 (y − X M θ
σ | {z }
SSE(M )
50 Linear Regression

(up to a constant), i.e. that AIC(2) is equal to Cp up to a constant (whose

exact value is irrelevant for model selection) when both use the same σ.
If furthermore σ is replaced by its maximum likelihood estimate
SSE(M )
b2 (M ) =
σ ,
n
Akaike’s criterion is

σ 2 (M )) + α |M |
AIC(α) = n log(b

(again up to a constant). Using the Taylor expansion of the logarithm at σ 2 ,

we conclude that
SSE(M )
AIC(2) ≈ n log(σ 2 ) + − n + 2|M |.
σ2
Thus we see that even in this more general case the AIC is very similar to Cp
– at least for those models, in which SSE(M )/n lies near the estimate of σ 2
used to compute Cp .

1.9 The Gauss-Markov theorem

This theorem states the “optimality” of the least-squares estimate in a certain

sense. There is a version under the assumption of normality, and one without
this assumption. The difference between these two is not only in the assump-
tions; there are also significant differences between the statements they make!
First we give the result that does not assume normality.
Theorem 1.9.1 (Gauss–Markov). Let

Y = Xθ + ε E [ε] = 0 Cov[ε] = σ 2 I rank[X] = p.

Furthermore, let c be an arbitrary p-dimensional vector, and θ

b the least squares
T
estimator. Then c θ has minimal variance amongst all linear unbiased estima-
b
tors of cT θ.
Because of this, we also call least squares estimators “BLUE” (“best linear
unbiased estimators”).
The version assuming normality is as follows:
Theorem 1.9.2. Let furthermore ε be normally distributed. Then cT θ
b has
T
minimal variance amongst all unbiased estimators of c θ.
Thus we say that least squares estimators are “UMVU”(“uniformly minimum
variance unbiased”). Uniformly, as this holds for arbitrary values of θ and σ 2 .)
It is this theorem (not the Gauss–Markov theorem) that best justifies using
least squares (besides the simplicity of this method). However, unbiasedness is
not always absolutely necessary, and does not hold e.g. in Bayesian regression
and in “ridge regression”. We shall not discuss these methods here, though.
1.9 The Gauss-Markov theorem 51

The omission of the normality assumption in the Gauss–Markov theorem is

paid for by the restriction of this result to linear estimators. Crucially, though,
all linear estimators are less efficient (i.e. have a higher variance), even if
the deviation from normality is relatively small.

Illustrative examples: If we attempt to approximate the true error distri-

butions of high-quality data by a t distribution with ν degrees of freedom, we
find ν = 5 − 9 degrees of freedom for data appearing to be “quite normal”,
but often only ν = 3 degrees of freedom – and even ν = 1 (Cauchy distribu-
tion) can occasionally occur. Now the asymptotic efficiency (which is e.g. the
inverse ratio of the required sample sizes needed to obtain the same precision)
of the least squares estimator compared to an asymptotically best estimator
(such as maximum likelihood) under a t distribution with ν degrees of freedom
is = 1 − 6/(ν(ν + 1)) (for ν ≥ 2). Thus the actual efficiency of the least squares
method is 80-93% for “quite normal” data (t5 − t9 ), and ca. 50% for “fairly
b2 , the situation is much worse still !
normal” data (t3 ). If we look at σ

Proof of the Gauss–Markov theorem Let a be an n–dimensional vector

and a0 a constant for which aT y + a0 constitutes an unbiased estimate of cT θ.
Then we have
E aT y + a0 = aT Xθ + a0 = cT θ
 

for all θ. From this we conclude that a0 = 0 and aT X = cT ⇔ X T a = c.

The vector aKQ = X(X T X)−1 c from the least-squares estimate is a spe-
cial solution of X T a = c. Furthermore, aKQ is orthogonal to all solutions
ah of the homogeneous system of linear equations, X T a = 0, as aKQ T ah =
cT (X T X)−1 X T ah = 0. As Cov[Y] = σ 2 I, we thus can see that

Var((aKQ + ah )T Y) = Var(aKQ T Y) + Var(ah T Y) ≥ Var(aKQ T Y).

Proof of the variant Gauss–Markov theorem (using the multidimensional

Cramér-Rao inequality, which is proven in the class on Mathematical Statistics):

We regard the following general situation: Let (fη (y)) be a parametric family
of strictly positive densities in Rn . Let η be a variable parameter with values
in an open subset of Rk , and let fη (y) be differentiable wrt η. Our parameter
of interest is g(η), where g is an arbitrary real-valued function of η. Then we
have
Theorem 1.9.3 (Cramér-Rao). If T (y) is an arbitrary unbiased estimate of
g(η), i.e.
Eη [T (y)] = g(η) ∀η,

then g is differentiable and

∂g T ∂g
Varη (T (y)) ≥ I(η)−1 ,
∂η ∂η
52 Linear Regression

where I(η) denotes the so-called Fisher information matrix:

" #
∂ log fη (Y) ∂ log fη (Y) T
I(η) = Eη .
∂η ∂η

We now apply this to

1
fη (y) = (2π)−n/2 (σ 2 )−n/2 exp(− (y − Xθ)T (y − Xθ)),
2σ 2

η = (σ 2 , θ T )T
and
g(η) = cT θ.

After some calculations, we obtain the Fisher information

 n 
2σ 4 0
I(η) = 1 ,
0 σ2
XT X

i.e. the least squares estimate cT θ

b attains the Cramér-Rao lower bound. It
thus obviously has minimal variance. 2
Chapter 2

Nonlinear and nonparametric

methods

2.1 Robust methods

We generally call a statistical procedure for a parametric model robust if its

properties do not change much when slight deviations from the model occur.
For the linear model, the method of least squares is not robust, as heavy-tailed
error distributions have serious ramifications ! The distribution of the estimates
is fairly stable (cf. Section 1.4.3), but for comparatively small deviations from
normality, there are better estimators than least squares. Consequently, the
level of tests and confidence intervals based on least squares is robust, but their
power is not.
One related effect is that least squares estimates and the subsequent tests and
confidence intervals are extremely sensitive to individual outliers – and obser-
vations that look like outliers are produced quite often when the underlying
error distribution is heavy-tailed.
We will now first take a closer look at the effects of an outlier on the least
squares method; then we shall discuss more robust alternatives.

2.1.1 Influence of individual observations on the LSE

First we examine the effect that omitting or adding observations has on the
least squares estimate. The following lemma of Gauss will be useful here:
Lemma 2.1.1. Let A be an invertible matrix of dimension p × p, and let a and
b be two vectors of dimension p for which bT A−1 a 6= −1. Then A + abT is
also invertible, and we find that

1
(A − abT )−1 = A−1 + A−1 abT A−1 .
1− bT A−1 a

53
54 Nonlinear and nonparametric methods

Proof: exercise.

We denote by θb(−i) the least squares estimate when the i-th observation is
omitted, and we furthermore use the shorthand notation
n
X n
X
T
A=X X= xj xTj , c=X y= T
yj xTj .
j=1 j=1

Then the above lemma tells us that

b(−i) = (A − xi xT )−1 (c − yi xi )
θ i
1
= A−1 c − yi A−1 xi + A−1 xi xTi A−1 (c − yi xi )
1 − xTi A−1 xi
b − yi A−1 xi (1 + xTi A−1 xi 1
= θ )+ xTi θA−1 xi ,
1 − xTi A−1 xi 1− xTi A−1 xi

and thus

b(−i) − θ 1
θ b = − (X T X)−1 xi (yi − xTi θ)
b
1 − xTi (X T X)−1 xi
ri
= − (X T X)−1 xi .
1 − Pii

We see that the influence of the i-th observation depends on the i-th residual and
also on the diagonal entry Pii in the “hat matrix” (when using all observations,
including the i-th). Thus a plot of the residuals ri against the Pii is often used
to detect influential observations.
The difference of the parameter estimates is somewhat difficult to interpret, as
it is a whole vector that moreover depends on how the covariates are scaled.
One scalar invariant can be obtained by computing the length of θb(−i) − θ
b with
respect to the metric defined by the estimated covariance matrix of θ:
b

(−i) (−i)
(θ
b b T (X T X)(θ
− θ) b − θ)
b 1 ri2 Pii
Di = 2
= 2
.
pb
σ pσb (1 − Pii ) 1 − Pii

We refer to Di as Cook’s distance.

√ It is a simple function of Pii and the square of
σ 1 − Pii ), and is frequently used as a diagnostic
the studentized residual ri /(b
tool. Observations yielding a much larger value of Di than the others should
be looked at more closely or omitted altogether.
We can likewise regard the change that occurs when an observation is added at
an arbitrary location (y, x):

1
∆θ
b= (X T X)−1 x(y − xT θ).
b
1+ xT (X T X)−1 x

We see that the least squares estimate can be changed arbitrarily much by a sin-
gle new observation – in other words, the least squares estimator is not robust.
2.1 Robust methods 55

Moreover, this effect depends quite strongly on the location of the new obser-
vation. This formula can be made a little clearer if we assume that the xi are
chosen randomly and are i.i.d. . We then obtain the first-order approximation
1 
(E xi xTi )−1 x(y − xT θ)

b∼
∆θ
n
for n → ∞.
The discovery – and subsequent special treatment – of influential observations
using Cook’s distance does have two drawbacks, however: First, the effect of
omitting two or more observations is not merely the sum of the individual effects
(as one influential observation can mask others). Secondly, omitting influential
observations leads to question marks over the validity of tests and confidence
intervals based on the remaining data.

2.1.2 Huber and L1 regression

The reason for the large influence that individual observations can have on
the least squares estimate is that large residuals have a high weight when a
quadratic criterion is being used. To avoid this, we can instead look at the L1
estimator:
Xn
θ = arg min
b |yi − xTi θ|.
θ
i=1

Historically speaking, this method is even older than least squares: it was
suggested and discussed by Boscovich in 1760 and by Laplace in 1789 !
In the location model, i.e. when p = 1 and xi ≡ 1, the solution is the median of
the data, an estimator which for normally distributed data is substantially less
precise than the arithmetic mean (i.e. the least squares estimator): to reach
the same precision, the median requires 50% more observations.
One compromise between minimizing the L2 distance and minimizing the L1
distance is given by Huber regression:
n
X
θ
b = arg min ρc (yi − xTi θ),
θ
i=1

where
1 c
ρc (u) = u2 (|u| ≤ c), ρc (u) = c(|u| − ) (|u| ≥ c),
2 2
cf. Figure 2.1. Choosing c = 0 leads to L1 regression. If we compute derivatives
and set them to zero, we obtain the equations
n
X
ψc (yi − xTi θ)xi = 0,
i=1

where ψc (u) = ρc (u)0 = sign(u) min(|u|, c).

56 Nonlinear and nonparametric methods

Figure 2.1: The Huber function and its derivative for various values of c

However, Huber regression only makes sense if the “corner” at c is chosen in

relation to the variance of the residuals. Thus the estimators generally used
are:
n
!
X yi − xTi θ
b
ψc xi = 0,
σ
b
i=1
n
!
X yi − xTi θ
b
χ = 0.
σ
b
i=1

The function χ(u) is chosen either to be ψc (u)2 −β or to be χ(u) = sign(|u|−β).

The constant β is fixed by the condition
Z
χ(u) exp(−u2 /2)du = 0,

which ensures that for normally distributed errors σ b is a consistent estimate

of the standard deviation. The first choice, χ, is Huber’s so-called Proposal 2,
while the second choice is 1/β times the median of the absolute residuals.
Closed-form computation of the L1 and Huber estimators is no longer possible,
but efficient algorithms for them are now known. The computational problem
of L1 regression can even be reduced to a linear optimization problem that “in-
terior point” methods solve more quickly than least squares. Huber regression
is performed by iterating the weighted least squares procedure using the weights
ψc ((yi − xTi θ)/b
b σ) cb
σ
wi ∝ ∝ min(1, )
T
yi − x θ b |yi − xT θ|
b
i i
until they stabilize.
Furthermore, no closed-form expression exists for the distribution of the L1 or
Huber regression estimators. Asymptotic methods are thus used to show that
when the covariates xi are random, independent and identically distributed, the
√ b
standardized vector n(θ − θ) approximately has a normal distribution with
expectation zero and covariance matrix
E ψc (εi /σ)2 2 
 
−1
2
σ E xi xTi .
P [|εi | ≤ cσ]
2.1 Robust methods 57

Up to a constant factor, this covariance matrix is the same for least squares.
When c lies in the interval [1, 1.5], this factor is much less than 1 for heavy-
tailed distributions and only slightly larger than 1 for the normal distribution.
When comparing numerical estimates, however, we should remember that for
non-normal errors εi , the parameter σ is no longer the error standard deviation,
but the solution of
E [χ(εi /σ)] = 0.

This asymptotic approximation is also the basis of tests and confidence intervals,
whose details we omit here.
Unfortunately Huber regression does not solve all the problems we have with
influential observations. The exact effect of adding or omitting an observation
can no longer be specified. We can, however, approximate the difference in the
estimator caused by the addition of an observation at (x, y) by

1 y − xT θ
(E xi xTi )−1 x ψc (
 
b∼
∆θ )σ.
nP [|εi | ≤ cσ] σ

Thus the influence of large values of y is bounded when x is fixed, but by varying
x we can nonetheless increase this influence arbitrarily. More refined methods
lacking this weakness will be discussed in the next two sections.
Huber regression is frequently replaced by estimators that assume ψ to be
an odd function and χ an even function. Such estimators are known as M -
estimators. Particularly popular choices of ψ include those that converge to
zero for |r| large, as they then remove large outliers entirely. However, this
usually comes at the price of non-unique solutions to the defining equations.
The solution actually found then depends on the algorithm – and especially on
the choice of starting value.

2.1.3 Regression estimators with restrictions on influence

To restrict the influence of both y and x, we regard estimators defined by the

following type of equations:
n
!
X yi − xTi θ
b
η xi , xi = 0.
σ
b
i=1

To compute σ b, we use an equation similar to that in Huber regression. Several

kinds of η(x, r) are used here, including
a
η(x, r) = min(1, )ψc (r) (Mallows)
||Ax||

and
1
η(x, r) = ψc ((||Ax||)r) (Schweppe).
||Ax||
58 Nonlinear and nonparametric methods

The matrix A is to be chosen in such a way that ||Ax|| expresses the deviation
of x from the set of covariates (xi )1≤i≤n . This can be achieved e.g. by setting

||Ax||2 = const · xT (X T X)−1 x,

though this choice in turn can be influenced quite strongly

P byT an observation
T
that has an unusual xi , since X X can be rewritten as i xi xi . Thus X T X is
replaced by a similar, yet robust, quantity. This leads to additional equations
for A which we shall not discuss any further here.
Mallows’ choice of η always chooses lower weights for observations with strongly
deviating explanatory variables. To better understand Schweppe’s choice, it
helps to look at the identity ψc (dr)/d = ψc/d (r). Thus Schweppe’s suggestion
merely lowers the corner in the Huber function, thus enabling an observation
with strongly deviant covariates to nonetheless have full weight, if the corre-
sponding residual is close to zero.
For this procedure, the difference in estimates after adding an observation at
(x, y) is approximately given by
−1
y − xT θ
   
1 ∂ εi
b∼
∆θ E η(xi , )xi xTi x η x, σ.
n ∂r σ σ

Thus the influence of both x and y is restricted under both choices of η.

However, these estimators are still not satisfactory, as their “breaking point” is
no greater than 1/p. The breaking point is defined as the maximum proportion
of outliers an estimator can withstand without diverging.

2.1.4 Regression estimators with high breaking point

We can obtain estimators whose breaking point does not depend on the dimen-
sion by replacing the arithmetic mean by the median in
n n
X 1X
arg min (yi − xTi θ)2 = arg min (yi − xTi θ)2
θ θ n
i=1 i=1

to obtain
b = arg min median((yi − xT θ)2 )
θ i
θ

(Least median of squares; Hampel 1975, Rousseeuw 1984). In other words:

amongst pairs of parallel hyperplanes sandwiching 50% of all observations (yi , xi ),
we look for the pair whose distance along the y-axis is minimal. This procedure
is illustrated in the left half of Figure 2.2.
It is intuitively clear (and provable) that this procedure can tolerate outliers in
roughly 50% of all observations without diverging. Computing this estimator
is a much greater problem, however, as the target function median((yi − xTi θ)2 )
generally has many local minima (cf. Figure 2.2, right side). Thus we need to
2.2 Nonlinear least squares 59

Figure 2.2: Least median of squares for simple regression. On the left, two
lines with a band containing 50% of the observations. The intercept is chosen so
as to minimize the band’s diameter along the y-axis. On the right, the diameter
of the band as a function of its slope.

search the entire space to find the global minimum, and this quickly becomes
too costly as we move to higher dimensions. In general, stochastic algorithms
are used to select p + 1 data points at random, fit a plane through them and
then compute the value of the target function for the corresponding θ.
One further drawback of this method is its lack of efficiency in the normal case:
the estimator then only converges at the rate n−1/3 . A better convergence
rate can be obtained e.g. by replacing the median by a truncated mean of the
((yi − xTi θ)2 ), omitting the αn (for α < 0.5) largest squared residuals. Usually
this estimator is used as the starting point from which to perform a Newton
iteration to solve the estimator equations for an M -estimator using a function
ψ that converges to zero. The result is then called the M M -estimator.
Developing robust regression estimators that exhibit good statistical and algo-
rithmic properties is still an on-going topic for research.

2.2 Nonlinear least squares

In this chapter, we discuss methods for estimating θ in models of the form

yi = f (xi , θ) + εi .

Here f is a known function of the experimental conditions and of the parameters.

The key assumptions on f that we shall make are that it is nonlinear in the
parameters θ, and that we cannot (or do not want to) transform it into a
linear model. Now p, the dimension of the parameter, need no longer be the
dimension of the explanatory variables. For the vector ε of errors, we make the
same assumptions as in the linear model, i.e. E [ε] = 0 and Cov(ε) = σ 2 I.
60 Nonlinear and nonparametric methods

Many applied problems are of this type, and the shape of f usually follows
from the theory of the science in which a particular application arises. As an
example, the description of the cumulative oxygen usage y of microorganisms
in samples of river water as a function of incubation time x is usually performed
by the model
f (x, θ) = θ1 (1 − exp(−θ2 x)).
Thus the parameter θ1 is the saturation point, and θ1 · θ2 is the slope at x = 0.
Some sample data and a possible regression function are shown in Figure 2.3.

Figure 2.3: Oxygen usage data as a function of incubation time (left), and a
typical regression function (right).

Another, similar, example is given by the so-called Michaelis-Menten model,

which describes the dependence of reaction speed y on substrate concentration
x. In this model, we have
θ1 x
f (x, θ) = .
θ2 + x
The transformation y → 1/y, x → 1/x turns this into a linear regression model.
However, the data shown in Figure 2.4 no longer have constant variance once
this transformation has been applied; thus nonlinear least squares provide a
much better fit.
In many applications the quantities xi are times or locations at which a vari-
able is observed whose development satisfies an ordinary or partial differential
equation. In such a case, the parameters θ are the parameters of the differential
equation (including boundary conditions, if required), and f (x, θ) denotes the
solution of the differential equation for parameter θ at x.
The least squares estimate is defined to be

θ
b = arg min S(θ),
θ

where
n
X
S(θ) = (yi − f (xi , θ))2 .
i=1
This carries the following geometric interpretation: If we vary θ, the points
(f (x1 , θ), . . . , f (xn , θ))T describe a p-dimensional curved surface in Rn , the so-
2.2 Nonlinear least squares 61

Figure 2.4: Reaction time as a function of substrate. On the left, a fitted

line after a linearizing transformation; on the right, fitted functions on the
untransformed scale estimated in the untransformed variables (solid line) and
the transformed variables (–·–), respectively.

called response surface. We now seek the point on the response surface that is
closest to the observation (y1 , . . . , yn )T .
The solution does not generally have a closed form, and thus iterative meth-
ods (Gauss-Newton, Levenberg-Marquardt) are used. In practice, having good
starting values is crucial.
The error variance σ 2 can be estimated in a way similar to that for linear
regression:
n
1 X
b2 =
σ b 2.
(yi − f (xi , θ))
n−p
i=1

2.2.1 Asymptotic confidence intervals and tests

In non-linear models, even the assumption of normal errors is not enough to

allow exact tests and confidence intervals. We must thus rely on asymptotics.
The asymptotic properties of θ
b can be obtained by means of a Taylor approx-
imation around the true parameter θ 0 . Thus we approximate the response
surface in ia neighbourhood of the true parameter by a plane:

f (xi , θ) ≈ f (xi , θ 0 ) + a(θ 0 )Ti (θ − θ 0 ),

where  T
∂
a(θ)i = f (xi , θ); j = 1, . . . p .
∂θ j
This means that we approximately have a linear model with explanatory vari-
ables a(θ 0 )i near the true parameter. Furthermore, it can be shown that (under
certain technical conditions) the distribution of θ
b is asymptotically the same as
in the approximating model, i.e.

b asymptotically
θ ∼ N (θ 0 , σ 2 (A(θ 0 )T A(θ 0 ))−1 ).
62 Nonlinear and nonparametric methods

Here the matrix A(θ) consists of the rows a(θ)Ti .

These approximated distributions do not yet allow the construction of tests

and confidence intervals, as both σ and A(θ 0 ) are still unknown. However, we
can plug in σb and A(θ)
b instead. As in the linear model, we generally use the
t distribution instead of the normal distribution here, and the F distribution
instead of the chi-squared distribution. Thus the confidence interval for θk is
q
θbk ± tn−p;1−α/2 se(θbk ), se(θbk ) = σ
b ((A(θ) b −1 )kk .
b T A(θ))

2.2.2 More precise tests and confidence intervals

In a similar way to the F test used in linear regression, we can test two nested
models by means of the differences in sums of squared deviations. To test the
null hypothesis Bθ = b, we also need the least squares estimate under the null
hypothesis:
θ
b0 = arg min S(θ).
θ;Bθ=b

Then we compute the test statistic

b0 ) − S(θ))/q
(S(θ b
T = ,
S(θ)/(n
b − p)

where q denotes the rank of B. In the linear model, this statistic was identical
to the one computed from joint normal distribution of θ b and had an Fq,n−p
distribution under the null hypothesis. In the nonlinear setup, these two statis-
tics differ, and they are only approximately F distributed. This approximation,
however, is often significantly better the the normal approximation of θ.
b

In particular, we can use this to test the null hypothesis θk = θk∗ for any arbitrary
but fixed value of θk∗ . The test statistic then becomes

(−k) (−k)
S(θ
b ) − S(θ)
b S(θ
b ) − S(θ)
b
Tk (θk∗ ) = = ,
S(θ)/(n
b − p) b2
σ

where
b(−k) = θ
θ b(−k) (θ∗ ) = arg min S(θ).
k ∗ θ;θk =θk

(Thus θ b(−k) is the least squares estimate under the null hypothesis, i.e. it is
equal to θb0 in the above notation). Since an F distribution with one degree
of freedom in the denominator is simply the distribution of the square of a
t-distributed variable, we can also regard the test statistic

sign(θk∗ − θbk )
q
τk (θk∗ ) = b(−k) ) − S(θ).
S(θ b
σ
b
2.2 Nonlinear least squares 63

Under the null hypothesis, this follows a t distribution with n − p degrees of

freedom. By reversing the test, we obtain the following confidence interval for
θk :  q 
(−k)
θk∗ | S(θ
b ) − S(θ)
b ≤ tn−p;1−α/2 σ
b .

If we apply a monotone transformation to θk , this interval can simply be

transformed in the same way, an advantage not shared by the interval θbk ±
tn−p;1−α/2 se(θbk ). The differences between these two intervals show us the ef-
fects of the nonlinearity in the model. Thus we often plot τk against θk∗ and
use the evident curvature to gain an impression of the degree of nonlinearity
present.
When p = 2, we can of course also plot the contours of S(θ1 , θ2 ). In the linear
case, these are ellipses; thus the degree of nonlinearity can be seen in the amount
by which these contours differ from elliptical shapes. Simultaneous confidence
sets are bounded precisely by these contour lines. The further away they are
from being circular, the greater the dependence between the two parameters
being estimated. In the same figure we can plot the so-called profile traces
(−1) (−2)
θb2 (θ1∗ ) versus θ1∗ and θb1 (θ2∗ ) versus θ2∗ . They intersect at θ,
b and their angle
shows how strongly the estimated parameters depend on each other. Further-
more, the profile traces intersect the contours in points where the latter have
horizontal or vertical tangents, respectively, since the gradient is perpendicular
to the contours. If the model is a linear one, S is quadratic and the profile
traces are straight lines (cf. Lemma 1.5.1). Thus we obtain a diagram like that
in Figure 1.13 (and the two regression lines there are exactly the two profile
traces). The contours and profile traces for the data in Figures 2.3 and 2.4 are
shown in Figure 2.5. We can see that the effects of nonlinearity are compara-
tively harmless in the reaction speed example, while they are fairly extreme in
the oxygen use example.

Figure 2.5: Contours and profile traces for the examples of reaction speed
(left) and oxygen use (right). The underlying data are shown in Figures 2.4
and 2.3.

Once p > 2, it is much more difficult to visualize the effect of nonlinearities and
of the dependence between the parameters. We could in principle compute the
64 Nonlinear and nonparametric methods

contours of the so-called profile likelihood

min S(θ)
θ;θj =θj∗ ,θk =θk∗

for all pairs (j, k), but this is often too expensive computationally. Instead,
(−k)
we can at least plot pairs of profile traces, plotting θbj (θk∗ ) against θk∗ and
(−j)
θb (θ∗ ) against θ∗ for all j < k. The curvature of these profile traces indicates
k j j
the extent of nonlinearity, while their angle shows the level of dependence be-
tween the corresponding parameter estimates. Furthermore, the profile traces
once more intersect the contours of the profile likelihood at points where the
tangents are horizontal or vertical, respectively. This gives us an impression of
where the contour lines should be.

2.3 Generalized linear models

2.3.1 Logistic regression

Many applications have a binary response variable Y (e.g. in medicine, patients

may be healed or dead), and the success probability P [Y = 1] depends on
explanatory variables. Models that have this probability as a linear function
of the explanatory variables are rarely meaningful, as the success probability
must always lie between 0 and 1. Logistic regression uses the model
  p
Pθ [Yi = 1] X
log = xij θj = xTi θ
Pθ [Yi = 0]
j=1

(if the covariates xi are also random, the LHS contains conditional probabilities
given xi ). Solving for Pθ [Yi = 1] gives us

exp(xTi θ)
Pθ [Yi = 1] = = P [U ≥ −xTi θ],
1 + exp(xTi θ)

where U has a so-called logistic distribution:

u
et
Z
exp(u)
P [U ≤ u] = P [U ≥ −u] = = dt.
1 + exp(u) −∞ (1 + et )2

This can be interpreted as a model in which latent variables Zi satisfy a linear

relationship with logistic errors εi :

Zi = xTi θ + εi .

Instead of Zi , it is the indicator Yi = 1[Zi ≥0] which is observed. Taking normally

distributed errors instead of logistically distributed ones gives us the so-called
probit model. These models differ only very slightly, but the logistic model is
computationally simpler.
2.3 Generalized linear models 65

The parameters in such a model are nearly always estimated using maximum
likelihood. It is easy to see that for any y ∈ {0, 1}, we have

Pθ [Yi = 1] y
 
Pθ [Yi = y] = Pθ [Yi = 0] = exp(y · xTi θ − log(1 + exp(xTi θ))).
Pθ [Yi = 0]

Thus (assuming the independence of the observations) we have the log-likelihood

n
X
`(θ) = (yi xTi θ − log(1 + exp(xTi θ))).
i=1

This is a concave function whose maximum is given by

n
X
(yi − Pθb [Yi = 1])xi = 0.
i=1

These equations are generally solved numerically using iterative methods. If the
P response values (yij ; 1 ≤ j ≤ ni ),
experimental condition xi has several observed
then `(θ) only depends on their sum yi+ = j yij and the total ni .
Confidence intervals and tests in this model depend on the asymptotic normal
approximation
b asymptotically
θ ∼ N (θ, V (θ)).
The asymptotic covariance matrix V (θ) of θ
b is the inverse of the Fisher infor-
mation (cf. Section 1.9 or Mathematical Statistics):
n n
−1
X  X exp(xTi θ)
xi xTi E 2
xi xTi

V (θ) = I(θ) = (yi − Pθ [Yi = 1]) = .
i=1 i=1
(1 + exp(xTi θ))2

Comparisons of two nested models with dimensions p and q < p are possible by
using double the log-likelihood quotient
(p)
2(`(θ
b b(q) )),
) − `(θ

which is known to asymptotically follow a chi-squared distribution with (p − q)

degrees of freedom.

2.3.2 General case

In the general case we have independent observations Yi and a density or dis-

tribution function of the form

pβi (yi ) = exp(yi βi + c(βi ))h(yi )

(a so-called exponential family). For such a model, we have

E [Yi ] = µ(βi ) = −c0 (βi )

66 Nonlinear and nonparametric methods

R
(which can be seen by taking the derivative of pβ (y)dy = 1 with respect
to β and switching the integration and differentiation steps). Many common
distributions take on the form of such an exponential family, including the
normal, binomial and Poisson distributions. For the normal distribution, we
have
µ µ2 y2
p(y) = exp(y 2 − 2 − 2 ).
σ 2σ 2σ
Thus if σ is known, we have an exponential family with parameters
µ 1
β= , c(β) = − σ 2 β 2 .
σ2 2
The binomial(n, p) distribution can be written as
 y  
p n n
p(y) = (1 − p) .
1−p y
This is an exponential family with
 
p
β = log , c(β) = −n log(1 + eβ ).
1−p

Similarly, the Poisson(λ) distributions form an exponential family with

β = log λ, c(β) = −eβ .

In a generalized linear model, the effect of the explanatory variables xi on the

observation Yi is described by a link between xi and the parameter βi for the
i-th observation. That is, we want there to be a function g for which

g(µ(βi )) = xTi θ,

i.e. a suitable transformation of βi is to be linear in the explanatory variables.

The function g is called a link function. If g is exactly µ−1 , we call it the canon-
ical link function. The linear model with normal errors, and logistic regression,
are examples of generalized linear models with a canonical link function.
The general case can be treated in much the same way as the logistic case, i.e.
taking the maximum likelihood estimate and using the asymptotic normality of
this estimate (or the asymptotic chi-squared distribution of the likelihood quo-
tient) to construct tests and confidence intervals. For details, see the literature.

2.4 Cox regression

In medical and technical applications, the response variable is often a survival

or failure time. We could in principle use a generalized linear model with
exponential or Gamma distributions. For practical purposes, though, the Cox
model has emerged as the front runner – and it is not bound to any specific
family of distributions. Let F be a distribution on the positive real numbers,
2.4 Cox regression 67

and denote its density by f . Then the failure rate (also known as hazard or
risk function) is defined as

1 f (t) d
λ(t) = lim P [t ≤ T ≤ t + h|T ≥ t] = = − log(1 − F (t)).
h→0 h 1 − F (t) dt

Thus the failure rate determines the distribution, namely by

Z t
F (t) = 1 − exp(− λ(u)du).
0

The Cox model now postulates that the shape of the i-th failure rate as a
function of the explanatory variables xi is

λi (t) = exp(xTi θ)λ0 (t),

where λ0 is a base rate. Of course this model does not allow an intercept to be
present, as a constant could be absorbed into λ0 . Raising the j-th component
of the explanatory variables by one unit leads to multiplying the failure rate
by the factor exp(θ j ) (homogeneously for all times). Because of this property,
this model is also known as “proportional hazard model”. It is easy to see that
any strictly monotonous and differentiable transformation of the survival times
transforms a Cox model into another Cox model with the same parameters and
a different λ0 . If λ0 is not specified, this merely means that the choice of a
timescale remains open.
The base rate λ0 naturally appears in the likelihood function, too, which means
that we cannot simply employ the maximum likelihood estimator. We shall
instead estimate the parameter θ by means of the so-called partial likelihood,
defined as
n
Y exp(xTi θ)
P T
.
i=1 j;tj ≥ti exp(xj θ)

The i-th factor is the conditional probability of failure for the i-th observed
unit in the interval [ti , ti + dt), given the failure of one of the units that was still
working just before time ti . Thus the only information on failure times used is
the order in which the failures occur.
Nearly all data of this type suffer from the added complication of censored ob-
servations, whose failure times Ti are not known exactly, but it is only known
that they are greater that some observed censoring time Ci . The reasons for this
include the termination of the study before the failure of all units, or patients
moving away or dying of other causes. In such situations the partial likelihood
can be defined in a similar way: We compute the product over all uncensored
observations, but in the denominator, we only add up the uncensored observa-
tions with tj ≥ ti and the censored observations with cj ≥ ti .
Once again, tests and confidence intervals depend on asymptotic considerations
that we shall not elaborate here.
68 Nonlinear and nonparametric methods

2.5 Nonparametric regression

Over the past 20-30 years much attention has been focused on procedures that
make no parametric assumptions about the shape of f in the model

yi = f (xi ) + εi ,

and only ask for f to be a smooth function. Weakening the assumptions in

this way is of course a great advantage, especially in explorative investigation,
where we aim at being as open as possible towards the data.
We separate the deterministic case where the xi are fixed and the stochastic
case where the xi are realizations of random variables that follow a distribution
G, irrespectively of the properties of the errors εi . A more general model
assumes that the data (Xi , Yi ) are i.i.d random vectors with some arbitrary
joint distribution, and that our intention is to estimate

f (x) = E [Yi | Xi = x] .

We can show that for any function g, the inequality

E (Yi − f (Xi ))2 ≤ E (Yi − g(Xi ))2

   

holds (as long as the second moments do in fact exist). This means that f (Xi )
is the best prediction of Yi from Xi in terms of mean squared error. The
conditional error variance

E (Yi − f (xi ))2 | Xi = xi ,

 

which corresponds to Var[εi ], then depends on xi in general.

All nonparametric estimates of f (x) essentially average the responses yi for
those i whose xi are in close proximity to x. There are, however, large differ-
ences in the method of averaging and in how to determine proximity to x.

2.5.1 Some procedures for the one-dimensional case

Kernel estimators:

As an estimate of f (x), we take the weighted mean of the yi in such a way that
the weight of yi is monotonely decreasing in the distance |x − xi |. The weights
are determined by a so-called kernel K and a bandwidth h > 0. A kernel is a
probability density which is symmetric around 0 and either has support [−1, 1]
or is very rapidly decreasing (as is e.g. the normal density). There are two
versions that differ considerably in the stochastic case. The Nadaraya-Watson
version is defined as
Pn
yi K((x − xi )/h)
f (x) = Pi=1
b n .
i=1 K((x − xi )/h)
2.5 Nonparametric regression 69

Figure 2.6: Weights of the Nadaraya-Watson (solid lines) and Gasser-Müller

(dashed lines) versions. The locations xi of the observations are marked by
circles on the x-axis, and the boundary points si by the dotted lines. All
weights have been scaled by a constant factor.

For the Gasser-Müller version we assume that we have ordered covariates

0 ≤ x1 < x2 < . . . < xn ≤ 1,

and we define s0 = −∞, si = (xi + xi+1 )/2 for 0 < i < n and sn = +∞. Our
estimate is then
n Z si
X 1
f (x) =
b yi K((x − u)/h)du.
si−1 h
i=1

The differences between these two estimators are quite important if the the
covariates xi are distributed unevenly, as we can see in Figure 2.6. We can
see that the Gasser-Müller version gives higher weights to observations whose
corresponding xi have an isolated location.
In both versions, the bandwidth h regulates the smoothness of fb: The larger
h is chosen to be, the smoother the estimate, but the less it can adapt to the
data.
If the covariates xi are chosen at random, the number of observations having a
significantly non-zero weight can vary quite strongly as x varies. For the first
estimator (Nadaraya-Watson), there may even be none of them. To avoid such
an occurrence, we can choose a variable bandwidth, e.g. to fix the number of
observations for which x − h ≤ xi ≤ x + h. The latter procedure leads to a
so-called nearest neighbour estimator.
Kernel estimators especially struggle to estimate f at the fringes, i.e. for x < h
and x > 1 − h. As for such points the averaged observations nearly all lie on
one side of x, systematic errors can occur.
70 Nonlinear and nonparametric methods

Local polynomial regression:

Here we no longer assume the function f to be locally constant. Instead, we

locally use a polynomial of degree p > 0, which is then fitted using weighted
least squares. This leads to
fb(x) = θb0 (x),
where
n
X p
X
θ(x)
b = arg min K((x − xi )/h)(yi − θj (xi − x)j )2 .
θ
i=1 j=0

It turns out that p is best chosen to be odd. In practice, p = 1 and p = 3 are

used most often. The advantages are seen quite clearly at the fringes.
Instead of a fixed bandwidth, we can once more use the nearest neighbour
version that ensures a fixed number of observations in [x − h, x + h]. This is the
key idea of the function loess in the statistical software packages S-Plus and
R.

Smoothing splines:

A smoothing spline is defined implicitly to be the solution of the minimization

problem
Xn Z 1
arg min 2
(yi − f (xi )) + λ f 00 (x)2 dx.
f 0
i=1

The first term, a sum over i, indicates the goodness of fit of f to the observa-
tions. The second term, the L2 norm of the second derivative, measures the
smoothness of f , and the parameter λ regulates the compromise that governs
the trade-off of the opposite aims of making both terms small simultaneously.
It can be shown that the solution of the minimization problem is a cubical spline
with nodes at the xi , and that this solution is furthermore linear on the fringe
intervals [0, x1 ] and [xn , 1]. If we take λ → 0, we obtain the spline that just
interpolates the data, and for λ → ∞ we obtain the least squares line. Thus
the role of λ corresponds to that of the bandwidth h.
To compute the smoothing spline,we choose a basis of the vector space of splines
with nodes xi . Computing the coefficients of the solution with respect to this
basis then leads to minimization of a quadratic function. If we are applying a
numerically stable and fast procedure, a good choice of basis is crucial. Here
the so-called B-splines have a good track record.

2.5.2 Bias-variance tradeoff

All non-parametric procedures contain a smoothing parameter that has a strong

influence on their behavior. For splines, this is λ, and for kernel estimators, it is
the bandwidth h. The purpose of this smoothing parameter becomes apparent
2.5 Nonparametric regression 71

if we look at the bias and variance of such an estimator. We can show that for
local polynomials with odd p,
h i
E fb(x) − f (x) ∼ const(K, p) hp+1 f (p+1) (x)

and −1
σ2

1 X
Var[fb(x)] ∼ const(K, p) ε K((x − xi )/h) .
nh nh
Here we assume the xi to be fixed (or if they are not, we condition on them),
and that the bandwidth h = hn satisfies hn → 0 and nhn → ∞. From the
above formulæ, we see that h should be small to obtain a small (absolute) bias,
and h should be large to obtain a small variance. Hence the title of this section.
One quantity that considers bias and variance simultaneously is the mean
squared error:
h i h i 1
E (fb(x) − f (x))2 = Var[fb(x)] + (E fb(x) − f (x))2 = O( ) + O(h2(p+1) ).
nh
The order of this is minimal if both its summands have the same order, i.e.
h = O(n−1/(2p+3) ). Such a choice leads to a mean squared error of the order
O(n−(2p+2)/(2p+3) ). Looking closely, we would prefer to have p as large as
possible; in practice, however, this is not quite right, as a larger value of p
also enlarges the constants and requires stronger assumptions on f . In any
case, the constants are the greatest problem when applying these results: they
contain unknown terms such as the derivatives of f , and the optimal choice of
h depends on the location x. Thus the data-dependent yet optimal choice of
bandwidth is quite a difficult problem.
We could also find the optimal kernel K (by minimizing const(K, p)). It turns
out, however, that the choice of K is of secondary importance, as almost all
continuous kernels are nearly as good as each other.
The bias of the Nadaraya-Watson kernel estimators has a more complicated
form, but their variance is asymptotically equivalent to that of the local poly-
nomials with p = 1. For the Gasser-Müller version, the bias is the same as for
local polynomials with p = 1, but instead the constant involved in the variance
is 1.5 the size as the one obtained for local polynomials of degree p = 1. We
even know the bias and variance for the smoothing splines; it is similar to that
of local polynomials of degree p = 3.

2.5.3 Curse of dimensionality

We can in principle generalize kernel estimators and local polynomials to more

than one dimension. In practice, though, they usually fail from dimensions 3 or
4 upwards. The reason for this is that high-dimensional spaces contain a large
amount of space and are not well covered by a finite set of points. In other words:
any two xi nearly always lie very far apart, making a reasonable compromise
on the bias-variance dilemma impossible. This “curse of dimensionality” is
72 Nonlinear and nonparametric methods

quite well illustrated by the following example: If the points xi are equally
distributed inside the cube [−1, 1]p , the proportion of points lying in the unit
sphere {x; ||x|| ≤ 1} is approximately the same as the probability of a xi lying
in this unit sphere, i.e. the volume of the unit sphere times 2−p . If p = 2, this
probability is 79%, for p = 5 it is still 16%, and for p = 10 only 0.25% ! On
√
the other hand, the ratio of diameters of sphere and cube is 1 : p, thus for
p = 10, we roughly have the ratio 1 : 3. So the assumption that f is constant or
linear on the unit sphere does not differ significantly from the assumption that
f is constant or linear on the entire cube – yet we still nearly always have too
few observations inside the unit sphere to estimate a constant or linear function
well.
Appendix A

Results from probability

theory

A.1 Computation of moments

Let the following be given:

• Random variables X, Y, Z
• Constants a, b
Computation of means:
• E [aX + b] = aE [X] + b.
• E [X + Y ] = E [X] + E [Y ].
• E [XY ] = E [X] · E [Y ], if X and Y are independent or at least uncorre-
lated.
Computation of variance and covariance:
• Var[aX + b] = a2 Var[X], and thus: σ[aX + b] = |a|σ[X]
• Cov[X, Y ] := E [(X − E [X])(Y − E [Y ])] = E [XY ]−E [X] E [Y ] = Cov[Y, X].
• Var[X + Y ] = Var X + Var Y + 2 Cov[X, Y ]
• In particular: Var[X +Y ] = Var[X]+Var[Y ] , if X and Y are independent
or at least uncorrelated.
• Cov[aX + bY, Z] = a Cov[X, Z] + b Cov[Y, Z].

1
2 Results from probability theory

Moments of random vectors:

Let Y be a n × 1-dimensional random vector. We define E [Y] to be the vector

with components E [yi ], and the covariance matrix of Y to be
 
Var[y1 ] Cov[y1 , y2 ] ... Cov[y1 , yn ]
 Cov[y2 , y1 ] Var[y2 ] ... Cov[y2 , yn ] 
Cov[Y] = (Cov(yi , yj )ij ) = 
 ...

... ... ... 
Cov[yn , y1 ] Cov[yn , y2 ] ... Var[yn ]

Thus a covariance matrix is always symmetric.

Let A be a fixed matrix of dimension m × n, and let b be a fixed vector of

dimension m × 1. Then as in the scalar case, we have:

• Cov[Y] = E YYT − E [Y] E [Y]T .

 

• E [AY + b] = AE [Y] + b,

• Cov[AY + b] = A · Cov[Y] · AT .

In particular, the latter rule tells us that for an arbitrary vector a of dimension
n × 1:
0 ≤ Var[aT Y] = Cov[aT Y] = aT Cov[Y]a.

In other words, every covariance matrix is positive semidefinite (and usually

even positive definite).

If Σ is an arbitrary positive semidefinite matrix, there exist matrices A satis-

fying AAT = Σ. We call each such matrix A a root of Σ, and write A = Σ1/2 .
Note that Σ1/2 is only defined up to multiplication by an orthogonal matrix.
Numerically speaking, the easiest way to find Σ1/2 uses the Cholesky decom-
position, which yields a lower triangular matrix. Some care must also be taken
when using Σ1/2 in computations, as e.g. (Σ−1 )1/2 = (Σ1/2 )−T .

If Y is a random vector with covariance matrix Σ and A = Σ1/2 , the rules

listed above imply that the covariance matrix of the random vector X = A−1 Y
is precisely the identity matrix, i.e. that the components of X are uncorrelated
and have variance 1.

A.2 The normal distribution

A.2.1 Univariate normal distribution

a) Density of the “standard normal distribution” N (0, 1):

1 1 2
ϕ(x) := √ e− 2 x
2π
A.2 The normal distribution 3

(This is a probability density, i.e. the integral of ϕ(x) over the real axis
is 1). The cumulative density function
Z x
Φ(x) = ϕ(y)dy
−∞

has no closed-form description, but it can be looked up in tables.

b) 3 reasons why the normal distribution is important (as an ideal
model):
(i) Simplicity (and the beauty of the ensuing theory).
(ii) Central limit theorem (and the elementary error hypothesis):
An error in measurement consists of many tiny yet independent el-
ementary errors which combine additively; thus their sum is ap-
proximately normally distributed. A multiplicative combina-
tion of elementary errors gives us a logarithmic normal distribu-
tion, which transforms into a normal distribution when logarithms
are taken.
(iii) Experience Many data sets are approximately normally distributed
(possibly only after a suitable transformation has been applied).
c) Why is the normal distribution so simple in some sense?
If f (x) is a density on R that satisfies

d log f (x) f0
= (x) = ax + b ,
dx f
1 2
it follows that f (x) = e 2 ax +bx+c . This quantity f 0 /f is quite important
in statistics. Thus in a certain sense, the normal distribution really is the
most simple continuous distribution on the whole real line!
d) Linear transformations of the standard normal distribution N (0, 1) :
Regard X ∼ N (0, 1) and the transformation x 7−→ y := µ + σx, taking
an arbitrary µ and some σ > 0. Then the distribution of Y := µ + σX is
the general normal distribution N (µ, σ 2 ).
Computing the density of Y : From

fX (x)dx = fY (y)dy with dy = σdx,

we conclude that
2
1 1 2 1 1 (y−µ)
√ e− 2 x dx = √ e− 2 σ2 dy.
2π 2π · σ
| {z }
density of Y :=µ+σX

e) Moments of the normal distribution:

4 Results from probability theory

Standard normal distribution: X ∼ N (0, 1).

Z ∞
E [X] = xϕ(x)dx = 0
Z−∞
∞ Z ∞
+∞
Var[X] = x · xϕ(x)dx = −xϕ(x) |−∞ + ϕ(x)dx = 1
−∞ | {z } −∞
0 | {z }
1
E X 3 = 0 (symmetry)
 
Z ∞ Z ∞
+∞
 4 3 3
E X = x · xϕ(x)dx = −x ϕ(x) |−∞ + 3x2 ϕ(x)dx = 3.
−∞ −∞

General normal distribution: Y ∼ N (µ, σ 2 ).

Mean: E [Y ] = µ
Variance: Var[Y ] = σ 2
E (Y − µ)3  /σ 3 = 0

Skewness (standardized 3rd moment): γ1 =
Kurtosis, excess (stand. 4th moment): γ2 = E (Y − µ)4 /σ 4 − 3 = 0

f) Shape of the normal distribution

Figure A.1 gives a graphical representation of the density of the standard
normal distribution N (0, 1).

Figure A.1: Density of the standard normal distribution

Some values of the cumulative distribution function:

x 0 0.6745 1 1.64 1.96 2.58 3.3
1 3 5
Φ(x) 2 4 84 % ≈ 6 95 % 97.5 99.5 99.95 %
The density f (x) of the normal distribution goes to 0 more and more
quickly in relative terms when x −→ ±∞, i.e.|f 0 |/f −→ ∞. Although the
normal distribution has a positive density all the way from −∞ to +∞,
it is in practice a very “short-tailed” distribution that all but disappears
beyond µ ± 3σ or µ ± 4σ. This contradicts most empirical distributions
of measurement errors. Thus the “normal approximation” is usually
A.2 The normal distribution 5

only helpful in the middle of the distribution, e.g. up to µ ± 2σ or

µ ± 2.5σ.

A.2.2 Multivariate normal distribution

Let Y1 , Y2 . . . Yn be independent standard normal random variables. Then their

joint density is simply the product of their individual densities, that is,

fY (y) = (2π)−n/2 exp(− 21 yT y).

This is the n-dimensional standard normal distribution. It exhibits spherical

symmetry. The general n-dimensional normal distribution is defined as the
distribution of a vector X that is obtained by a linear transformation of a
standard normally distributed n-dimensional vector:

X = AY + µ.

Here µ is an (n×1)-dimensional vector, and A is an (n×n)-dimensional matrix.

If A is a singular matrix, we call the resulting distribution degenerate. We shall
first discuss the non-degenerate case, however. In this case the distribution of
X has the density

fX (x) = (2π)−n/2 (| det A|)−1 exp(− 21 (A−1 (x − µ))T (A−1 (x − µ)))

= (2π)−n/2 (det Σ)−1/2 exp(− 21 (x − µ)T Σ−1 (x − µ)),

where Σ = AAT . According to the computational laws for means, we have

E [X] = AE [Y] + µ = µ

and

Cov[X] = E (X − µ)(X − µ)T = E AYYT AT = AE YYT AT = AAT = Σ.

     

Thus as in the univariate case, the normal distribution is determined entirely

by the first two moments, and we therefore write N n (µ, Σ). The covariance
matrix Σ can be an arbitrary positive definite matrix. The joint density has its
maximum at µ and is constant on similar ellipsoids centred on µ (the principal
axes of ellipsoids are given by the eigenvectors of Σ).
The following results both follow straight from this definition:
Theorem A.2.1. Uncorrelated and jointly normally distributed random vari-
ables are independent. More precisely: if X ∼ N n (µ, Σ) and if Σij = 0 for all
i ∈ I and j ∈ J for two disjoint index sets I, J ⊂ {1, . . . , n}, then the random
vectors (Xi , i ∈ I) and (Xj , j ∈ J) are independent of each other.
Proof The joint density can be decomposed as the product of the individual
densities.
Theorem A.2.2. Standard normally distributed random variables keep their
distribution when subjected to orthogonal transformations.
6 Results from probability theory

With some more effort, we can even show the following:

Theorem A.2.3. Linear combinations of jointly normally distributed random
variables are themselves jointly normally distributed.
Proof For linear combinations AX+b that use an n×n-dimensional matrix A,
the claim follows immediately from the definition of the n-dimensional normal
[Link] less than n linear combinations are involved, we assume without
loss of generality that X has a standard normal distribution. For the sake of
simplicity, we take a single linear combination aT X = ni=1 ai Xi and assume
P
the length of a is 1. Then we choose an orthogonal matrix A whose first row is
equal to aT , i.e. aT X is the first component of AX. This makes it clear that
aT X is normally distributed, as the components of AX are independent and
normally distributed.
Thus we can understand the structure of the degenerate normal distribution:
Theorem A.2.4. If X has a degenerate n-dimensional normal distribution,
then it has r components (Xi1 , Xi2 , . . . Xir ) which have a non-degenerate r-
dimensional distribution, while the remaining n − r components can be written
as linear combinations of these. Here r is the rank of Σ.
We can furthermore use Theorem A.2.3 to show that in particular:
Corollary A.2.1. The n marginal distributions of an n-dimensional normally
distributed random variable themselves all follow normal distributions.
The converse of this corollary is not true, however. Even if all marginal dis-
tributions of a random variable are normal, the joint distribution need not
necessarily be a normal distribution!

Figure A.2: Despite normal marginal distributions, the joint distribution need
not be a normal one.

Consider the following example:

Let U be a univariate standard normal random variable. Let X = Y = U , which
implies that the entire distribution of (X, Y ) is concentrated on the diagonal
and the marginal distributions are certainly normal. If we now cut away part of
the joint distribution in a symmetric way (as illustrated in Figure A.2) and place
the cut-off piece on the other diagonal, we no longer have a joint normal
distribution, while the projections that gives us the marginal distributions are
still the same, i.e. they are still normally distributed!
We should furthermore note that
Cut-off point −→ ∞ =⇒ Correlation +1
Cut-off point −→ 0 =⇒ Correlation −1,
A.2 The normal distribution 7

and between these two extremes the correlation varies in a continuous manner,
including taking on the value 0 at some point. However, the two variables here
are always dependent. Thus we also see that without joint normality, lack of
correlation does not mean independence.

A.2.3 Chi-squared, t and F distributions

These distributions are derived from the normal distribution and have an im-
portant role in various tests in regression. Let
X1 , X2 , . . . , Xm , Y1 , Y2 , . . . , Yn independent ∼ N (0, 1).

Then the distribution of

m
X
Zm = Xi2
i=1
is called the chi-squared distribution with m degrees of freedom (written
χ2m ). In particular, we have E [Zm ] = m, V ar[Zm ] = 2m and
 
Zm − m m→∞
L √ −→ N (0, 1).
2m
(L stands for (distribution) law, and by convergence we mean weak conver-
gence.)
The distribution of
X1
Vn = q P
1 n 2
n i=1 Yi

is called t distribution with n degrees of freedom (written tn ). In particular,

t1 is the (standard) “Cauchy distribution”, and
n→∞
L(Vn ) −→ N (0, 1).

The distribution of Pm
1
X2
Wm,n = m
1 Pni=1 2i
n i=1 Yi
is called F distribution with m degrees of freedom in the enumerator and n
degrees of freedom in the denominator (and is written Fm,n ). In particular, we
have
1
L(Wm,n ) → χ2m (m fest, n → ∞)
m
and
L(Wm,n ) → 1, (m → ∞, n → ∞).

On the computation of the densities of these distributions

We start with χ21 ; that is, let X ∼ N (0, 1) be given and L(X 2 ) sought. The
solution:
√ √ √ √ √
P [X 2 ≤ c] = P [|X| ≤ c] = P [− c ≤ X ≤ c] = Φ( c) − Φ(− c)
8 Results from probability theory

Differentiating gives us the corresponding density.

The formula for the density of χ2m can be obtained by repeated application of
the convolution formula (density of the sum of independent random variables).
The formula for the t and F distributions is based on the following line of
thought: let U and V > 0 be two independent random variables with densities
fU and fV . Then we have the identity
Z Z Z ∞
U
P [ ≤ x] = fU (u)fV (v)dudv = fV (v)FU (xv)dv.
V {u≤xv} 0

Taking derivatives by x on the right-hand side, we then obtain the density of

U/V : Z ∞
fU/V (x) = fV (v)fU (xv)vdv.
0

As actually computing these integrals is a time-consuming process that provides

no new insights, we omit it here.
Finally we mention a result than can sometimes be quite useful:
Lemma A.2.1. If a random vector X has a N n (µ, Σ) distribution, then if
follows that (X − µ)T Σ−1 (X − µ) has a chi-squared distribution with n degrees
of freedom.
Proof We can write X as µ + AY, where Y follows a N n (0, 1n ) distribution
and AAT = Σ. Then the quadratic form in X is exactly YT Y =
P 2
Yi , and
the claim thus follows by the definition of the chi-squared distribution.
Appendix B

Literature

Agresti, A. (1990). Categorical Data Analysis, Wiley, N. Y. Generalized linear

models.

Bates, D. M. and Watts, D. G. (1988). Nonlinear regression analysis and its

applications, Wiley N. Y. Standard book on nonlinear regression.

Bowman, A. W. and Azzalini, A. (1997) Applied Smoothing Techniques for

Data Analysis, Oxford University Press. Nonparametric regression.

Chambers, J. M. and Hastie, T. J. (1992). Statistical Models in S, Wadsworth

& Brooks/Cole. Discussion of regression and other models with a view to
implementation in S and R.

Cook, R. D. and Weisberg, S. (1994). An Introduction to Regression Graphics,

Wiley, N. Y. Graphical and exploratory methods for regression models.

Cox, D. R. and Oakes, D. (1984). Analysis of Survival Data, Chapman and

Hall, London. Survival times, Cox regression.

Daniel, C. and Wood, F. S. (1980). Fitting Equations to Data, 2nd edn, Wiley,
N. Y. A classic for data analysis and applications.

Dobson, A. J. (1991). An Introduction to Generalized Linear Models, Chapman

and Hall, London. Short and clear introduction to linear and generalized
linear models.

Fuller, W. A. (1987). Measurement Error Models, Wiley, N. Y. Specialist

literature on the topic of “errors in variables”.

Hampel, F. R., Ronchetti, E. M., Rousseeuw, P. J. and Stahel, W. A. (1986).

Robust Statistics: The Approach Based on Influence Functions, Wiley,
N. Y. Robustness in general, with 2 chapters on robust regression.

Lindsey, J. K. (1997). Applying Generalized Linear Models, Springer, N. Y.

Generalized linear models, with more complex applications.

9
10 Literature

McCullagh, P. and Nelder, J. A. (1989, 1997). Generalized Linear Models, 2nd

edn, Chapman-Hall, London. Standard book on generalized linear models.
Rousseeuw, P. J. and Leroy, A. M. (1987). Robust Regression & Outlier De-
tection, Wiley, N. Y. Robust regression, easier to read than Hampel et al.
Ryan, T. P. (1997). Modern Regression Methods, Wiley, N. Y. Especially for
reference purposes.

Seber, G. A. F. (1980). The Linear Hypothesis: A General Theory, 2nd edn,

Charles Griffin, London. Mathematically elegant and concise treatment.
Seber, G. A. F. and Wild, C. J. (1989). Nonlinear Regression, Wiley, N. Y.
Somewhat more mathematical than Bates and Watts, quite extensive.

Sen, A. and Srivastava, M. (1990). Regression Analysis; Theory, Methods, and

Applications, Springer-Verlag, N. Y. Textbook and reference work, similar
to this class.

Simonoff, J. S. (1996). Smoothing Methods in Statistics, Springer Series in

Statistics, Springer-Verlag, N. Y. Nonparametric regression.
Venables, W. N. and Ripley, B. D. (1999). Modern Applied Statistics with S-
Plus, 3rd edn, Springer-Verlag, N. Y. Overview of all prevalent statistical
procedures and their implementation in S and R.

Weisberg, S. (1990). Applied Linear Regression, 2nd edn, Wiley, N. Y. Text-

book, underlines applications.