Notes 1

1 Simple Linear Regression
Text: RPD, Chapter 1

Problems:
1.1 Statistical Model

In simple linear regression, the model contains a random dependent (or response or outcome or
end point) variable Y , that is hypothesized to be associated with an independent (or predictor
or explanatory) variable X. The simple linear regression model specifies that the mean, or
expected value of Y is a linear function of the level of X. Further, X is presumed to be set by
the experimenter (as in controlled experiments) or known in advance to the activity generating the
response Y . The experiment consists of obtaining a sample of n pairs (Xi , Yi ) from a population
of such pairs (or nature). The model with respect to the mean is:
E[Yi ] = β 0 + β 1 Xi
where β0 is the mean of when when X = 0 (assuming this is a reasonable level of X), or more
generally the Y –intercept of the regression line; β1 is the change in the mean of Y as X increases
by a single unit, or the slope of the regression line. Note that in practice β0 and β1 are unknown
parameters that will be estimated from sample data.
Individual measurements are assumed to be independent, and normally distributed around the
mean at their corresponding X level, with standard deviation σ. This can be stated as below:
Yi = β 0 + β 1 Xi + ε i εi ∼ N ID(0, σ 2 )
where εi is a random error term and N ID(0, σ 2 ) means normally and independently distributed
with mean 0, and variance σ 2 .
1.1.1 Examples
The following two examples are based on applications of regression in pharmacodynamics and
microeconomics.
Example 1 – Pharmacodynamics of LSD
The following data were published by J.G. Wagner, et al, in the 1968 article: “Correlation
of Performance Test Scores with ‘Tissue Concentration’ of Lysergic Acid Diethylamide in Human
Subjects,” (Clinical Pharmacology & Therapeutics, 9:635–638).
Y — Mean score on math test (relative to control) for a group of five male volunteers.
X — Mean tissue concentration of LSD among the volunteers.
A sample of n = 7 points were selected, with Xi and Yi being measured at each point in time.
These 7 observations are treated as a sample from all possible realizations from this experiment.
The parameter β1 represents the systematic change in mean score as tissue concentration increases
by one unit, and β0 represents the true mean score when the concentration is 0. The data are given
in Table 1.1.1.
i Xi Yi
1 1.17 78.93
2 2.97 58.20
3 3.26 67.47
4 4.69 37.47
5 5.83 45.65
6 6.00 32.92
7 6.41 29.97
Table 1: LSD concentrations and math scores – Wagner, et al (1968)

.
Example 2 – Estimating Cost Functions of a Hosiery Mill
The following (approximate) data were published by Joel Dean, in the 1941 article: “Statistical
Cost Functions of a Hosiery Mill,” (Studies in Business Administration, vol. 14, no. 3).
Y — Monthly total production cost (in $1000s).
X — Monthly output (in thousands of dozens produced).
A sample of n = 48 months of data were used, with Xi and Yi being measured for each month.
The parameter β1 represents the change in mean cost per unit increase in output (unit variable
cost), and β0 represents the true mean cost when the output is 0, without shutting plant (fixed
cost). The data are given in Table 1.1.1 (the order is arbitrary as the data are printed in table
form, and were obtained from visual inspection/approximation of plot).
1.1.2 Generating Data from the Model

To generate data from the model using a computer program, use the following steps:
1. Specify the model parameters: β0 , β1 , σ
2. Specify the levels of Xi , i = 1, . . . , n. This can be done easily with do loops or by brute force.
3. Obtain n standard normal errors Zi ∼ N (0, 1), i = 1, . . . , n. Statistical routines have them
built in, or transformations of uniform random variates can be obtained.
4. Obtain random response Yi = β0 + (β1 Xi ) + (σZi ), i = 1, . . . , n.
5. For the case of random Xi , these steps are first completed for Xi in 2), then continued for Yi .
The Zi used for Xi must be independent of that used for Yi .
1.2 Least Squares Estimation

The parameters β0 and β1 can take on any values in the range (∞, ∞), and σ can take on any values
in the range [0, ∞) (if it takes on 0, then the model is deterministic, and not probabilistic). The
most common choice of estimated regression equation (line in the case of simple linear regression),
is to choose the line that minimizes the sum of squared vertical distances between the observed
i Xi Yi i Xi Yi i Xi Yi
1 46.75 92.64 17 36.54 91.56 33 32.26 66.71
2 42.18 88.81 18 37.03 84.12 34 30.97 64.37
3 41.86 86.44 19 36.60 81.22 35 28.20 56.09
4 43.29 88.80 20 37.58 83.35 36 24.58 50.25
5 42.12 86.38 21 36.48 82.29 37 20.25 43.65
6 41.78 89.87 22 38.25 80.92 38 17.09 38.01
7 41.47 88.53 23 37.26 76.92 39 14.35 31.40
8 42.21 91.11 24 38.59 78.35 40 13.11 29.45
9 41.03 81.22 25 40.89 74.57 41 9.50 29.02
10 39.84 83.72 26 37.66 71.60 42 9.74 19.05
11 39.15 84.54 27 38.79 65.64 43 9.34 20.36
12 39.20 85.66 28 38.78 62.09 44 7.51 17.68
13 39.52 85.87 29 36.70 61.66 45 8.35 19.23
14 38.05 85.23 30 35.10 77.14 46 6.25 14.92
15 39.16 87.75 31 33.75 75.47 47 5.45 11.44
16 38.59 92.62 32 34.29 70.37 48 3.79 12.69
Table 2: Production costs and Output – Dean (1941)

.
responses (Yi ) and the fitted regression line (Ŷi = β̂0 + β̂1 Xi ), where β̂0 and β̂1 are sample based
estimates of β0 and β1 , respectively.
Mathematically, we can label the error sum of squares as the sum of squared distances between
the observed data and their mean values based on the model:
n
X n
X
Q = (Yi − E(Yi ))2 = (Yi − (β0 + β1 Xi ))2
i=1 i=1
The least squares estimates of β0 and β1 that minimize Q, which are obtained by taking derivatives,
setting them equal to 0, and solving for β̂0 and β̂1 .
n
X
∂Q
= 2 (Yi − β0 − β1 Xi )(−1) = 0
∂β0 i=1
n
X n
X n
X
⇒ (Yi − β0 − β1 Xi ) = Yi − nβ0 − β1 Xi = 0 (1)
i=1 i=1 i=1
n
X
∂Q
= 2 (Yi − β0 − β1 Xi )(−Xi ) = 0
∂β1 i=1
n
X n
X n
X n
X
⇒ (Yi − β0 − β1 Xi )Xi = X i Yi − β 0 Xi − β 1 Xi2 = 0 (2)
i=1 i=1 i=1 i=1
From equations (1) and (2) we obtain the so–called “normal equations”:
n
X n
X
nβ̂0 + β̂1 Xi = Yi (3)
i=1 i=1
n
X n
X n
X
β̂0 Xi + β̂1 Xi2 = X i Yi (4)
i=1 i=1 i=1
Pn
Multipliying equation (3) by i=1 Xi and equation (4) by n, we obtain the following two equations:
n
X n
X n
X n
X
2
nβ̂0 Xi + β̂1 ( Xi ) = ( Xi )( Yi ) (5)
i=1 i=1 i=1 i=1
n
X n
X n
X
nβ̂0 Xi + nβ̂1 Xi2 = n X i Yi (6)
i=1 i=1 i=1
Subtracting equation (5) from (6), we get:
n
X n
X n
X n
X n
X
β̂1 (n Xi2 − ( X i )2 ) = n X i Yi − ( Xi )( Yi )
i=1 i=1 i=1 i=1 i=1
Pn Pn
Pn ( Xi )( i=1 Yi ) Pn
i=1 Xi Yi −
i=1
Pn n 2 i=1 (Xi − X)(Yi − Y)
⇒ β̂1 = Pn = Pn 2
(7)
( i=1 Xi )
2 i=1 (Xi − X)
i=1 Xi − n
Now, from equation (1), we get:
n
X n
X
nβ̂0 = Yi − β̂1 Xi ⇒ β̂0 = Y − β̂1 X (8)
i=1 i=1
and the estimated (or fitted or prediction) equation (Ŷi ):
Ŷi = β̂0 + β̂1 Xi i = 1, . . . , n (9)
The residuals are defined as the difference between the observed responses (Yi ) and their predicted
values (Ŷi ), where the residuals are denoted as ei (they are estimates of εi ):
ei = Yi − Ŷi = Yi − (β̂0 + β̂1 Xi ) i = 1, . . . , n (10)
The residuals sum to 0 for this model:

ei = (Yi − Ŷi ) = Yi − (β̂0 + β̂1 Xi )
= Yi − {[Y − β̂1 X] + β̂1 Xi ]}

n
X n
X n
X n
X n
X n
X n
X
⇒ ei = Yi − {nY − nβ̂1 X + β̂1 Xi } = Yi − { Yi − β̂1 Xi + β̂1 Xi } = 0
i=1 i=1 i=1 i=1 i=1 i=1 i=1
1.2.1 Examples
Numerical results for the two examples desribed before are given below.
This dataset has n = 7 observations with a mean LSD tissue content of X = 4.3329, and a
mean math score of Y = 50.0871.
n
X n
X n
X n
X n
X
Xi = 30.33 Xi2 = 153.8905 Yi = 350.61 Yi2 = 19639.2365 Xi Yi = 1316.6558
i=1 i=1 i=1 i=1 i=1
Pn Pn
Pn ( Xi )( i=1 Yi ) (30.33)(350.61)
i=1 Xi Yi −
i=1
1316.6558 −
β̂1 = Pn n = 7
=
Pn ( i=1 Xi )2 (30.33)2
2
i=1 Xi −
153.8905 − 7
n
1316.6558 − 1519.1430 −202.4872

= = = −9.0091
153.8905 − 131.4146 22.4759
β̂0 = Y − β̂1 X = 50.0871 − (−9.0091)(4.3329) = 89.1226
Ŷi = β̂0 + β̂1 Xi = 89.1226 − 9.0091Xi i = 1, . . . , 7
ei = Yi − Ŷi = Yi − (89.1226 − 9.0091Xi ) i = 1, . . . , 7
Table 1.2.1 gives the raw data, their fitted values, and residuals.
i Xi Yi Ŷi ei
1 1.17 78.93 78.5820 0.3480
2 2.97 58.20 62.3656 -4.1656
3 3.26 67.47 59.7529 7.7171
4 4.69 37.47 46.8699 -9.3999
5 5.83 45.65 36.5995 9.0505
6 6.00 32.92 35.0680 -2.1480
7 6.41 29.97 31.3743 -1.4043
Table 3: LSD concentrations, math scores, fitted values and residuals – Wagner, et al (1968)
.
A plot of the data and regression line are given in Figure 1.
Example 2 – Estimating Cost Function of a Hosiery Mill
This dataset has n = 48 observations with a mean output (in 1000s of dozens) of X = 31.0673,
and a mean monthly cost (in $1000s) of Y = 65.4329.
n
X n
X n
X n
X n
X
Xi = 1491.23 Xi2 = 54067.42 Yi = 3140.78 Yi2 = 238424.46 Xi Yi = 113095.80
i=1 i=1 i=1 i=1 i=1
Pn Pn
Pn ( Xi )( i=1 Yi ) (1491.23)(3140.78)
i=1 Xi Yi −
i=1
113095.80 −
β̂1 = Pn n 2 = 48
=
Pn ( i=1 Xi ) (1491.23)2
2
i=1 Xi −
54067.42 − 48
n
113095.80 − 97575.53 15520.27

= = = 2.0055
54067.42 − 46328.48 7738.94
score Y=Math Score X=Tissue LSD Concentration
80
score
80
70
70
60
60
50
50
40
40
33 00
20
1 2 3 4 5 6 7
conc
Figure 1: Regression of math score on LSD concentration (Wagner, et al, 1968)

β̂0 = Y − β̂1 X = 65.4329 − (2.0055)(31.0673) = 3.1274
Ŷi = β̂0 + β̂1 Xi = 3.1274 + 2.0055Xi i = 1, . . . , 48
ei = Yi − Ŷi = Yi − (3.1274 + 2.0055Xi ) i = 1, . . . , 48
Table 1.2.1 gives the raw data, their fitted values, and residuals.
A plot of the data and regression line are given in Figure 2.
cost Y=Math Score X=Tissue LSD Concentration

100
score
80
90
8 07 0
70
60
60
50
50
40
40
30
30
20
1 02 0
01 21 0 3 20 4 30 5 40 6 5 07
scioznec
Figure 2: Estimated cost function for hosiery mill (Dean, 1941)
1.3 Analysis of Variance

The total variation in the response (Y ) can be partitioned into parts that are attributable to various
sources. The response Yi can be written as follows:
Yi = Ŷi + ei i = 1, . . . , n
We start with the total (uncorrected) sum of squares for Y:

n
X n
X n
X n
X n
X
SS(TOTAL UNCORRECTED) = Yi2 = (Ŷi + ei )2 = Ŷi2 + e2i + 2 Ŷi ei
i=1 i=1 i=1 i=1 i=1
i Xi Yi Ŷi ei
1 46.75 92.64 96.88 -4.24
2 42.18 88.81 87.72 1.09
3 41.86 86.44 87.08 -0.64
4 43.29 88.80 89.95 -1.15
5 42.12 86.38 87.60 -1.22
6 41.78 89.87 86.92 2.95
7 41.47 88.53 86.30 2.23
8 42.21 91.11 87.78 3.33
9 41.03 81.22 85.41 -4.19
10 39.84 83.72 83.03 0.69
11 39.15 84.54 81.64 2.90
12 39.20 85.66 81.74 3.92
13 39.52 85.87 82.38 3.49
14 38.05 85.23 79.44 5.79
15 39.16 87.75 81.66 6.09
16 38.59 92.62 80.52 12.10
17 36.54 91.56 76.41 15.15
18 37.03 84.12 77.39 6.73
19 36.60 81.22 76.53 4.69
20 37.58 83.35 78.49 4.86
21 36.48 82.29 76.29 6.00
22 38.25 80.92 79.84 1.08
23 37.26 76.92 77.85 -0.93
24 38.59 78.35 80.52 -2.17
25 40.89 74.57 85.13 -10.56
26 37.66 71.60 78.65 -7.05
27 38.79 65.64 80.92 -15.28
28 38.78 62.09 80.90 -18.81
29 36.70 61.66 76.73 -15.07
30 35.10 77.14 73.52 3.62
31 33.75 75.47 70.81 4.66
32 34.29 70.37 71.90 -1.53
33 32.26 66.71 67.82 -1.11
34 30.97 64.37 65.24 -0.87
35 28.20 56.09 59.68 -3.59
36 24.58 50.25 52.42 -2.17
37 20.25 43.65 43.74 -0.09
38 17.09 38.01 37.40 0.61
39 14.35 31.40 31.91 -0.51
40 13.11 29.45 29.42 0.03
41 9.50 29.02 22.18 6.84
42 9.74 19.05 22.66 -3.61
43 9.34 20.36 21.86 -1.50
44 7.51 17.68 18.19 -0.51
45 8.35 19.23 19.87 -0.64
46 6.25 14.92 15.66 -0.74
47 5.45 11.44 14.06 -2.62
48 3.79 12.69 10.73 1.96
Table 4: Approximated Monthly Outputs, total costs, fitted values and residuals – Dean (1941)
.
Here is a proof that the final term on the right-hand side is 0 (which is very easy in matrix
algebra):
Pn Pn
i=1 (Xi − X)(Yi − Y) i=1 (Xi − X)(Yi − Y)
Ŷi = β̂0 + β̂1 Xi = (Y − Pn 2
X) + Xi Pn 2
(11)
i=1 (Xi − X) i=1 (Xi − X)
Pn Pn
i=1 (Xi − X)(Yi − Y) i=1 (Xi − X)(Yi − Y)
ei = Yi − Ŷi = Yi − (Y − Pn 2
X) + Xi Pn 2
(12)
i=1 (Xi − X) i=1 (Xi − X)
Combining equations (11) and (12), we get:

Pn Pn
i=1 (Xi − X)(Yi − Y ) i=1 (Xi − X)(Yi − Y )
ei Ŷi = Yi Y − X Pn 2
+ X Y
i i P n 2
−
i=1 (Xi − X) i=1 (Xi − X)
Pn 2 Pn 2
i=1 (Xi − X)(Yi − Y ) 2 i=1 (Xi − X)(Yi − Y )
Y −X Pn 2
− X i P n 2
−
i=1 (Xi − X) i=1 (Xi − X)
Pn Pn
i=1 (Xi − X)(Yi − Y ) i=1 (Xi − X)(Yi − Y )
2Xi Y − X Pn 2
Pn 2
i=1 (Xi − X) i=1 (Xi − X)
Pn Pn
i=1 (Xi − X)(Yi − Y ) i=1 (Xi − X)(Yi − Y ) 2
= Yi Y − X Pn 2
+ X Y
i i P n 2
−Y −
i=1 (X i − X) i=1 (X i − X)
Pn 2 Pn
2 i=1 (Xi − X)(Yi − Y) i=1 (Xi − X)(Yi − Y)
X Pn 2
+ 2Y X Pn 2
−
i=1 (Xi − X) i=1 (Xi − X)
Pn 2 Pn 2
i=1 (Xi − X)(Yi − Y) i=1 (Xi − X)(Yi − Y)
Xi2 Pn 2
− 2Xi X Pn 2
i=1 (Xi − X) i=1 (Xi − X)
Pn 2
i=1 (Xi − X)(Yi − Y) 2
= Pn 2
(2Xi X − Xi2 − X )+
i=1 (Xi − X)
Pn
i=1 (Xi − X)(Yi − Y) 2
Pn 2
(−Yi X + Xi Yi + 2Y X − 2Xi Y ) + Yi Y − Y
i=1 (Xi − X)
Now summing ei Ŷi over all observations:

X
n Pn 2 X
n X
n
i=1 (Xi − X)(Yi − Y) 2
ei Ŷi = Pn 2
(2X Xi − Xi2 − nX )+
i=1 i=1 (Xi − X) i=1 i=1
Pn X n Xn Xn Xn
i=1 (Xi − X)(Yi − Y ) 2
Pn 2
( X i Y i − X Y i + 2nXY − 2Y X i ) + Y Yi − nY
i=1 (Xi − X) i=1 i=1 i=1 i=1
Pn 2 n
X
i=1 (Xi − X)(Yi − Y) 2 2
= Pn 2
(2nX − Xi2 − nX )+
i=1 (Xi − X) i=1
Pn X n
i=1 (Xi − X)(Yi − Y ) 2 2
Pn 2
( Xi Yi − nXY + 2nXY − 2nXY ) + nY − nY
i=1 (X i − X) i=1
Pn 2 n
X Pn Xn
i=1 (Xi − X)(Yi − Y) 2 i=1 (Xi − X)(Yi − Y)
= Pn 2
(nX − Xi2 ) + Pn 2
( Xi Yi − nXY )
i=1 (Xi − X) i=1 i=1 (Xi − X) i=1
Now making use of the following two facts:

X
n X
n
2 X
n X
n
2
(Xi − X)2 = Xi2 + nX − 2X Xi = Xi2 − nX
i=1 i=1 i=1 i=1
n
X n
X n
X n
X n
X
(Xi − X)(Yi − Y ) = Xi Yi + nXY − X Yi − Y Xi = Xi Yi − nXY
i=1 i=1 i=1 i=1 i=1
We obtain the desired result:

n
X Pn 2 n
X Pn X
n
i=1 (Xi − X)(Yi − Y) 2 i=1 (Xi − X)(Yi − Y)
ei Ŷi = Pn 2
(− (Xi −X) )+ Pn 2
( (Xi −X)(Yi −Y ))
i=1 i=1 (Xi − X) i=1 i=1 (Xi − X) i=1
Pn Pn
( i=1 (Xi − X)(Yi − Y ))2 ( i=1 (Xi − X)(Yi − Y ))2
= − Pn 2
+ Pn 2
= 0
i=1 (Xi − X) i=1 (Xi − X)
Thus we can partition the total (uncorrected) sum of squares into the sum of squares of the
predicted values (the model sum of squares) and the sum of squares of the errors (the residual sum
of squares).
n
X n
X n
X
Yi2 = Ŷi2 + e2i
i=1 i=1 i=1
SS(TOTAL UNCORRECTED) = SS(MODEL) + SS(RESIDUAL)
The computational formulas are obtained as follows:

n
X n
X
SS(Model) = Ŷi2 = (β̂0 + β̂1 Xi )2
i=1 i=1
X
n X
n X
n X
n
= nβ̂02 + 2β̂0 β̂1 Xi + β̂12 Xi2 = 2
n(Y − β̂1 X) + 2(Y − β̂1 X)β̂1 Xi + β̂12 Xi2
i=1 i=1 i=1 i=1
n
X
2 2 2
= nY + nβ̂12 X − 2nβ̂1 Y X + 2nβ̂1 Y X − 2β̂12 nX + β̂12 Xi2
i=1
2 2 X
n
2 X
n
2
= nY − nβ̂12 X + β̂12 Xi2 = nY + β̂12 ( Xi2 − nX )
i=1 i=1
n
X
2
= nY + β̂12 (Xi − X)2
i=1
n
X n
X
SS(RESIDUAL) = e2i = (Yi − Ŷi )2 = SS(TOTAL UNCORRECTED) − SS(MODEL)
i=1 i=1
The total (uncorrected) sum of squares is of little interest by itself, since it depends on the
level of the data, but not the variability (around the mean). The total (corrected) sum of squares
measures the sum of the squared deviations around the mean.
n
X n
X 2
SS(TOTAL CORRECTED) = (Yi − Y )2 = Yi2 − nY
i=1 i=1
n
X n
X
2
= (SS(MODEL) − nY ) + SS(RESIDUAL) = β̂12 (Xi − X)2 + e2i
i=1 i=1
= SS(REGRESSION) + SS(RESIDUAL)
Note that the model sum of squares considers both β0 and β1 , while the regression sum of
squares considers only the slope β1 .
For a general regression model, with p independent variables, we have the following model:
Yi = β0 + β1 Xi1 + · · · βp Xip i = 1, . . . , n
which contains p0 = p + 1 model parameters. The Analysis of Variance is given in Table 1.3,
which contains sources of variation, their degrees of freedom, and sums of squares.
Source of Degrees of Sum of

Variation Freedom Squares
Pn 2
Total (Uncorrected) n i=1 Yi
2
Correction Factor 1 nY
Pn 2 Pn 2 2
Total (Corrected) n−1 i=1 (Yi − Y ) = i=1 Yi − nY
Pn
Model p0 = p + 1 i=1 Ŷi
2
2
Correction Factor 1 nY
Pn P 2
Regression p i=1 (Ŷi − Y = ni=1 Ŷi2 − nY
)2
Pn 2
Residual n-p’ i=1 (Yi − Ŷi )
Table 5: The Analysis of Variance

.
The mean squares for regression and residuals are the corresponding sums of squares divided
by their respective freedoms:
SS(REGRESSION)
M S(REGRESSION) =
p
SS(RESIDUAL)
M S(RESIDUAL) =
n − p0
The expected value of the mean squares are given below (for the case where there is a single
independent variable), the proof is given later. These are based on the assumption that the model
is fit is the correct model.
n
X
E[M S(REGRESSION)] = σ 2 + β12 (Xi − X)2
i=1
E[M S(RESIDUAL)] = σ2
The Coefficient of Determination (R2 ) is the ratio of the regression sum of squares to the
total (corrected) sum of squares, and represents the fraction of the variation in Y that is “explained”
by the set of independent variables X1 , . . . , Xp .
SS(REGRESSION) SS(RESIDUAL)
R2 = = 1−
SS(TOTAL CORRECTED) SS(TOTAL CORRECTED)
1.3.1 Examples
Numerical results for the two examples desribed before are given below.
The Analysis of Variance for the LSD/math score data are given in Table 1.3.1. Here, n = 7,
p = 1, and p0 = 2. All relevant sums are obtained from previous examples.

Pn 2
Total (Uncorrected) n=7 i=1 Yi = 19639.24
2
Correction Factor 1 nY = 17561.02
Pn 2 Pn 2 2
Total (Corrected) n−1=6 i=1 (Yi − Y ) = i=1 Yi − nY = 2078.22
Pn
Model p0 = p + 1 = 2 i=1Ŷi2 = 19385.32
2
Pn P 2
Regression p=1 i=1 (Ŷi −PY ) = ni=1 Ŷi2 − nY = 1824.30
2
n 2
Residual n-p’=5 i=1 (Yi − Ŷi ) = 253.92
Table 6: The Analysis of Variance for the LSD/math score data

.
The mean squares for regression and residuals are as follow:
SS(REGRESSION) 1824.30
M S(REGRESSION) = = = 1824.30
p 1
SS(RESIDUAL) 253.92
M S(RESIDUAL) = = = 50.78
n − p0 5
The coefficient of determination for this data is:
R2 = = = 0.8778
SS(TOTAL CORRECTED) 2078.22
Approximately 88% of the variation in math scores is “explained” by the linear relation between
math scores and LSD concentration.
The Analysis of Variance for the output/cost data are given in Table 1.3.1. Here, n = 48, p = 1,
and p0 = 2. All relevant sums are obtained from previous examples and computer output.

Pn 2
Total (Uncorrected) n = 48 i=1 Yi = 238424.46
2
Pn 2 Pn 2 2
Total (Corrected) n − 1 = 47 i=1 (Yi − Y ) = i=1 Yi − nY = 32914.17
Pn
Model p0 = p + 1 = 2 i=1 Ŷi
2
= 236636.27
2
Pn P 2
Regression p=1 i=1 (Ŷi −PY )2 = ni=1 Ŷi2 − nY = 31125.98
n 2
Residual n-p’=46 i=1 (Yi − Ŷi ) = 1788.19
Table 7: The Analysis of Variance for the hosiery mill cost function data
.
The mean squares for regression and residuals are as follow:
M S(REGRESSION) = = = 31125.98
p 1
SS(RESIDUAL) 1788.19
M S(RESIDUAL) = = = 38.87
n − p0 46
The coefficient of determination for this data is:
R2 = = = 0.9457
SS(TOTAL CORRECTED) 32914.17
Approximately 95% of the variation in math scores is “explained” by the linear relation between
math scores and LSD concentration.
1.4 Precision and Distribution of Estimates

Important results from mathematical statistics regarding linear functions of random variables. Let
P
U = ni=1 ai Yi , where a1 , . . . , an are fixed constants and Yi are random variables with E(Yi ) = µi ,
V ar(Yi ) = σ 2 , and Cov(Yi , Yj ) = 0, i 6= j:
X
n X
n X
n
E[U ] = E[ ai Yi ] = ai E[Yi ] = ai µi (13)
i=1 i=1 i=1
n
X n
X n
X
V ar[U ] = V ar[ ai Yi ] = a2i V ar[Yi ] = ai σi2 (14)
i=1 i=1 i=1
Pn n
Y h ∗ i
E[etU ] = E[et i=1
a i Yi
] = E [eta1 Y1 · · · etan Yn = E eti Yi t∗i = ai t (15)
i=1
1.4.1 Distribution of β̂1
Consider β̂1 as a linear function of Y1 , . . . , Yn :
Pn Pn P n n
i=1 (Xi − X)(Yi − Y) i=1 (Xi − X)Yi − Y ni=1 (Xi − X) X Xi − X X
β̂1 = Pn 2
= Pn 2
= Pn 2
Y i = ai Yi
i=1 (Xi − X) i=1 (Xi − X) i=1 i=1 (Xi − X) i=1
Under the simple linear regression model:
E[Yi ] = µi = β 0 + β 1 Xi V ar[Yi ] = σi2 = σ2
From equation (13):

n
X n
X Xi − X
E[β̂1 ] = ai E[Yi ] = Pn 2
(β0 + β1 Xi )
i=1 i=1 i=1 (Xi − X)
Xn Xn n
X
β0 β1 β1
= Pn 2
(X i − X) + Pn 2
(X i − X)X i = Pn 2
(Xi − X)Xi
i=1 (Xi − X) i=1 i=1 (Xi − X) i=1 i=1 (Xi − X) i=1
" n n
# " n #
β1 X X β 1
X 2
= Pn 2
Xi2 − X Xi = Pn 2
Xi2 − nX
i=1 (Xi − X) i=1 i=1 i=1 (X i − X) i=1
n
X
β1
= Pn 2
(Xi − X)2 = β1
i=1 (Xi − X) i=1
From Equation (14):

n
X n
X n
X 2
Xi − X
V ar[β̂1 ] = a2i V ar[Yi ] = a2i σ 2 = σ2 Pn 2
i=1 i=1 i=1 i=1 (Xi − X)
2 X
n
1 σ2
= Pn 2
σ2 (Xi − X)2 = Pn 2
i=1 (X i − X) i=1 i=1 (Xi − X)
With the further assumption that Yi (and, more specifically, εi ) being normally distributed, we
can obtain the specific distribution of β̂1 .
 P P 
h i n
t
n
Pn Xi −X Yi h Pn ∗ i
i=1 i=1 (Xi −X)2
E etβ̂1 = E e i=1  = E e i=1 ti Yi
where t∗i = t Pn X(X

i −X
−X)2
. If Y ∼ N (µ, σ 2 ), then the moment generating function for Y is:
i=1 i
σ 2 t2
mY (t) = E etY = eµt+ 2
( 2 )
h ∗ i Xi − X Xi − X σ 2 t2
⇒ E eti Yi = exp (β0 + β1 Xi ) Pn t+ Pn
i=1 (Xi − X)
2
i=1 (Xi − X)
2 2
By independence of the Yi , we get:

 P 
h i n
t Pn Xi −X Yi n
Y h ∗ i
i=1 (Xi −X)2
E etβ̂1 = E e i=1  = E eti Yi
i=1
n
Y 2
Xi − X Xi − X σ 2 t2
exp{(β0 + β1 Xi ) Pn 2
t+ Pn 2
}
i=1 i=1 (Xi − X) i=1 (Xi − X)
2
( n n 2 )
X Xi − X X Xi − X σ 2 t2
= exp (β0 + β1 Xi ) Pn t+ Pn (16)
i=1 i=1 (Xi − X)
2
i=1 i=1 (Xi − X)
2 2
Expanding the first term in the exponent in equation (16), we get:

n
X Xn Xn
Xi − X t
(β0 + β1 Xi ) Pn 2
t = Pn 2
{β 0 (X i − X) + β 1 Xi (Xi − X)}
i=1 i=1 (Xi − X) i=1 (Xi − X) i=1 i=1
Xn
t
= Pn 2
{0 + β 1 (Xi − X)2 } = β1 t (17)
i=1 (X i − X) i=1
Expanding the second term in the exponent in equation (16), we get:
n
X 2 Xn
Xi − X σ 2 t2 σ 2 t2 σ 2 t2
Pn 2
= Pn (Xi − X)2 = Pn (18)
i=1 i=1 (Xi − X)
2 2( i=1 (X i − X) 2 )2
i=1
2 i=1 (Xi − X)
2
Putting equations (17) and (18) back into equation 1(16), we get:
h i σ 2 t2
mβ̂1 (t) = E etβ̂1 = exp{β1 t + Pn 2
2 i=1 (Xi − X)
which is the moment generating function of a normally distributed random variable with mean
P
β1 and variance σ 2 / ni=1 (Xi − X)2 . Thus, we have the complete sampling distribution of β̂1 under
the model’s assumptions.
1.4.2 Distribution of β̂0

Consider the following results from mathematical statistics:
n
X n
X
U = ai Yi W = di Yi
i=1 i=1
where {ai } and {di } are constants and {Yi } are random variables. Then:
n
X n X
X
Cov[U, W ] = ai di V (Yi ) + ai dj Cov[Yi , Yj ]
i=1 i=1 j6=i
Then, we can write β̂0 as two linear functions of the {Yi }:

Xn
1 Xn
X(Xi − X)
β̂0 = Y − β̂1 X = Yi − Pn Y
2 i
= U −W (19)
i=1
n i=1 i=1 (Xi − X)
The expected values of the the two linear functions of the Yi in equation (19) are as follow:
Xn n
1 1 1 X
E[U ] = (β0 + β1 Xi ) = (nβ0 ) + β1 Xi = β0 + β1 X
i=1
n n n i=1
n
X X(Xi − X)
E[W ] = Pn 2
(β0 + β1 Xi )
i=1 i=1 (Xi − X)
Xn
X
= Pn 2
[β0 Xi + β1 Xi2 − β0 X − β1 XXi ]
i=1 (X i − X) i=1
Xn Xn
X 2
= Pn [β
2 0
(X i − X) + β 1 ( Xi2 − nX )]
i=1 (X i − X) i=1 i=1
Xn
X
Pn 2
[0 + β 1 (Xi − X)2 ] = β1 X
i=1 (X i − X) i=1
Putting these together in equation (19):
E[β̂0 ] = E[U − W ] = E[U ] − E[W ] = β0 + β1 X − β1 X = β0
Now to get the variance of β̂0 (again assuming that Cov[Yi , Yj ] = 0f ori 6= j) :
V ar[U − W ] = V ar[U ] + V ar[W ] − 2Cov[U, W ]
Xn n 2
X 2
1 1 1 σ2
V ar[U ] = V ar[ Yi ] = V ar[Yi ] = n σ2 =
i=1
n i=1
n n n
2 2 σ2
V ar[W ] = V ar[β̂1 X] = X V ar[β̂1 ] = X Pn 2
i=1 (Xi − X)
Xn n
X
1 X(Xi − X) σ2 X
Cov[U, W ] = Pn V ar[Yi ] = Pn (Xi − X) = 0
i=1
n i=1 (Xi − X)
2 n 2
i=1 (Xi − X) i=1
2
" 2
#
σ2 X σ2 2 1 X
⇒ V ar[β̂0 ] = V ar[U ]+V ar[W ]−2Cov[U, W ] = + Pn = σ + Pn
n i=1 (Xi − X)
2 n i=1 (Xi − X)
2
Note that Cov[U, W ] = Cov[Y , β̂1 X] = 0, then Y and β̂1 are independent. We can also write
β̂0 as a single linear function of Yi , allowing use of the moment generating function method to
determine it’s normal sampling distribution:
" 2
#
Xn
1 X (Xi − X) X
n
β̂0 = − Pn Yi = ai Yi
i=1
n i=1 (Xi − X)
2
i=1
1.4.3 Distribution of Ŷi
The distibution of Ŷi , which is an estimate of the population mean of Yi at the level Xi of the
independent variable is obtained as follows:
Ŷi = β̂0 + β̂1 Xi = (Y − β̂1 X) + β̂1 Xi = Y + β̂1 (Xi − X)
E[Ŷi ] = E[Y ] + (Xi − X)E[β̂1 ] = β0 + β1 X + β1 (Xi − X) = β0 + β1 Xi
σ2 (Xi − X)2 σ 2
V ar[Ŷi ] = V ar[Y ] + (Xi − X)2 V ar[β̂1 ] + 2(Xi − X)Cov[Y , β̂1 ] = + Pn +0
n i=1 (Xi − X)
2

2 1 (Xi − X)2
= σ + Pn
n i=1 (Xi − X)
2
Further, since Ŷi is a linear function of the Yi , then Ŷi has a normal sampling distribution.
1.4.4 Prediction of future observation Y0 when X = X0

The predicted value of a future observation Y0 is Ŷpred0 = β̂0 + β̂1 X0 . The prediction error is Y0 − Ŷ0 ,
and the quantity E[(Ŷ0 − Y0 )2 ] is referred to as the mean square error of prediction. Assuming the
model is correct:
E[Ŷ0 − Y0 ] = 0

1 (X0 − X)2 1 (X0 − X)2
V ar[Ŷpred0 ] = V ar[Ŷ0 −Y0 ] = σ2 + Pn +σ 2 = σ 2 1 + + Pn
n i=1 (X i − X) 2 n i=1 (Xi − X)
2
1.4.5 Estimated Variances

For all of the sampling distributions derived above, the unknown observation variance σ 2 appears
in the estimators’ variances. To obtain the estimated variance for each of the estimators, we replace
σ 2 with s2 = M S(RESIDUAL). It’s important to keep in mind that this estimator is unbiased for
σ 2 only if the model is correctly specified. The estimated variances for each estimator and predictor
are given below:
s2
• s2 (β̂1 ) = Pn (Xi −X)2
Estimated variance of β̂1
i=1

2
1 X
• s2 (β̂ 0) = s2 n + Pn (X 2
Estimated variance of β̂0
i=1 i −X)

2
• s2 (Ŷ i) = s2 1
+ Pn(Xi −X) 2 Estimated variance of estimated mean at Xi
n (X i −X)
i=1

2
• s2 (Ŷ pred0 ) = s2 1+ 1
n + Pn(X0 −X)(Xi −X)2
Estimated variance of prediction at X0
i=1
1.4.6 Examples
Estimated variances are computed for both of the previous examples.
Here we obtain estimated variances for β̂1 , β̂0 , the true mean, and a future score when the tissue
concentration is 5.0:
s2 50.78
s2 (β̂1 ) = Pn 2
= = 2.26
i=1 (Xi − X)
22.48
" 2
#
2 2 1 X 1 4.33292
s (β̂0 ) = s + Pn = 50.78 + = 49.66
n i=1 (Xi − X)
2 7 22.48

1 (Xi − X)2 1 (5 − 4.3329)2
s2 (Ŷ5 ) = s2 + Pn = 50.78 + = 8.26
n i=1 (Xi − X)
2 7 22.48

2 2 1 (X0 − X)2 1 (5 − 4.3329)2
s (Ŷpred0 ) = s 1 + + Pn = 50.78 1 + + = 59.04
n i=1 (Xi − X)
2 7 22.48
Here we obtain estimated variances for β̂1 , β̂0 , the true mean, and a future cost when the
production output is 30.0:
s2 38.87
s2 (β̂1 ) = Pn 2
= = 0.0050
i=1 (Xi − X)
7738.94
" 2
#
2 2 1 X 1 31.06732
s (β̂0 ) = s + Pn = 38.87 + = 5.66
n i=1 (Xi − X)
2 48 7738.94

1 (Xi − X)2 1 (30 − 31.0673)2
s2 (Ŷ30 ) = s2 + Pn = 38.87 + = 0.82
n i=1 (Xi − X)
2 48 7738.94

2 2 1 (X0 − X)2 1 (30 − 31.0673)2
s (Ŷpred0 ) = s 1 + + Pn = 38.87 1 + + = 39.69
n i=1 (Xi − X)
2 48 7738.94
1.5 Tests of Significance and Confidence Intervals

Under the model assumptions of independence, normality and constant error variance; we can make
inferences concerning model parameters. We can conduct t–tests, F –tests, and obtain confidence
intervals regarding the unknown parameters.
1.5.1 Tests of Significance
The t-test can be used to test hypotheses regarding β0 or β1 , and can be used for 1–sided or 2–sided
alternative hypotheses. The form of the test is as follows, and can be conducted regarding any of
the regression coefficients:
• H0 : βi = m (m specified, usually 0 when testing β1 )
• (1) Ha : βi 6= m
(2) Ha : β1 > m
(3) Ha : β1 < m
β̂i −m
• T S : t0 =
s(β̂i )
• (1) RR : |t0 | ≥ t(α/2,n−p0 ) (p0 = 2 for simple regression)

(2) RR : t0 ≥ t(α,n−p0 ) (p0 = 2 for simple regression)
(3) RR : t0 ≤ −t(α,n−p0 ) (p0 = 2 for simple regression)
• (1) P –value: 2 · P (t ≥ |t0 |)
(2) P –value: P (t ≥ t0 )
(3) P –value: P (t ≤ t0 )
Using tables, we can only place bounds on these p–values, but statistical computing packages
will print them directly.
A second test is available to test whether the slope parameter is 0 (no linear association exists
between Y and X). This is based on the Analysis of Variance and the F -distribution:
1. H0 : β1 = 0 HA : β1 6= 0 (This will always be a 2–sided test)
M S(REGRESSION)
2. T.S.: F0 = M S(RESIDUAL)
3. R.R.: F0 > F(α,1,n−p0 (p0 = 2 for simple regression)

4. p-value: P (F > F0 ) (You can only get bounds on this from tables, but computer outputs
report them exactly)
Under the null hypothesis, the test statistic should be near 1, as β1 moves away from 0, the test
statistic should increase.
1.5.2 Confidence Intervals

Confidence intervals for model parameters can be obtained under all the previously stated assump-
tions. The 100(1 − α)100% confidence intervals can be obtained as follows:
β0 : β̂0 ± t(α/2,n−p0 ) s(β̂0 )
β1 : β̂1 ± t(α/2,n−p0 ) s(β̂1 )
β0 + β1 Xi : Ŷi ± t(α/2,n−p0 ) s(Ŷi )

Prediction intervals for future observations at X = X0 can be obtained as well in an obvious
manner.
1.5.3 Examples
The previously described examples are continued here.
To determine whether there is a negative association between math scores and LSD q concentra-
tion, we conduct the following test at α = 0.05 significance level. Note that s(β̂1 ) = s2 (β̂1 ) =
√
2.26 = 1.50.
H0 : β1 = 0 Ha : β1 < 0
β̂1 − 0 −9.01
T S : t0 = = = −6.01 RR : t0 ≤ −t.05,5 = −2.015 P − val = P (t ≤ −6.01)
s(β̂1 ) 1.50
Next, we obtain a confidence interval for the true mean

q score when the tissue concentration is
√
X = 5.0. The estimated standard error of Ŷ5 is s(Ŷ5 ) = s2 (Ŷ5 ) = 8.26 = 2.87, and t(0.025,5) =
2.571. The 95% confidence interval for β0 + β1 (5) is:
Ŷ5 = 89.12 − 9.01(5) = 44.07 44.07 ± 2.571(2.87) ≡ 44.07 ± 7.38 ≡ (36.69, 51.45)
Here, we use the F -test to determine whether there is an association between product costs and
the production output at α = 0.05 significance level.
H0 : β1 = 0 Ha : β1 6= 0
M S(REGRESSION) 31125.98
T S : F0 = = = 800.77 RR : F0 ≥ F(.05,1,46) ≈ 1.680 P −val = P (F ≥ 800.77)
M S(RESIDUAL) 38.87
Unit variable cost is the average increment in total production cost per unit increase in produc-
tion output (β1 ). We obtain a 95% confidence interval for this parameter:
√
β̂1 = 2.0055 s2 (β̂1 ) = .0050 s(β̂1 ) = .0050 = .0707 t(.025,46) ≈ 2.015
β̂1 ± t(.025,46) s(β̂1 ) 2.0055 ± 2.015(.0707) 2.0055 ± 0.1425 (1.8630, 2.1480)
As the production output increases by 1000 dozen pairs, we are very confident that mean costs
increase by between 1.86 and 2.15 $1000. The large sample size n = 48 makes our estimate very
precise.
1.6 Regression Through the Origin
In some practical situations, the regression line is expected (theoretically) to pass trough the origin.
It is important that X = 0 is a reasonable level of X in practice for this to be the case. For instance,
in the hosiery mill example, if firm knows in advance that production will be 0 they close plant
and have no costs if they are able to work in “short–run,” however most firms still have “long–run”
costs if they know they will produce in future. If a theory does imply that the mean response (Y )
is 0 when X = 0, we have a new model:
Yi = β 1 Xi + ε i i = 1, . . . , n
The least squares estimates are obtained by minimizing (over β1 ):

n
X
Q = (Yi − β1 Xi )2
i=1
This is obtained by taking the derivative of Q with respect to β1 , and setting it equal to 0. The
value β̂1 that solves that equality is the least squares estimate of β1 .
n
X
∂Q
= 2 (Yi − β1 Xi )(−Xi ) = 0
∂β1 i=1
n n Pn
X X X i Yi
⇒ Yi Xi = β̂1 Xi2 ⇒ β̂1 = Pi=1
n 2
i=1 i=1 i=1 Xi
The estimated regression equation and residuals are:
Ŷi = β̂1 Xi ei = Y − Ŷi = Y − β̂1 Xi
Note that for this model, the residuals do not necessarily sum to 0:
Pn
i=1 Xi Y i
ei = Y − Ŷi = Y − β̂1 Xi = Yi − Pn 2 Xi
i=1 Xi
n
X n
X Pn n
Xi Y i X
⇒ ei = Yi − Pi=1
n 2 Xi
i=1 i=1 i=1 Xi i=1
This last term is not necessarily (and will probably rarely, if ever, in practice be) 0.
The uncorrected sum of squares is:
n
X n
X n
X n
X
Yi2 = Ŷi2 + e2i + 2 Ŷi ei
i=1 i=1 i=1 i=1
The last term (the cross–product term) is still 0 under the no–interecept model:
n
X n
X Pn Pn
Pi=1 Xi Yi Pi=1 Xi Yi
ei Ŷi = (Yi − n 2 Xi )( n 2 Xi )
i=1 i=1 i=1 Xi i=1 Xi
X Pn Xn Pn Pn n Pn 2 X
Xi Y i X
n n
i=1 Xi Yi i=1 Xi Yi i=1 Xi Yi
= Yi ( P n 2 Xi )− ( P n 2 Xi ) 2 = Pi=1
n 2 Xi Y i − P n 2 Xi2 = 0
i=1 i=1 Xi i=1 i=1 Xi i=1 Xi i=1 i=1 Xi i=1
So, we obtain the same partitioning of the total sum of squares as before:
n
X n
X n
X
Yi2 = Ŷi2 + e2i
i=1 i=1 i=1
SS(TOTAL UNCORRECTED) = SS(MODEL) + SS(RESIDUAL)
The model sum of squares is based on only one parameter, so it is not broken into the components
of mean and regression as it was before. Similarly, the residual sum of squares has n − 1 degrees of
freedom. Assuming the model is correct:
n
X
E[M S(MODEL)] = σ 2 + β1 Xi2
i=1
E[M S(RESIDUAL)] = σ2
The variance of the estimator β̂1 is:
"P # n
n
X i Yi 1 X
V ar[β̂1 ] = V ar Pi=1
n 2 = Pn 2 2 Xi2 V ar[Yi ]
i=1 Xi ( i=1 Xi ) i=1
n
X
1 σ2
= Pn 2 2
( Xi2 )σ 2 = Pn 2
( i=1 Xi ) i=1 i=1 Xi
Similarly, the variance of Ŷ0 = X0 β̂1 is:
σ 2 X02
V ar[Ŷ0 ] = V ar[X0 β̂1 ] = X02 V ar[β̂1 ] = Pn 2
i=1 Xi
Estimates are obtained by replacing σ 2 with s2 = M S(RESIDUAL).
1.6.1 Example – Galton’s Height Measurements

In what is considered by many to be the first application of regression analysis, Sir Frances Galton
(1889, Natural Inheritance) obtained heights of n = 928 adult children (Y ) and the“midheight” of
their parents (X). Since the mean heights of adult children and their parents were approximately
the same (68.1” for adult children and 68.3” for their parents). Once both datasets have been
centered around their means, Galton found that adult chidrens heights were less extreme than
their parents. This phenomenon has been observed in many areas of science, and is referred to as
regression to the mean.
Here we fit a regression model through the origin, which for this centered data is the point
(68.1,68.3). We have the following quantities based on the centered data given in Galton’s table:
n
X n
X n
X
n = 928 Xi2 = 3044.92 Yi2 = 5992.48 Xi Yi = 1965.46
i=1 i=1 i=1
From this data, we obtain the following quantities:

Pn
X i Yi 1965.46
β̂1 = β̂1 = Pi=1
n 2 = = 0.6455
i=1 Xi 3044.92
n
X n
X n
X
SS(MODEL) = Ŷi2 = (β̂1 Xi )2 = β̂12 Xi2 = (0.6455)2 (3044.92) = 1268.73
i=1 i=1 i=1
n
X n
X
SS(RESIDUAL) = Yi2 − Ŷi2 = 5992.48 − 1268.73 = 4723.75
i=1 i=1
SS(RESIDUAL) 4723.75
s2 = M S(RESIDUAL) = = = 5.10
n−1 927
s2 5.10
s2 (β̂1 ) = Pn 2 = = .0017 s(β̂1 ) = .0409
i=1 Xi 3044.92
From this, we get a 95% confidence interval for β1 :
bh1 ± z(.025) s(β̂1 ) ≡ 0.6455 ± 1.96(.0409) ≡ 0.6455 ± 0.0802 ≡ (0.5653, 0.7257)
Note that there is a positive association between adult children’s heights and their parent’s
heights. However, as the parent’s height increases by 1”, the adult child’s height increases by
between 0.5633” and 0.7257” on average. This is an example of regression to the mean.
1.7 Models with Several Independent Variables

As was discussed in the section on the Analysis of Variance, models can be generalized to contain
p < n independent variables. However, the math to obtain estimates and their estimated variances
and standard errors is quite messier. This can be avoided by making use of matrix algebra, which
is introduced shortly. The general form of the multiple linear regression model is:
Yi = β0 + β1 Xi1 + · · · + βp Xip + εi ε ∼ N ID(0, σ 2 )
The least squares estimates β̂0 , β̂1 , . . . , β̂p are the values that minimize the residual sum of
squares:
n
X n
X
SS(RESIDUAL) = (Yi − Ŷi )2 = (Yi − β̂0 + β̂1 Xi1 + · · · + β̂p Xip )2
i=1 i=1
An unbiased estimate of σ 2 is:

SS(RESIDUAL)
s2 = SS(RESIDUAL) =
n − (p + 1)
We will obtain these estimates after we write the model in matrix notation.
1.8 SAS Programs and Output

In this section, SAS code and its corresponding output are given for the two examples in Rawlings,
Pantula, and Dickey (RPD).
2 Introduction to Matrices
Text: RPD, Sections 2.1-2.6
Problems:
In this section, important definitions and results from matrix algebra that are useful in regression
analysis are introduced. While all statements below regarding the columns of matrices can also be
said of rows, in regression applications we will typically be focusing on the columns.
A matrix is a rectangular array of numbers. The order or dimension of the matrix is the
number of rows and columns that make up the matrix. The rank of a matrix is the number of
linearly independent columns (or rows) in the matrix.
A subset of columns is said to be linearly independent if no column in the subset can
be written as a linear combination of the other columns in the subset. A matrix is full rank
(nonsingular) if there are no linear dependencies among its columns. The matrix is singular if
lineardependencies exist.
The column space of a matrix is the collection of all linear combinations of the columns of a
matrix.
The following are important types of matrices in regression:
Vector – Matrix with one row or column
Square Matrix – Matrix where number of rows equals number of columns
Diagonal Matrix – Square matrix where all elements off main diagonal are 0
Identity Matrix – Diagonal matrix with 1’s everywhere on main diagonal
Symmetric Matrix – Matrix where element aij = aji ∀i, j
Scalar – Matrix with one row and one column (single element)
The transpose of a matrix is the matrix generated by interchanging the rows and columns of
the matrix. If the original matrix is A, then its transpose is labelled A0 . For example:
 
" # 2 1
2 4 7  
A= ⇒ A0 =  4 7 
1 7 2
7 2
Matrix addition (subtraction) can be performed on two matrices as long as they are of
equal order (dimension). The new matrix is obtained by elementwise addition (subtraction) of the
two matrices. For example:
" # " # " #
2 4 7 1 3 0 3 7 7
A= B= ⇒ A+B=
1 7 2 2 4 8 3 11 10
Matrix multiplication can be performed on two matrices as long as the number of columns
of the first matrix equals the number of rows of the second matrix. The resulting has the same
number of rows as the first matrix and the same number of columns as the second matrix. If
C = AB and A has s columns and B has s rows, the element in the ith row and j th column of C,
which we denote cij is obtained as follows (with similar definitions for aij and bij ):
s
X
cij = ai1 b1j + ai2 b2j + · · · ais bsj = aik bkj
k=1
For example:  
" # 1 5 6
2 4 7  
A= B= 2 0 1  ⇒
1 7 2
3 3 3
" # " #
2(1) + 4(2) + 7(3) 2(5) + 4(0) + 7(3) 2(6) + 4(1) + 7(3) 31 31 37
C = AB = =
1(1) + 7(2) + 2(3) 1(5) + 7(0) + 2(3) 1(6) + 7(1) + 2(3) 21 11 19
Note that C has the same number of rows as A and the same number of columns as C. Note
that in general AB 6= BA; in fact, the second matrix may not exist due to dimensions of matrices.
However, the following equality does hold: (AB)0 = B0 A0 .
Scalar Multiplication can be performed between any scalar and any matrix. Each element
of the matrix is multiplied by the scalar. For example:
" # " #
2 4 7 4 8 14
A= ⇒ 2A =
1 7 2 2 14 4
The determinant is scalar computed from the elements of a matrix via well–defined (although
rather painful) rules. Determinants only exist for square matrices. The determinant of a matrix A
is denoted as |A|.
For a scalar (a 1 × 1 matrix): |A| = A.

For a 2 × 2 matrix: |A| = a11 a22 − a12 a21 .
For n × n matrices (n > 2):
1. Ars ≡ (n − 1) × (n − 1) matrix with row r and column s removed from A
2. |Ars | ≡ the minor of element ars
3. θrs = (−1)r+s |Ars | ≡ the cofactor of element ars
4. The determinant is obtained by summing the product of the elements and cofactors for any
P
row or column of A. By using row i of A, we get |A| = nj=1 aij θij
Example – Determinant of a 3 × 3 matrix

We compute the determinant of a 3 × 3 matrix, making use of its first row.
 
10 5 2
 
A =  6 8 0 
2 5 1
" #
8 0
a11 = 10 A11 = |A11 | = 8(1) − 0(5) = 8 θ11 = (−1)1+1 (8) = 8
5 1
" #
6 0
a12 = 5 A12 = |A11 | = 6(1) − 0(2) = 6 θ12 = (−1)1+2 (6) = −6
2 1
" #
6 8
a13 = 2 A13 = |A13 | = 6(5) − 8(2) = 14 θ13 = (−1)1+3 (14) = 14
2 5
Then the determinant of A is:

n
X
|A| = a1j θ1j = 10(8) + 5(−6) + 2(14) = 78
j=1
Note that we would have computed 78 regardless of which row and column we used.
An important result in linear algebra states that if |A| = 0, then A is singular, otherwise A is
nonsingular (full rank).
The inverse of a square matrix A, denoted A−1 , is a matrix such that A−1 A = I = AA−1
where I is the identity matrix of the same dimension as A. A unique inverse exists if A is square
and full rank.
The identity matrix, when multiplied by any matrix (such that matrix multiplication exists)
returns the same matrix. That is: AI = A and IA = A, as long as the dimensions of the matrices
conform to matrix multiplication.
For a scalar (a 1 × 1 matrix): A−1

" =1/A. #
−1 1 a22 −a12
For a 2 × 2 matrix: A = |A| .
−a21 a11
For n × n matrices (n > 2):
1. Replace each element with its cofactor (θrs )
2. Transpose the resulting matrix
3. Divide each element by the determinant of the original matrix
Example – Inverse of a 3 × 3 matrix

We compute the inverse of a 3 × 3 matrix (the same matrix as before).
 
10 5 2
 
A =  6 8 0  |A| = 78
2 5 1
|A11 | = 8 |A12 | = 6 |A13 | = 14

|A21 | = −5 |A22 | = 6 |A23 | = 40
|A31 | = −16 |A32 | = −12 |A33 | = 50
θ11 = 8 θ12 = −6 θ13 = 14 θ21 = 5 θ22 = 6 θ23 = −40 θ31 = −16 θ32 = 12 θ33 = 50
 
8 5 −16
1  
A−1 =  −6 6 12 
78
14 −40 50
As a check:
    
8 5 −16 10 5 2 1 0 0
1     
A−1 A =  −6 6 12   6 8 0  =  0 1 0  = I3
78
14 −40 50 2 5 1 0 0 1
To obtain the inverse of a diagonal matrix, simply compute the recipocal of each diagonal
element.
The following results are very useful for matrices A, B, C and scalar λ, as long as the matrices’
dimensions are conformable to the operations in use:
1. A + B = B+A
2. (A + B) + C = A + (B + C)
3. (AB)C = A(BC)
4. C(A + B) = CA + CB
5. λ(A + B) = λA + λB
6. (A0 )0 = A
7. (A + B)0 = A 0 + B0
8. (AB)0 = B 0 A0
9. (ABC)0 = C 0 B0 A0
10. (AB)−1 = B−1 A−1
11. (ABC)−1 = C−1 B−1 A−1
12. (A−1 )−1 = A
13. (A0 )−1 = (A−1 )0

The length of a column vector x and the distance between two column vectors u and v are:
√ q
l(x) = x0 x l((u − v)) = (u − v)0 (u − v)
Vectors x and w are orthogonal if x0 w = 0.

2.1 Linear Equations and Solutions
Suppose we have a system of r linear equations in s unknown variables. We can write this in matrix
notation as:
Ax = y
where x is a s × 1 vector of s unknowns; A is a r × s matrix of known coefficients of the s

unknowns; and y is a r × 1 vector of known constants on the right hand sides of the equations.
This set of equations may have:
• No solution
• A unique solution
• An infinite number of solutions
A set of linear equations is consistent if any linear dependencies among rows of A also appear
in the rows of y. For example, the following system is inconsistent:
    
1 2 3 x1 6
    
 2 4 6   x2  =  10 
3 3 3 x3 9
This is inconsistent because the coefficients in the second row of A are twice those in the first row,
but the element in the second row of y is not twice the element in the first row. There will be no
solution to this system of equations.
A set of equations is consistent if r(A) = r([Ay]) where [Ay] is the augmented matrix [A|y].
When r(A) equals the number of unknowns, and A is square:
x = A−1 y
2.2 Projection Matrices

The goal of regression is to transform a n-dimensional column vector Y onto a vector Ŷ in a
subspace (such as a straight line in 2-dimensional space) such that Ŷ is as close to Y as possible.
Linear transformation of Y to Ŷ, Ŷ = PY is said to be a projection iff P is idempotent and
symmetric, in which case P is said to be a projection matrix.
A square matrix A is idempotent if AA = A. If A is idempotent, then:

n
X
r(A) = aii = tr(A)
i=1
where tr(A) is the trace of A. The subspace of a projection is defined, or spanned, by the columns
or rows of the projection matrix P.
Ŷ = PY is the vector in the subspace spanned by P that is closest to Y in distance. That is:
SS(RESIDUAL) = (Y − Ŷ)0 (Y − Ŷ)

is at a minimum. Further:
e = (I − P)Y
is a projection onto a subspace orthogonal to the subspace defined by P.
Ŷ 0 e = (PY)0 (I − P)Y = Y 0 P0 (I − P)Y = Y 0 P(I − P)Y = Y 0 (P − P)Y = 0
Ŷ + e = PY + (I − P)Y = Y
2.3 Vector Differentiation

Let f be a function of x = [x1 , . . . , xp ]0 . We define:
 
 ∂f 
 ∂x1 
df  ∂f 
 
=  ∂x2 
dx  .. 
 . 
 
∂f
∂xp
From this, we get for p × 1 vector a and p × p symmetric matrix A:
d(a0 x) d(x0 Ax)

=a = 2Ax
dx dx
“Proof ” – Consider p = 3:
d(a0 x) d(a0 x)
a0 x = a1 x1 + a2 x2 + a3 x3 = ai ⇒ =a
dxi dx
 
x1
x0 Ax = x1 a11 + x2 a21 + x3 a31 x1 a12 + x2 a22 + x3 a32 x1 a13 + x2 a23 + x3 a33  x2  =
x3
= x21 a11 + x1 x2 a21 + x1 x3 a31 + x1 x2 a12 + x22 a22 + x2 x3 a32 + x1 x3 a13 + x2 x3 a23 + x23 a33
∂x0 Ax X
⇒ = 2aii xi + 2 aij xj (aij = aji )
∂xi
j6=i
 0   
∂x Ax
∂x1 2a11 x1 + 2a12 x2 + 2a13 x3
∂x0 Ax  ∂x0 Ax  
⇒ = ∂x2  = 2a21 x1 + 2a22 x2 + 2a23 x3  = 2Ax
∂x ∂x0 Ax 2a31 x1 + 2a32 x2 + 2a33 x3
∂x3
2.4 SAS Programs and Output

In this section, SAS code and its corresponding output are given for the examples in Rawlings,
Pantula, and Dickey (RPD).
3 Multiple Regression in Matrix Notation
Multiple linear regression model:
Yi = β0 + β1 Xi1 + β2 Xi2 + · · · βp Xip + εi i = 1, . . . , n
where i represents the observational unit, the second subscript on X represents the independent
variable number, p is the number of independent variables, and p0 = p + 1 is the number of
model parameters (including the intercept term). For the model to have a unique set of regression
coefficients, n > p0 .
We can re–formulate the model in matrix notation:
Y — n × 1 column vector of observations on the dependent variable Y
X — n × p0 model matrix containing a column of 1’s and p columns of levels of the independent
variables X1 , . . . , Xp
β — p0 × 1 column vector of regression coefficients (parameters)
ε — n × 1 column vector of random errors
Y = X β+ ε
      
Y1 1 X11 ··· X1p β0 ε1
 Y2   1 X21 ··· X2p  β1   ε2 
      
 ..  =  .. .. .. ..  .. + .. 
 .   . . . .  .   . 
Yn 1 Xn1 ··· Xnp βp εn
For our models, X will be of full column rank, meaning r(X) = p0 .
The elements of β, are referred to as partial regression coefficients, βj represents the change
in E(Y ) as the j th independent variable is increased by 1 unit, while all other variables are held
constant. The terms “controlling for all other variables” and “ceteris parabis” are also used to
describe the effect.
We will be working with many different models (that is, many different sets of independent
variables). Often we will need to be more specific of which independent variables are in our model.
We denote the partial regression coefficient for X2 in a model containing X1 , X2 , and X3 as β2·13 .
3.1 Distributional Properties

We still have the same assumption on the error terms as before:
εi ∼ N ID(0, σ 2 ) i = 1, . . . , n
This implies that the joint probability density function for the random errors is:
n n
" ( )# ( Pn )
Y Y −ε2i − 2
i=1 εi
f (ε1 , ε2 , . . . , εn ) = fi (εi ) = (2π)−1/2 σ −1 exp = (2π)−n/2 σ −n exp
i=1 i=1
2σ 2 2σ 2
In terms of the observed responses Y1 , . . . , Yn , we have:
Yi ∼ N ID(β0 + β1 Xi1 + · · · + βp Xip , σ 2 ) ⇒ Cov(Yi , Yj ) = 0 ∀i 6= j

From this, the joint probability density function for Y1 , . . . , Yn is:
n
Y Yn
−1/2 −1 −(yi − (β0 + β1 Xi1 + · · · + βp Xip ))2
f (y1 , y2 , . . . , yn ) = fi (yi ) = (2π) σ exp =
i=1 i=1
2σ 2
Pn
−n/2 −n − i=1 (yi − (β0 + β1 Xi1 + · · · + βp Xip ))2
= (2π) σ exp
2σ 2
The least squares estimates are Best Linear Unbiased Estimates (B.L.U.E.). Under normality
assumption, maximum likelihood estimates are Minimum Variance Unbiased Estimates (M.V.U.E.).
In either event, the estimate of β is:
For the LSD concentration/math score example, we have the following model for Y = X β + ε:
     
78.93 1 1.17 ε1
 58.20   1 2.97   ε2 
     
 67.47   1 3.26   ε3 
    β0  
 37.47  =  1 4.69   
    β1 +  ε4 
 45.65   1 5.83   ε5 
     
 32.92   1 6.00   ε6 
29.97 1 6.41 ε7
Note that β and ε are unobservable and must be estimated.
3.2 Normal Equations and Least Squares Estimates

Consider the matrices X0 X and X0 Y:
    Pn Pn 
1 1 ··· 1 1 X11 ··· X1p Pn n Pi=1
n
Xi1 ··· Pn i=1 Xip
 X11 X21 ··· Xn1  1 X21 ··· X2p   Xi1 2
i=1 Xi1 ··· i=1 Xi1 Xip

    i=1 
X0 X =  .. .. .. ..  .. .. .. .. = .. .. .. .. 
 .  .   . 
. . . . . . Pn . Pn . Pn . 2
X1p X2p ··· Xnp 1 Xn1 ··· Xnp i=1 Xip i=1 Xip Xi1 ··· i=1 Xip
For least squares estimation, we minimize Q( β), the error sum of squares with respect to β:
Q( β) = (Y − X β)0 (Y − X β) = Y 0 Y − Y 0 X β − β 0 X0 Y + β 0 X0 X β = Y 0 Y−2Y 0 X β + β 0 X0 X β
By taking the derivative of Q with respect to β, and setting this to 0, we get:

dQ( β)
= 0 − 2X0 Y + 2X0 X β̂ = 0 ⇒ X0 X β̂ = X0 Y
dβ
This leads to the normal equations and the least squares estimates (when the X matrix is of
full column rank).
X0 X β̂ = X0 Y ⇒ β̂ = (X0 X)−1 X0 Y
For the LSD concentration/math score example, we have the following normal equations and
least squares estimates:
" #" # " #
0 0 7 30.33 β̂0 350.61
X X β̂ = XY ⇒ =
30.33 153.8905 β̂1 1316.6558
" #
1 153.8905 −30.33
(X0 X)−1 =
7(153.8905) − (30.33)2 −30.33 7
" #" # " # " #
1 153.8905 −30.33 350.61 1 14021.3778 89.129
β̂ = = =
7(157.3246 −30.33 7 1316.6558 7(157.3246 −1417.4607 −9.0095
3.3 Fitted and Predicted Vectors

The vector of fitted (or predicted) values Ŷ is obtained as follows:
    
Ŷ1   1 X11 · · · X1p β̂0
  β̂0 + β̂1 X11 + · · · β̂p X1p   
 Ŷ2     1 X21 · · · X2p  β̂1 
Ŷi = 
 ..  =  β̂0 + β̂1 X21 + · · · β̂p X2p  = 
    .. .. .. .. 
 .. 

 .  ..  . . . .  . 
.β̂0 + β̂1 Xn1 + · · · β̂p Xnp
Ŷn 1 Xn1 · · · Xnp β̂p
= X β̂ = X(X0 X)−1 X0 Y = PY
Here, P is the projection of hat matrix, and is of dimension n × n. The hat matrix is symmetric
and idempotent:
P = X(X0 X)−1 X0 = (X(X0 X)−1 X0 )0 = P0 ⇒ Symmetric
PP = X(X0 X)−1 X0 X(X0 X)−1 X0 = X(X0 X)−1 X0 = P ⇒ Idempotent

For the LSD concentration/math score example, we have the following hat matrix (generated
in a computer matrix language):
 
1 1.17
 1 2.97 
 
 1 3.26 
  1 153.8905 −30.33 350.61 1 1 1 1 1 1 1
P=
 1 4.69 
 7(157.3246 =
  −30.33 7 1316.6558 1.17 2.97 3.26 4.69 5.83 6.00 6.41
 1 5.83 
 1 6.00 
1 6.41
 
0.58796 0.33465 0.29384 0.09260 −0.06783 −0.09176 −0.14946
 0.33465 0.22550 0.20791 0.12120 0.05207 0.04176 0.01690 
 
 0.29384 0.20791 0.19407 0.12581 0.07139 0.06327 0.04370 
 
=
 0.09260 0.12120 0.12581 0.14853 0.16665 0.16935 0.17586 
 −0.06783 0.05207 0.07139 0.16665 0.24259 0.25391 0.28122 
 
 −0.09176 0.04176 0.06327 0.16935 0.25391 0.26652 0.29694 
−0.14946 0.01690 0.04370 0.17586 0.28122 0.29694 0.33483
The vector of residuals, e is the vector generated by elementwise subtraction between the data
vector Y and the fitted vector Ŷ . It can be written as follows:
e = Y − Ŷ = Y − PY = (I − P)Y
Also, note:
Ŷ + e = PY + (I − P)Y = (P + I − P)Y = Y
For the LSD concentration/math score example, we have the following fitted and residual vec-
tors:
   
89.12 − 9.01(1.17) = 78.58 78.93 − 78.58 = 0.35
 89.12 − 9.01(2.97) = 62.36   58.20 − 62.36 = −4.16 
   
 89.12 − 9.01(3.26) = 59.75   67.47 − 59.75 = 7.72 
   
Ŷ = X β̂ = 
 89.12 − 9.01(4.69) = 46.86 
 e = Y − Ŷ = 
 37.47 − 46.86 = −9.39 

 89.12 − 9.01(5.83) = 36.59   45.65 − 36.59 = 9.06 
   
 89.12 − 9.01(6.00) = 35.06   32.92 − 35.06 = −2.14 
89.12 − 9.01(6.41) = 31.37 29.97 − 31.37 = −1.40
3.4 Properties of Linear Functions of Random Vectors

Note that β̂, Ŷ , and e are all linear functions of the data vector Y, and can be written as AY:
• β̂ = (X0 X)−1 X0 Y ⇒ A = (X0 X)−1 X0

• Ŷ = PY ⇒ A = P = X(X0 X)−1 X0
• e = (I − P)Y ⇒ A = I − P = I − X(X0 X)−1 X0
Consider a general vector Z that is of dimension 3 × 1. This can be easily expanded to n × 1,

but all useful results can be seen in the simpler case.
 
z1
 
Z =  z2 
z3
The expectation vector is the vector made up of the elementwise expected values of the
elements of the random vector.
   
E(z1 ) µ1
   
E[Z] =  E(z2 )  =  µ2  = µz
E(z3 ) µ3
Note that the matrix (Z − µz )(Z − µz )0 is 3 × 3, and can be written as:

 
(z1 − µ1 )2 (z1 − µ1 )(z2 − µ2 ) (z1 − µ1 )(z3 − µ3 )
 
(z
 2 − µ 2 )(z 1 − µ 1 ) (z2 − µ2 )2 (z2 − µ2 )(z3 − µ3 ) 
(z3 − µ3 )(z1 − µ1 ) (z3 − µ3 )(z2 − µ2 ) (z3 − µ3 )2
The variance-covariance matrix is the 3 × 3 matrix made up of variances (on the main
diagonal) and the covariances (off diagonal) of the elements of Z.
 
V ar(z1 ) Cov(z1 , z2 ) Cov(z1 , z3 )
 
Var[Z] =  Cov(z2 , z1 ) V ar(z2 ) Cov(z2 , z3 )  = Vz =
Cov(z3 , z1 ) Cov(z3 , z2 ) V ar(z3 )
 
E[(z1 − µ1 )2 ] E[(z1 − µ1 )(z2 − µ2 )] E[(z1 − µ1 )(z3 − µ3 )]
 
 E[(z2 − µ2 )(z1 − µ1 )] E[(z2 − µ2 )2 ] E[(z2 − µ2 )(z3 − µ3 )]  =
E[(z3 − µ3 )(z1 − µ1 )] E[(z3 − µ3 )(z2 − µ2 )] 2
E[(z3 − µ3 ) ]
= E[(Z − µz )(Z − µz )0 ] = Vz
Now let A be a k × n matrix of constants and z be a n × 1 random vector with mean vector µz ,
and variance-covariance matrix Vz . Suppose further that we can write A and U = Az as follow:
 
a01
 
 a02 
A =  .. 
 
 . 
a0k
   
a01 z u1
   
 a02 z   u2 
U = Az = 
 .. =
  .. 

 .   . 
a0k z uk
where a0i is a 1 × n row vector of constants.
To obtain E[U] = µu , consider each element of U, namely ui and it’s expectation E(ui ).
E[ui ] = E[a0i z] = E[ai1 z1 +ai2 z2 +· · ·+ain zn ] = ai1 E[z1 ]+ai2 E[z2 ]+· · ·+ain E[zn ] = a0i µz i = 1, . . . , k
Piecing these together, we get:

     
E[u1 ] a01 µz a01
     
 E[u2 ]   a02 µz   a02 
E[U] =  .. = .. = ..  µz = Aµz = µu
     
 .   .   . 
E[uk ] a0k µz a0k
To obtain the variance covariance matrix of U = Az, first consider the definition of V[U], then
write in terms of U = Az:
Var[U] = Vu = E[(U − µu )(U − µu )0 ] =
= E[(Az − Aµz )(Az − Aµz )0 ] = E{[A(z − µz )][A(z − µz )]0 } =

= E[A(z − µz )(z − µz )0 A0 ] = AE[(z − µz )(z − µz )]A0 = AVz A0
Note that if Vz = σ 2 I, then Vu = σ 2 AA0 .
3.5 Applications of Linear Functions of Random Variables

In this section, we consider two applications, each assuming independent observations (Cov[Yi , Yj ] =
0 i 6= j).
Case 1 – Sampling from a single population
E[Yi ] = µ i = 1, . . . , n V ar[Yi ] = σ 2 i = 1, . . . , n
Let Y be the n × 1 vector made up of elements Y1 , . . . , Yn . Then:

   
E[Y1 ] µ
   
 E[Y2 ]   µ 
E[Y] = 
 .. =
  ..  = µ1

 .   . 
E[Yn ] µ
where 1 is a n × 1 column vector of 10 s.
   
V ar[Y1 ] Cov[Y1 , Y2 ] · · · Cov[Y1 , Yn ] σ2 0 · · · 0
   
 Cov[Y2 , Y1 ] V ar[Y2 ] · · · Cov[Y2 , Yn ]   0 σ2 · · · 0 
Var(Y) = 
 .. .. .. .. =
  .. .. . . ..  = σ2 I

 . . . .   . . . . 
Cov[Yn , Y1 ] Cov[Yn , Y2 ] · · · V ar[Yn ] 0 0 ··· 0
Now consider the estimator Y :

Pn h i
i=1 Yi
Y = = 1
n
1
n ··· 1
n Y = a0 Y
n
h i
1 1 1
Now we can obtain the mean and variance of Y from these rules, with a0 = n n ··· n :
 
µ
n

 µ 
 X 1 1
E[Y ] = a0 E[Y] = 1
n
1
n ··· 1
n  ..  = a‘µ1 = µ=n µ=µ
 .  i=1
n n
µ
 1

n n
X 2 2
 
1 1 1 σ2
V ar[Y ] = a0 Var[Y]a = 1 1
··· 1
σ2 I 

 = σ 2 a0 a = σ 2
n = σ2 n =
n n n ···  n n n
1 i=1
n
Case 2 – Multiple Linear Regression Model
E[Yi ] = β0 + β1 Xi1 + · · · + βp Xip i = 1, . . . , n V ar[Yi ] = σ 2 i = 1, . . . , n
Let Y be the n × 1 vector made up of elements Y1 , . . . , Yn . Then:

 
E[Y1 ]
 
 E[Y2 ] 
E[Y] = 
 .. =Xβ

 . 
E[Yn ]
Var[Y] = σ2 I
Now consider the least squares estimator β̂ of the regression coefficient parameter vector β.
β̂ = (X0 X)−1 X0 Y = A0 Y ⇒ A0 = (X0 X)−1 X0
The mean and variance of β̂ are:
E[ β̂] = E[A0 Y] = A0 E[Y] = (X0 X)−1 X0 X β = β
Var[ β̂] = Var[A0 Y] = A0 Var[Y]A = (X0 X)−1 X0 σ 2 IX(X0 X)−1 = σ 2 (X0 X)−1 (X0 X)(X0 X)−1 = σ 2 (X0 X)−1
3.6 Multivariate Normal Distribution
Suppose a n × 1 random vector Z has a multivariate normal distribution with mean vector 0 and
variance-covariance matrix σ 2 I. This would occur if we generated n independent standard normal
random variables and put them together in vector form. The density function for Z, evaluated any
fixed point z is:
1
Z ∼ N ID(0, σ 2 I) ⇒ fZ (z) = (2π)−n/2 |σ 2 I|−1/2 exp{− z0 (σ 2 I)−1 z}
2
More generally, let U = AZ + b with A being a k × n matrix of constants, and b being a n × 1

vector of constants. Then:
E[U] = AE[Z] + b = b = µU V[U] = AV[Z]A0 = σ 2 AA0 = VU
The density function for U, evaluated any fixed point u is:

1
U ∼ N ID(µU , σ 2 VU ) ⇒ fU (u) = (2π)−k/2 |VU |−1/2 exp{− (u − µU )0 (VU )−1 (u − µU )}
2
That is, any linear function of a normal random vector is normal.
3.6.1 Properties of Regression Estimates

Under the traditional normal theory linear regression model, we have:
Y=X β+ ε ε ∼ N(0, σ 2 I) ⇒ Y ∼ N(X β, σ 2 I)
Then the density function of Y evaluated is:

1
fY (y) = (2π)−n/2 σ −n exp{− (y − X β)0 (y − X β)}
2σ 2
Assuming the model is correct, we’ve already obtained the mean and variance of β̂. We further
know that its distibution is multivariate normal: β̂ ∼ N (β, σ 2 (X0 X)−1 ).
The vector of fitted values Ŷ = X β̂ = X(X0 X)−1 X0 Y = PY is also a linear function of Y and
thus also normally distributed with the following mean vector and variance-covariance matrix:
E[Ŷ] = X(X0 X)−1 X0 E[Y] = X(X0 X)−1 X0 X β = X β
Var[Ŷ] = X(X0 X)−1 X0 σ 2 IX(X0 X)−1 X0 = σ 2 X(X0 X)−1 X0 X(X0 X)−1 X0 = σ 2 X(X0 X)−1 X0 = σ 2 P
That is Ŷ ∼ N (X β, σ 2 P).
The vector of residuals e = Y − Ŷ = (I − P)Y is also normal with mean vector and variance-
covariance matrix:
E[e] = (I − P)X β = (X − PX) β = (X − X(X0 X)−1 X0 X) β = (X − X) β = 0
Var[e] = (I − P)σ 2 I(I − P)0 = σ 2 (I − P)(I − P)0 = σ 2 (I − P)(I − P) = σ 2 (I − P)
Note the differences between the distributions of ε and e:
ε ∼ N ( 0, σ 2 I) ⇒ e ∼ N ( 0, σ 2 (I − P))
Often the hgoal is to predict a ifuture outcome when the set of independent levels are at a given
setting, x00 = 1 x01 · · · x0p . The future observation Y0 and its predicted value based on the
estimated regression equation are:
Y0 = x00 β + ε0 ε0 ∼ N ID(0, σ 2 )
Ŷ0 = x00 β̂ Ŷ0 ∼ N (x00 β, σ 2 x00 (X0 X)−1 x0 )
It is assumed ε0 ∼ N (0, σ 2 ) and is independent from the errors in the observations used to fit
the regression model (ε1 , . . . , εn ).
The prediction error is:
Y0 − Ŷ0 = x00 β + ε0 − x00 β̂ = x00 ( β − β̂) + ε0
which is normal with mean and variance:
E[Y0 − Ŷ0 ] = x00 ( β − E[ β̂]) + E[ε0 ] = x00 ( β − β) + 0 = 0
V [Y0 − Ŷ0 ] = V [Y0 ] + V [Ŷ0 ] = σ 2 + σ 2 x00 (X0 X)−1 x0 = σ 2 [1 + x00 (X0 X)−1 x0 ]
and we have that:
Y0 − Ŷ0 ∼ N (0, σ 2 [1 + x00 (X0 X)−1 x0 ])

4 Analysis of Variance and Quadratic Forms
The sums of square in the Analysis of Variance can be written as quadratic forms in Y. The
form we use is Y 0 AY where A is a matrix of coefficients, referred to as the defining matrix.
The following facts are important and particularly useful in regression models (for a very detailed
discussion, see Linear Models (1971), by S.R. Searle).
1. Any sum of squares can be written as Y 0 AY where A is a square, symmetric nonnegative

definite matrix
2. The degrees of freedom associated with any quadratic form is equal to the rank of the defining
matrix, which is equal to its trace when the defining matrix is idempotent.
3. Two quadratic forms are orthogonal if the product of their defining matices is 0
4.1 The Analysis of Variance

Now consider the Analysis of Variance.
n
X
0
Y = Ŷ + e YY= Yi2 = SS(TOTAL UNCORRECTED)
i=1
Note that Y 0 Y = Y 0 IY, so that I is the defining matrix, which is symmetric and idempotent. The
degrees of freedom for SS(TOTAL UNCORRECTED) is then the rank of I, which is its trace, or
n.
Now, we decompose the Total uncorrected sum of squares into it’s model and error components.
Y0 Y = (Ŷ + e)0 (Ŷ + e) = Ŷ0 Ŷ + Ŷ0 e + e0 Ŷ + e0 e =
= (PY)0 (PY) + (PY)0 (I − P)Y + [(I − P)Y]0 PY + [(I − P)Y]0 [(I − P)Y] =
= Y0 P0 PY + Y0 P0 (I − P)Y + Y0 (I − P)0 PY + Y0 (I − P)0 (I − P)Y =
= Y0 PPY + (Y0 PY − Y0 PPY) + (Y0 PY − Y0 PPY) + (Y0 IY − Y0 IPY − Y0 PIY + Y0 PPY) =
= Y0 PY + (Y0 PY − Y0 PY) + (Y0 PY − Y0 PY) + (Y0 Y − Y0 PY − Y0 PY + Y0 PY) =
= Y0 PY + 0 + 0 + (Y0 Y − Y0 PY) = Y0 PY + Y0 (I − P)Y =
= Y0 P0 PY + Y0 (I − P)0 (I − P)Y = Ŷ0 Ŷ + e0 e
We obtain the degrees of freedom as follow, making use of the following identities regarding the
trace of matrices:
tr(AB) = tr(BA) tr(A + B) = tr(A) + tr(B)
SS(MODEL) = Ŷ 0 Ŷ = Y 0 PY
df (MODEL) = tr(P) = tr(X(X0 X)−1 X0 ) = tr((X0 X)−1 X0 X) = tr(Ip0 ) = p0 = p + 1
SS(RESIDUAL) = e0 e = Y 0 (I − P)Y
df (RESIDUAL) = tr(I − P) = tr(In ) − tr(P) = n − p0 = n − p − 1
Table 8 gives the Analysis of Variance including degrees of freedom, and sums of squares (both
definitional and computational forms).
Source of Degrees of Sum of Squares
Variation Freedom Definitional Computational
TOTAL(UNCORRECTED) n YY0 Y0 Y
0 0
MODEL p =p+1 Ŷ 0 Ŷ = Y 0 PY β̂ X0 Y
0
ERROR n − p0 e0 e = Y 0 (I − P)Y Y 0 Y− β̂ X0 Y
Table 8: Analysis of Variance in Matrix form
We obtain the Analysis of Variance in matrix form:

 
78.93
 58.20 
 
 67.47 
 
SS(TOTAL UNCORRECTED) = Y Y = 78.93 58.20 67.47 37.47 45.65 32.92 29.97 
0
 37.47 =

 45.65 
 
 32.92 
29.97
n
X
= Yi2 = 19639.2365 df (TOTAL UNCORRECTED) = n = 7
i=1

0 0
350.61
SS(MODEL) = β̂ X Y = 89.1239 −9.0095 =
1316.6558
89.1239(350.61) + (−9.0095)(1316.6558) = 19385.3201 df (MODEL) = p0 = 2
0
SS(RESIDUAL) = Y0 Y− β̂ X0 Y = 19639.2365−19385.3201 = 253.9164 df (RESIDUAL) = n−p0 = 7−2 = 5
The total uncorrected sum of squares represents variation (in Y ) around 0. We are usually
interested in variation around the sample mean Y . We will partition the model sum of squares into
two components: SS(REGRESSION) and SS(µ). The first sum of squares is associated with β1
and the second one is associated with β0 .
Model with only the mean µ = β0 (β1 = 0)
Consider the following model, we obtain the least squares estimates and model sum of squares.
Yi = β 0 + εi = µ + εi
 
1
 1 
 
Y=X β+ ε X= .. =1 β = [µ] = [β0 ]
 . 
1
X
n
β̂ = (X0 X)−1 X0 Y = (10 1)−1 10 Y 10 1 = n 1 0 Y = Yi
i=1
Pn
Yi
⇒ β̂ = (10 1)−1 10 Y = i=1
=Y
n
n
X
0 2
SS(µ) = β̂ X0 Y = Y ( Yi ) = nY
i=1
1
= Y0 1(10 1)−1 10 Y = Y0 ( 110 )Y =
n
 1 1 
n n · · · n1
 1 1 ··· 1 
 n n n  0 1
Y0  . . . ..  Y = Y ( J)Y
.
 . . . . . .  n
1 1 1
n n ··· n
where J is a n × n matrix of 1s . Note that (1/n)J is an idempotent matrix:

 
n n ··· n
2 2  n n ··· n 
1 1 1 1   1
( J)( J) = JJ =  .. .. . . .. = J
n n n n  . . . .  n
n n ··· n
By subtraction, we get SS(REGRESSION) = SS(MODEL) − SS(µ):

1 1
SS(REGRESSION) = SS(MODEL) − SS(µ) = Y0 PY − Y0 ( J)Y = Y0 (P− J)Y
n n
To demonstrate that the defining matrix for SS(REGRESSION) is idempotent and that the
three sum of squares are orthogonal, consider the following algebra where X∗ is the matrix made
up of the columns of X associated with the p independent variables and not the column for the
intercept.
X = [1|X∗ ] PX = P[1|X∗ ] = X = [1|X∗ ]
⇒ P1 = 1 ⇒ PJ = J
0 ∗ 0
X = [1|X ] X P = [1|X ] P = X0 = [1|X∗ ]0
0 ∗ 0
10 P = 10 ⇒ JP = J
⇒
1 1 1 1 1 1 1 1 1 1
(P− J)(P− J) = PP − P( J)− JP + ( J)( J) = P − ( J) − ( J) + ( J) = P− J
n n n n n n n n n n
Summarizing what we have obtained so far (where all defining matrices are idempotent):
SS(TOTAL UNCORRECTED) = Y0 IY = Y0 Y df (TOTAL UNCORRECTED) = tr(In ) = n
1 1 1
SS(µ) = Y0 ( J)Y df (µ) = tr( J) = (n) = 1
n n n
1 1 1
SS(REGRESSION) = Y0 (P− J)Y df (REGRESSION) = tr((P− J) = tr(P)−tr( J) = p0 −1 = p+1−1 = p
n n n
SS(RESIDUAL) = Y0 (I − P)Y df (RESIDUAL) = tr(I − P) = tr(I) − tr(P) = n − p0
To show that the sums of squares for the mean, regression, and residual are pairwise orthogonal,
consider the products of their defining matrices: First for SS(µ) and SS(REGRESSION):
1 1 1 1 1 1 1
( J)(P− J) = JP − ( J)( J) = J − J = 0
n n n n n n n
Next for SS(µ) and SS(RESIDUAL):

1 1 1 1 1
( J)(I − P) = JI − JP = J − J = 0
n n n n n
Finally for SS(REGRESSION) and SS(RESIDUAL):

1 1 1 1 1
(P− J)(I − P) = PI − PP− JI+ JP = P − P − J+ J = 0
n n n n n
For the LSD concentration/math score example, we have the ANOVA in Table 9.
Source of Degrees of Sum of Mean

Variation Freedom Squares Square
TOTAL(UNCORRECTED) 7 19639.24 —
MEAN 1 17561.02 —
TOTAL (CORRECTED) 6 2078.22 —
REGRESSION 1 1824.30 1824.30
RESIDUAL 5 253.92 50.78
Table 9: Analysis of Variance for LSD data
A summary of key points regarding quadratic forms:
• The rank, r(X) is the number of linearly independent columns in X
• The model is full rank if r(X) = p0 assuming n > p0
• A unique least squares solution exists iff the model is full rank.
• All defining matrices in the Analysis of Variance are idempotent.
• The defining matrices for the mean, regression, and residual are pairwise orthogonal and
sum to I. Thus they partition the total uncorrected sum of squares into orthogonal sums of
squares.
• Degrees of freedom for quadratic forms are the ranks of their defining matrices; when idem-
potent, the trace of a matrix is its rank.
4.2 Expectations of Quadratic Forms

In this section we obtain the expectations of the sums of squares in the Analysis of Variance, making
use of general results in quadratic forms. The proofs are given in Searle (1971). Suppose we have
a random vector Y with the following mean vector and variance-covariance matrix:
E[Y] = µ Var[Y] = VY = Vσ 2
Then, the expectation of a quadratic form Y 0 AY is:
E[Y 0 AY] =tr(AVY ) + µ0 A µ = σ 2 tr(AV) + µ0 A µ
Under the ordinary least squares assumptions, we have:
E[Y] = X β Var[Y] =σ 2 In
Source of Variation “A” Matrix

TOTAL UNCORRECTED I
MODEL P = X(X0 X)−1 X0
REGRESSION P− n1 J
RESIDUAL I−P
Now applying the rules on expectations of quadratic forms:
E[SS(MODEL)] = E[Y0 PY] = σ 2 tr(PI) + β 0 X0 PX β =
= σ 2 tr(P) + β0 X0 X(X0 X)−1 X0 X β = σ 2 p0 + βX0 X β
1 1 1
E[SS(REGRESSION)] = E[Y0 (P− J)Y] = σ 2 tr(P− J) + β 0 X0 (P− J)X β =
n n n
1 1
σ 2 (p0 − 1) + β 0 X0 X β − β 0 X0 JX β = pσ 2 + β 0 X0 (I− J)X β
n n
This last matrix can be seen to involve the regression coefficients: β1 , . . . , βp , and not β0 as
follows:
1 1
X0 (I− J)X = X0 X − X0 JX
n n
 Pn Pn 
Pn n Pi=1
n
Xi1 · · · P i=1 Xip
n
 Xi1 2
i=1 Xi1 ··· i=1 Xi1 Xip

 i=1 
X0 X =  .. .. .. .. 
 . 
Pn . Pn . Pn . 2
i=1 Xip i=1 Xip Xi1 ··· i=1 Xip
 
Pn n Pn n ··· Pn n
1 0 1
 i=1 Xi1 i=1 Xi1 ··· i=1 Xi1 

X JX =  .. .. .. .. X =
n n . 
Pn . Pn . Pn .
i=1 Xip i=1 Xip ··· i=1 Xip
 Pn Pn 
n2 n i=1 Xi1 ··· n i=1 Xip
Pn Pn 2 Pn Pn
1
 n i=1 Xi1 ( i=1 Xi1 ) ··· ( i=1 Xi1 ) ( i=1 Xip )


=  .. .. .. .. =
n . . . . 
Pn Pn Pn Pn 2
n i=1 Xip ( i=1 Xip ) ( i=1 Xi1 ) ··· ( i=1 Xip )
 Pn Pn 
n Pi=1 Xi1 ··· Pn Xip
i=1 P
 Pn n
( i=1 Xi1 )
2
( Xi1 )(
n
Xip ) 
 
i=1 Xi1 ··· i=1 i=1
 n n 
= .. .. .. .. 
 . .P . 
 Pn Pn . 
Pn ( Xip )(
n
Xi1 ) ( Xip )
2
i=1 Xip ···

i=1 i=1 i=1
n n
1 1
⇒ X0 (I− J)X = X0 X − X0 JX =
n n
 Pn Pn 
 Pn Pn  n Xi1 ··· Xip
Pi=1 Pn i=1 P
Pn n Pi=1
n
Xi1 ··· Pn i=1 Xip  Pn n
( i=1 Xi1 )
2
( Xi1 )(
n
Xip ) 
 Xi1 2
i=1 Xi1 ··· i=1 Xi1 Xip
  
i=1 Xi1 ··· i=1 i=1
 i=1   n n 
= .. .. .. .. − .. .. .. .. 
 .   . 
Pn . Pn . Pn . 2  . Pn .P Pn . 
i=1 Xip i=1 Xip Xi1 ··· i=1 Xip
Pn ( Xip )(
n
Xi1 ) ( Xip )
2
i=1 Xip ···

i=1 i=1 i=1
n n
 
0 0P ··· 0
P Pn
 Pn n
( i=1 Xi1 )
2
Pn (
n
Xi1 )( Xip ) 
 0 X 2
− ··· Xi1 Xip − i=1 i=1 
 i=1 i1 n i=1 n 
= . .. .. .. =
 .. . . . P 
 P Pn 
Pn n
( i=1 Xip )( i=1 Xi1 ) Pn 2
n
( i=1 Xip )
2
0 i=1 Xip Xi1 − n ··· i=1 Xip − n
 
0 Pn 0 ··· Pn 0
 0 − X 1 )2 
 i=1 (Xi1 ··· i=1 (Xi1 − X 1 )(Xip − X p ) 
= .. .. .. .. =
 . . . . 
Pn Pn 2
0 i=1 (X ip − X p )(Xi1 − X 1 ) ··· i=1 (X ip − X p )
Thus, E[SS(REGRESSION)] involves a quadratic form in β1 , . . . , βp , and not in β0 since the first
row and column of the previous matrix is made up of 0s . Now, we return to E[SS(RESIDUAL)]:
E[SS(RESIDUAL)] = E[Y0 (I − P)Y] = σ 2 tr(I − P) + β0 X0 (I − P)X β =
σ 2 (n − p0 ) + βX0 X β − βX0 PX β = σ 2 (n − p0 ) + βX0 X β − βX0 X β = σ 2 (n − p0 )
Now we can obtain the expected values of the mean squares from the Analysis of Variance:
SS(REGRESSION) 1 0 0 1
M S(REGRESSION) = ⇒ E[M S(REGRESSION)] = σ 2 + β X (I− J)X β
p p n
SS(RESIDUAL)
M S(RESIDUAL = ⇒ E[M S(RESIDUAL)] = σ 2
n − p0
Note that the second term in E[M S(REGRESSION)] is a quadratic form in β, if any βi 6= 0
(i = 1, . . . , p), then E[M S(REGRESSION)] > E[M S(RESIDUAL)], otherwise they are equal.
4.2.1 The Case of Misspecified Models

The above statements presume that the model is correctly specified. Suppose:
Y =X β+Z γ+ ε ε ∼ N (0, σ 2 I)
but we fit the model:

Y=X β+ ε
Then E[M S(RESIDUAL)] can be written as:
E[SS(RESIDUAL)] = σ 2 tr(I − P) + (X β + Z γ)0 (I − P)(X β + Z γ) =
= σ 2 (n − p0 ) + β0 X0 (I − P)X β + β 0 X0 (I − P)Z γ + γ 0 Z0 (I − P)X β + γ 0 Z0 (I − P)Z γ =
= σ 2 (n − p0 ) + 0 + 0 + 0 + γ 0 Z0 (I − P)Z γ since X0 (I − P) = (I − P)X = 0
E[SS(RESIDUAL)] 1
E[M S(RESIDUAL)] = 0
= σ2 + γ 0 Z0 (I − P)Z γ
n−p n − p0
which is larger than σ 2 if the elements of γ are not all equal to 0 (which would make our fitted
model correct).
Theoretical Estimated
Estimator Variance Variance
β̂ σ 2 (X0 X)−1 s (X0 X)−1
2
Ŷ σ 2 X(X0 X)−1 X0 = σ 2 P s2 X(X0 X)−1 X0 = s2 P

e σ 2 (I − P) s2 (I − P)
Table 10: Theoretical and estimated variances of regression estimators in Matrix form
4.2.2 Estimated Variances

Recall the variance-covariance matrices of β̂, Ŷ, and e. Each of these depended on σ 2 , which is
in practice unknown. Unbiased estimators of these variances can be obtained by replacing σ 2 with
an unbiased estimate:
1
s2 = M S(RESIDUAL) = Y 0 (I − P)Y
n − p0
assuming the model is correct. Table 10 gives the true and estimated variances for these estimators.
4.3 Distribution of Quadratic Forms

We have obtained means of quadratic forms, but need their distributions to make statistical infer-
ences. The assumptions for the traditional inferences to be made is that ε and Y are normally
distributed, otherwise tests are approximate.
The following results are referred to as Cochran’s Theorem, see Searle (1971) for proofs. Suppose
Y is distributed as follows with nonsingular matrix V:
Y ∼ N (µ, Vσ 2 ) r(V) = n
then:

1
1. Y 0 σ2 A Y is distributed noncentral χ2 with:
(a) Degrees of freedom = r(A)

1
(b) Noncentrality parameter = Ω = σ2
µ0 A µ if AV is idempotent
2. Y 0 AY and Y 0 BY are independent if AVB = 0
3. Y 0 AY and linear function BY are independent if BVA = 0
4.3.1 Applications to Normal Multiple Regression Model

The sums of squares for the Analysis of Variance are all based on idempotent defining matrices:
For the Model sum of squares:

SS(MODEL) 1
= Y0 P Y AV = PI = P AVAV = PP = P
σ2 σ2
df (MODEL) = r(A) = r(P) = p0

1 1
Ω(MODEL) = 2
β 0 X0 PX β = β 0 X0 X β
2σ 2σ 2
For the Mean sum of squares:

SS(µ) 1 1 1 1 1 1 1
= Y0 J Y AV = JI = J AVAV = J J = J
σ2 σ2 n n n n n n
n
1 1X 1
df (MEAN) = r( J) = 1= n=1
n n i=1 n
1 1 1 (10 X β)2
Ω(MEAN) = 2
β 0 X0 JX β =
2σ n 2σ 2 n
The last equality is obtained by recalling that J = 110 , and:
β 0 X0 JX β = β 0 X0 110 X β = ( β 0 X0 1)(10 X β) = (10 X β)2
For the Regression sum of squares:

SS(REGRESSION) 1 1 1
= Y0 (P− J) Y AV = P− J)I
σ2 σ2 n n
1 1 1 1 1
AVAV = PP − P J − JP + J J = P− J
n n n n n
1 1
df (REGRESSION) = r(P− J)) = r(P) − r( J) = p0 − 1
n n
1 1 1 1 1 1
Ω(REGRESSION) = 2 β 0 X0 (P− J)X β = 2 β 0 X0 (P− J)X β = 2 β 0 X0 (I− J)X β
2σ n 2σ n 2σ n
For the Residual sum of squares:

SS(RESIDUAL) 1
= Y0 (I − P) Y AV = (I − P)I = (I − P) AVAV = (I − P)(I − P) = (I − P)
σ2 σ2
df (RESIDUAL) = r(A) = r((I − P)) = r(I) − r(P) = n − p0

1 1
Ω(RESIDUAL) = β0 X0 (I − P)X β = ( β0 X0 X β − β0 X0 X β) = 0
2σ 2 2σ 2
Since we have already shown that the quadratic forms for SS(µ), SS(REGRESSION), and
SS(RESIDUAL) are all pairwise orthogonal, and in our current model V = I, then these sums of
squares are all independent due to the second part of Cochran’s Theorem.
Consider any linear function K0 β̂ = K0 (X0 X)−1 X0 Y = BY. Then, by the last part of Cochran’s
Theorem, K0 β̂ is independent of SS(RESIDUAL):
B = K0 (X0 X)−1 X0 V=I A = (I − P)
⇒ BVA = K0 (X0 X)−1 X0 − K0 (X0 X)−1 X0 P =
K0 (X0 X)−1 X0 − K0 (X0 X)−1 X0 X(X0 X)−1 X0 = K0 (X0 X)−1 X0 − K0 (X0 X)−1 X0 = 0
Consider the following random variable F :
X12 /ν1
F =
X22 /ν2
where X12 is distibuted noncentral χ2 with ν1 degrees of freedom and noncentrality parameter Ω1 ,
and X22 is distibuted central χ2 with ν2 degrees of freedom. Further, assume that X12 and X22 are
independent. Then, F is distrbuted noncental F with ν1 numerator, ν2 denominator degrees of
freedom, and noncentrality parameter Ω1 .
This applies as follows for the F -test in the Analysis of Variance.

SS(REGRESSION) β 0 X0 (P− n1 J)X β
• σ2 ∼ noncentral-χ2 with df = p and Ω = 2σ 2
SS(RESIDUAL)
• σ2 ∼ central-χ2 with df = n − p0
• SS(REGRESSION) and SS(RESIDUAL) are independent

SS(REGRESSION)
/p
σ 2
S(REGRESSION)
• = MM S(RESIDUAL)
∼ noncentral-F with p numerator and n − p0 denom-
SS(RESIDUAL) 0
/(n−p )
σ2
β 0 X0 (P− n1 J)X β
inator degrees’ of freedom, and noncentrality parameter Ω = 2σ 2
• The noncentrality parameter for SS(REGRESSION) does not involve β0 , and for full rank X, Ω =
0 ⇐⇒ β1 = β2 = · · · = βp = 0, otherwise Ω > 0
This theory leads to the F -test to determine whether the set of p regression coefficients β1 , β2 , . . . , βp
are all equal to 0:
 
β1
 
 β2 
• H0 : β = 0 where β = 
∗

∗
.. 

 . 
βp
• HA : β ∗ 6= 0
SS(REGRESSION)/p M S(REGRESSION)
• T S : F0 = M S(RESIDUAL)/(n−p0 )
= M S(RESIDUAL)
• RR : F0 :≥ F(α,p,n−p0 )
• P -value: P r{F ≥ F0 } where F ∼ Fp,n−p0
• The power of the test under a specific alternative can be found by finding the area under
the relevant noncentral-F distribution to the right of the critical value defining the rejection
region.
Example 1 – LSD Pharmacodynamics
Suppose that the true parameter values are: β0 = 90, β1 = −10, and σ 2 = 50 (these are
consistent with the least squares estimates). Recall that the fact β0 6= 0 has no bearing on the
F -test, only that β1 6= 0. Then:
β 0 X0 (P− n1 J)X β
Ω= = 22.475
2σ 2
Figure 3 gives the central-F distribution (the distribution of the test statistic under the null hy-
pothesis) and the noncentral-F distribution (the distribution of the test statistic under this specific
alternative hypothesis). Further, the power of the test under these specific parameter levels is the
area under the noncentral-F distribution to the right of Fα,1,5 . Table 4.3.1 gives the power (the
probability we reject H0 under H0 and several sets of values in the alternative hypothesis) for three
levels of α, where F(.100,1,5) = 4.06, F(.050,1,5) = 6.61, and F(.010,1,5) = 16.26. The reason that
the column for the noncentrality parameter is 2Ω is that SAS’ function for returning
a tail area from a noncentral-F distribution is twice the noncentrality parameter we
use in this section’s notation.)
fho_act
0.020
0.018
0.016
0.014
0.012
0.010
0.008
0.006
0.004
0.002
0.000
0 10 20 30 40 50 60 70 80 90 100
f
Figure 3: Central and noncentral-F distributions for LSD example, β1 = −10, σ 2 = 50
Note that as the true slope parameter moves further away from 0 for a fixed α level, the power
increases. Also, as α (the size of the rejection region) decreases, so does the power of the test.
Under the null hypothesis (β1 = 0), the size of the rejection region is the power of the test (by
definition).
α
β1 2Ω 0.100 0.050 0.010
0 0 0.100 0.050 0.010
-2 1.80 0.317 0.195 0.053
-4 7.19 0.744 0.579 0.239
-6 16.18 0.962 0.890 0.557
-8 28.77 0.998 0.987 0.831
-10 44.95 1.000 0.999 0.959
Table 11: Power = Pr(Reject H0 ) under several configurations of Type I error rate (α) and slope
parameter (β1 ) for σ 2 = 50
4.4 General Tests of Hypotheses

Tests regarding linear functions of regression parameters are conducted as follow.
• Simple Hypothesis ⇒ One linear function
• Composite Hypothesis ⇒ Several linear functions
H 0 : K0 β = m HA : K0 β 6= m
where K0 is a k × p0 matrix of coefficients defining k linear functions of the βjs to be tested (k ≤ p0 ),

and m is a k × 1 column vector of constants (often, but not necessarily 0s ). The k linear functions
must be linearly independent, but need not be orthogonal. This insures that K0 will be full (row)
rank (that is r(K0 ) = k) and that H0 will be consistent ∀m.
Estimator and its Variance
Parameter – K0 β − m
Estimator – K0 β̂ − m E[K0 β̂ − m] = K0 β − m
Variance of Estimator – Var[K0 β̂ − m] = Var[K0 β̂] = K0 Var[ β̂]K = K0 σ 2 (X0 X)−1 K = Vσ 2
Sum of Squares for testing H0 : K 0 β = m

A quadratic form is created from the estimator K0 β̂ − m by using a defining matrix that is
the inverse of V. This is can be thought of as a matrix version of “squaring a t-statistic.
Q = (K0 β̂ − m)0 [K0 (X0 X)−1 K]−1 (K0 β̂ − m)
That is, Q is a quadratic form in K0 β̂ − m with A = [K0 (X0 X)−1 K]−1 = V−1 . Making use of the
earlier result, regarding expectations of quadratic forms, namely:
E[Y 0 AY] =tr(AVY ) + µ0 A µ = σ 2 tr(AV) + µ0 A µ

we get:
E[Q] = σ 2 tr[(K0 (X0 X)−1 K)−1 (K0 (X0 X)−1 K)] + (K0 β − m)0 [K0 (X0 X)−1 K]−1 (K0 β − m) =
= σ 2 tr[Ik ] + K0 β − m)0 [K0 (X0 X)−1 K]−1 (K0 β − m) = kσ 2 + (K0 β − m)0 [K0 (X0 X)−1 K]−1 (K0 β − m)
Now, AV = Ik is idempotent and r(A) = r(K) = k (with the restrictions on K0 stated above).
So as long as ε holds our usual assumptions (normality, constant variance, independent elements),
then Q/σ 2 is distributed noncentral-χ2 with k degrees of freedom and noncentrality parameter:
µ0 Aµ (K0 β − m)0 [K0 (X0 X)−1 K]−1 (K0 β − m)

ΩQ = =
2σ 2 2σ 2
where ΩQ = 0 ⇐⇒ K0 β = m
So, as before, for the test of β ∗ = 0, we have a sum of squares for a hypothesis that is
noncentral-χ2 , in this case having k degrees of freedom. Now, we show that Q is independent
of SS(RESIDUAL), for the case m = 0 (it holds regardless, but the math is messier otherwise).
Q = (K0 β̂)0 [K0 (X0 X)−1 K]−1 (K0 β̂) SS(RESIDUAL) = Y0 (I − P)Y
0
Q = β̂ K[K0 (X0 X)−1 K]−1 (K0 β̂) = Y0 X(X0 X)−1 K[K0 (X0 X)−1 K]−1 (K0 (X0 X)−1 X0 Y
Recall that Y 0 AY and Y 0 BY are independent if BVA = 0. Here, V = I.
BA = X(X0 X)−1 K[K0 (X0 X)−1 K]−1 (K0 (X0 X)−1 X0 (I − P) = 0
since X0 P = X0 X(X0 X)−1 X0 = X0 . Thus Q is independent of SS(RESIDUAL). This leads to the

F -test for the test.
• H0 : K0 β − m = 0
• HA : K0 β − m 6= 0
• T S : F0 = Q/k
= (K0 β̂ −m)0 [K0 (X0 X)−1 K]−1 (K0 β̂ −m)/k
s2 M S(RESIDUAL)
• RR : F0 ≥ F(α,k,n−p0 )
• P -value: P r(F ≥ F0 ) where F ∼ F(k,n−p0 )
4.4.1 Special Cases of the General Test

Case 1 - Testing a Simple Hypothesis (k = 1)
In this case, K0 (X0 X)−1 K is a scalar, with an inverse that is its recipocal. Also, K0 β̂ − m is
a scalar.
(K0 β̂ − m)2
Q = (K0 β̂ − m)0 [K0 (X0 X)−1 K]−1 (K0 β̂ − m) = 0 0 −1
K (X X) K
!2
(K0 β̂ − m)2 K0 β̂ − m
⇒ F0 = 2 0 0 −1 = = t20
s [K (X X) K] sqrts2 [K0 (X0 X)−1 K]
0
Case 2 - Testing k Specific βjs = 0
In this case, K0 β − m is simply a “subvector” of the vector β, and K0 (X0 X)−1 K is a “sub-
matrix” of (X0 X)−1 . Be careful of row and column labels because of β0 .
Suppose we wish to test that the last q < p elements of β are 0, controlling for the remaining
p − q independent variables:
H0 : βp−q+1 = βp−q+2 = · · · = βp = 0 HA : Not all βi = 0 (i = p − q + 1, . . . , p)
Here, K0 is a q × p0 matrix that can be written as K0 = [0|I], where 0 is a q × p0 − q matrix of 0s

and I is the q × q identity matrix. Then:
 
β̂p−q+1
 β̂p−q+2 
 
K0 β̂ =  ..  K0 (X0 X)−1 K =
 . 
β̂p
 
cp−q+1,p−q+1 cp−q+1,p−q+2 ··· cp−q+1,p
 cp−q+2,p−q+1 cp−q+2,p−q+2 ··· cp−q+2,p 
 
= .. .. .. .. 
 . . . . 
cp,p−q+1 cp,p−q+2 ··· cp,p
where ci,j is the element in the (i + 1)st row and (i + 1)st column of (X0 X)−1 . Then Q is:
 −1 
cp−q+1,p−q+1 cp−q+1,p−q+2 ··· cp−q+1,p β̂p
h i cp−q+2,p−q+1 cp−q+2,p−q+2 ··· cp−q+2,p   β̂p
..   
Q = (K0 β̂)0 [K0 (X0 X)−1 K]−1 (K0 β̂) = β̂p−q+1 β̂p−q+2 . β̂p  .. .. .. ..  
 . . . .  
cp,p−q+1 cp,p−q+2 ··· cp,p
Q/q
⇒ F0 =
s2
Case 3 – Testing a single βj = 0
This is a simplification of case 2, with K0 (X0 X)−1 K being the (j + 1)st element of (X0 X)−1 ,
and K0 β̂ = β̂j .
" #2
(β̂j )2 (β̂j )2 β̂j
Q= ⇒ F0 = 2 = p = t20
cjj s cjj s2 cjj
4.4.2 Computing Q from Differences in Sums of Squares

The sums of squares for a general test can be obtained by fitting various models, and taking
differences in Residual suns of squares.
H 0 : K0 β = m HA : K0 β 6= m
First, the Full Model is fit, that lets all parameters to be free (HA ), and the least squares estimate
is obtained. The residual sum of squares is obtained and labelled SS(RESIDUALFULL ). Under
the full model, with no restriction on the parameters, p0 parameters are estimated and this sum of
squares has n − p0 .
Next, the Reduced Model is fit, that places k ≤ p0 constraints on the parameters (H0 ). Any
remaining parameters are estimated by least squares. The residual sum of squares is obtained and
labelled SS(RESIDUALREDUCED ). Note that since we are forcing certain parameters to take
on specific values SS(RESIDUALREDUCED ) ≥ SS(RESIDUALFULL ), with the equality only
taking place if the estimates from the full model exactly equal the constrained values under H0 .
With the k constraints, only p0 − k parameters are being estimated and the residual sum of squares
has n − (p0 − k) degrees of freedom.
We obtain the sum of squares and degrees of frredom for the test by taking the difference in
the residual sums of squares and in their corresponding degrees’ of freedom:
Q = SS(RESIDUALREDUCED ) − SS(RESIDUALFULL ) df (Q) = (n − (p0 − k)) − (n − p0 ) = k
As before:
(SS(RESIDUAL )−SS(RESIDUAL
REDUCED FULL )
Q/k 0 0
(n−(p −k))−(n−p )
F = 2 =
s SS(RESIDUAL
FULL )
n−p0
Examples of Constraints and the Appropriate Reduced Models
Suppose that Yi = β0 + β1 Xi1 + β2 Xi2 + β3 Xi3 + β4 Xi4 + εi .
H0 : β1 = β2 (k = 1)
∗ ∗
Yi = β0 +β1 Xi1 +β1 Xi2 +β3 Xi3 +β4 Xi4 +εi == β0 +β1 (Xi1 +Xi2 )+β3 Xi3 +β4 Xi4 +εi = β0 +β1 Xi1 +β3 Xi3 +β4 Xi4 +εi Xi1 =X
H0 : β0 = 100 β1 = 5 (k = 2)
Yi = 100+5Xi1+β2 Xi2 +β3 Xi3 +β4 Xi4 +εi ⇒ Yi −100−5Xi1 = β2 Xi2 +β3 Xi3 +β4 Xi4 +εi Yi∗ = β2 Xi2 +β3 Xi3 +β4 Xi4 +ε
Some notes regarding computation of Q:

• Q can always be computed from differences in residual sums of squares.
• When β0 is in the model, and not involved in H0 : K0 β = 0 then we can use Q = SS(MODELFULL )−
SS(MODELREDUCED ).
• When β0 6= 0, is in the reduced model, you cannot use the difference in Regression sums of
squares, since SS(TOTAL UNCORRECTED)differs between the two models.
Best practice is always to use the error sums of squares.
4.4.3 R-Notation to Label Sums of Squares
Many times in practice, we wish to test that a subset of the partial regression coefficients are all
equal to 0. We can write the model sum of squares for a model containing β0 , β1 , . . . , βp as:
R(β0 , β1 , . . . , βp ) = SS(MODEL)
The logic is to include all βi in R(·) that are in the model being fit. Returning to the case of testing
the last q < p regression coefficients are 0:
H0 : βp−q+1 = βp−q+2 = · · · = βp = 0 HA : Not all βi = 0 (i = p − q + 1, . . . , p)
H0 : R(β0 , β1 , . . . , βp−q ) = SS(MODELREDUCED )

H0 : R(β0 , β1 , . . . , βp ) = SS(MODELFULL )
Q = SS(MODELFULL − SS(MODELREDUCED = R(β0 , β1 , . . . , βp ) − R(β0 , β1 , . . . , βp−q ) =
= R(βp−q+1 , . . . , βp |β0 , β1 , . . . , βp−q )
Special cases include:
SS(REGRESSION) = SS(MODEL) − SS(µ) = R(β0 , β1 , . . . , βp ) − R(β0 ) = R(β1 , . . . , βp |β0 )
Partial (TYPE III) Sums of Squares: R(β0 , . . . , βi−1 , βi , βi+1 , . . . , βp )−R(β0 , . . . , βi−1 , βi+1 , . . . , βp ) =
= R(βi |β0 , . . . , βi−1 , βi+1 , . . . , βp )
Sequential (TYPE I) Sums of Squares: R(β0 , . . . , βi−1 , βi ) − R(β0 , . . . , βi−1 ) = R(βi |β0 , . . . , βi−1 )
SS(REGRESSION) = R(β1 , . . . , βp |β0 ) = R(β1 |β0 ) + R(β2 |β1 , β0 ) + · · · + R(βp |β1 , . . . , βp−1 )
The last statement shows that sequential sums of squares (corrected for the mean) sum to the
regression sum of squares. The partial sums of squares do not sum to the regression sum of squares
unless the last p columns of X are mutually pairwise orthogonal, in which case the partial and
sequential sums of squares are identical.
4.5 Univariate and Joint Confidence Regions

In this section confidence regions for the regression parameters are given. See RPD for cool pictures.
Confidence Intervals for Partial Regression Coefficients and Intercept

Under the standard normality, constant variance, and independence assumptions; as well as the
independence of K0 β̂ and SSE(RESIDUAL), we have:
β̂j ∼ N (βj , σ 2 cjj ) where cjj is the (j + 1)st diagonal element of (X0 X)−1
q q
β̂j − βj
⇒ p ∼ t(n−p0 ) ⇒ P r{β̂j − t(α/2,n−p0 ) s2 cjj ≤ βj ≤ β̂j − t(α/2,n−p0 ) s2 cjj } = 1 − α
s2 cjj
q
⇒ (1 − α)100% Confidence Interval for βj : β̂j ± t(α/2,n−p0 ) s2 cjj
Confidence Interval for β0 + β1 X10 + · · · + βp Xp0 = x00 β
By a similar argument, we have a (1 − α)100% confidence interval for the mean at a given
combination of levels of the independent variables, where Ŷ0 = x00 β̂:
q
Ŷ0 ± t(α/2,n−p0 ) s2 x00 (X0 X)−1 x0
Prediction Interval for Future Observation Y0 when X1 = X10 , . . . , Xp = Xp0 (x0 )

For a (1 − α)100% prediction interval for a single outcome (future observation) at a given
combination of levels of the independent variables, where Ŷ0 = x00 β̂:
q
Ŷ0 ± t(α/2,n−p0 ) s2 [1 + x00 (X0 X)−1 x0 ]
Bonferroni’s Method for Simultaneous Confidence Statements
If we want to construct c confidence statements, with simultaneous conficence coefficient 1 − α,

then we can generate the c confidence intervals, each at level (1 − αc ). That is, each confidence
interval is more conservative (wider) than if they had been constructed one-at-a-time.
Joint Confidence Regions for β
From the section on the general linear tests, if we set K0 = Ip0 , we have the following distribu-
tional property:
( β̂ − β)0 [(X 0 X)−1 ]−1 ( β̂ − β)
∼ F(p0 ,n−p0 )
p0 s2
⇒ P r{( β̂ − β)0 (X0 X)( β̂ − β) ≤ p0 s2 F(1−α,p0 ,n−p0 ) } = 1 − α
Values of β in this set constitute a joint (1 − α)100% confidence region for β.
4.6 A Test for Model Fit
A key assumption for the model is that the relation between Y and X is linear (that is E[Y] = X β).
However, the relationship may be nonlinear. S-shaped functions are often seen in biological and
business applications, as well as the general notion of “diminishing marginal returns,” for instance.
A test can be conducted, when replicates are obtained at various combinations of levels of the
independent variables. Suppose we have c unique levels of x in our sample. It’s easiest to consider
the test when there is a single independent variable (but it generalizes straightforwardly). We
obtain the sample size (nj ), mean (Y j ) and variance (s2j ) at each unique level of X (Y is the overall
sample mean for Y ). We obtain the a partition of SS(TOTAL CORRECTED) to test:
H0 : E[Yi ] = β0 + β1 Xi HA : E[Yi ] = µi 6= β0 + β1 Xi
where µi is the mean of all observations at the level Xi and is not linear in Xi . The alternative
can be interpreted as a 1–way ANOVA where the means are not necessarily equal. The partition
is given in Table 12, with the following identities with respect to sums of squares:
SS(TOTAL CORR) = SS(REG) + SS(LF) + SS(PE)
where SS(RESIDUAL) = SS(LF) + SS(PE). Intuitively, these sums of squares and their degrees
of freedom can be written as:
X
n X
c
SS(LF) = (Y (i) − Ŷi )2 = nj (Y j − Ŷ(j) )2 dfLF = c − 2
i=1 j=1
n
X c
X
SS(PE) = (Yi − Y (i) )2 = (nj − 1)s2j dfLF = n − c
i=1 j=1
where Y (i) is the mean for the group of observations at the same level of X as observation i and
Ŷ(j) is the fitted value. The test for goodness-of-fit is conducted as follows:
Source df SS
Pn
Regression (REG) 1 (Ŷi − Y )2
Pc i=1
Lack of Fit (LF) c−2 n (Y − Ŷ(j) )2
Pc j j
j=1
2
Pure Error (PE) n−c j=1 (nj − 1)sj
Table 12: ANOVA for Lack of Fit F -test
• H0 : E[Yi ] = β0 + β1 Xi (Relation is linear)
• HA : E[Yi ] = µi 6= β0 + β1 Xi (Relationship is nonlinear)

M S(LF ) SS(LF )/(c−2)
• T S : F0 = M S(P E) = SS(P E)/(n−c)
• RR : F0 ≥ F(α,c−2,n−c)
Example – Building Costs

A home builder has 5 floor plans: 1000f t2 , 1500, 2000, 2500, and 3000. She knows that the
price to build individual houses varies, but believes that the mean price may be linearly related to
size in this size range. She samples from her files the records of n = 10 houses, and tabulates the
total building cost for each of the houses. She samples ni = 2 homes at each of the c = 5 levels of
X. Consider each of the following models:
Model 1: E[Y ] = 5000 + 10X σ = 500
Model 2: E[Y ] = −12500 + 30X − 0.05X 2 σ = 500

These are shown in Figure 12.
y
50000
40000 Model 2
30000
Model 1
20000
10000
-10000
-20000
1000 1500 2000 2500 3000
size
Figure 4: Models 1 and 2 for the building cost example
Data were generated from each of these models, and the simple linear regression model
was fit. The least squares estimates for model 1 (correctly specified) and for model 2 (incorrectly
fit) were obtained:
Model 1: Ŷ = 5007.22 + 10.0523X
Model 2: Ŷ = 4638.99 + 10.1252X
The observed values, group means, and fitted values from the simple linear regression model are
obtained for both the correct model (1) and the incorrect model (2) in Table 4.6.
The sums of squares for lack of fit are obtained by taking deviations between the group means
(which are estimates of E[Yi ] under HA ) and the fitted values (which are estimates of E[Yi ] under
H0 ).
n
X
Model 1 SS(LF) = nj (Y (i) − Ŷi )2 = 2(15303.96−15059.40)2 +· · ·+2(35496.89−35164.04)2 = 631872.71
i=1
Correct Model (1) Incorrect Model (2)
i Xi Yi Y (i) Ŷi Yi Y (i) Ŷi
1 1000 14836.70 15303.96 15059.49 12557.72 12429.03 14764.23
2 1000 15771.22 15303.96 15059.49 12300.33 12429.03 14764.23
3 1500 20129.51 19925.80 20085.63 21770.87 21045.46 19826.86
4 1500 19722.08 19925.80 20085.63 20320.05 21045.46 19826.86
5 2000 25389.36 25030.88 25111.77 27181.07 27242.41 24889.48
6 2000 24672.40 25030.88 25111.77 27303.75 27242.41 24889.48
7 2500 29988.95 29801.31 30137.90 31139.83 30931.27 29952.10
8 2500 29613.68 29801.31 30137.90 30722.71 30931.27 29952.10
9 3000 35362.12 35496.89 35164.04 33335.17 32799.24 35014.73
10 3000 35631.66 35496.89 35164.04 32263.31 32799.24 35014.73
Table 13: Observed, fitted, and group mean values for the lack of fit test
n
X
Model 2 SS(LF) = nj (Y (i) −Ŷi )2 = 2(12429.03−14764.23)2+· · ·+2(32799.24−35014.73)2 = 36683188.83
i=1
The sum of squares for pure error are obtained by taking deviations between the observed
outcomes and their group means (this is used as an unbiased estimate of σ 2 under HA , after dividing
through by its degrees of freedom).
n
X
Model 1 SS(PE) = (Yi − Y (i) )2 = (14836.70 − 15303.96)2 + · · · + (35631.66 − 35496.89)2 = 883418.93
i=1
n
X
Model 2 SS(PE) = (Yi − Y (i) )2 = (12557.72 − 12429.03)2 + · · · + (32263.31 − 32799.24)2 = 1754525.81
i=1
The F -statistics for testing between H0 (that the linear model is the true model) and HA (that
the true model is not linear) are:
M S(LF ) SS(LF )/(c − 2) 631872.71/(5 − 2)

Model 1 F0 = = = = 1.19
M S(P E) SS(P E)/(n − c) 883418.93/(10 − 5)
M S(LF ) SS(LF )/(c − 2) 36683188.83/(5 − 2)

Model 2 F0 = = = = 34.85
M S(P E) SS(P E)/(n − c) 1754525.81/(10 − 5)
The rejection region for these tests, based on α = 0.05 significance level is:
RR : F0 ≥ F(.05,3,5) = 5.41
Thus, we fail to reject the null hypothesis that the model is correct when the data were generated
from the correct model. Further, we do reject the null hypothesis when data were generated from
the incorrect model.

Notes 1

Uploaded by

Document Informationclick to expand document information

Copyright:

Available Formats

Notes 1

Uploaded by

Document Information

Original Description:

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Notes 1

Uploaded by

Copyright:

Available Formats

1 Simple Linear Regression

Text: RPD, Chapter 1

1.1 Statistical Model

Example 1 – Pharmacodynamics of LSD

X — Mean tissue concentration of LSD among the volunteers.

Table 1: LSD concentrations and math scores – Wagner, et al (1968)

Example 2 – Estimating Cost Functions of a Hosiery Mill

Y — Monthly total production cost (in $1000s).

X — Monthly output (in thousands of dozens produced).

1.1.2 Generating Data from the Model

1. Specify the model parameters: β0 , β1 , σ

4. Obtain random response Yi = β0 + (β1 Xi ) + (σZi ), i = 1, . . . , n.

1.2 Least Squares Estimation

Table 2: Production costs and Output – Dean (1941)

and the estimated (or fitted or prediction) equation (Ŷi ):

Ŷi = β̂0 + β̂1 Xi i = 1, . . . , n (9)

ei = Yi − Ŷi = Yi − (β̂0 + β̂1 Xi ) i = 1, . . . , n (10)

The residuals sum to 0 for this model:

= Yi − {[Y − β̂1 X] + β̂1 Xi ]}

Example 1 – Pharmacodynamics of LSD

1316.6558 − 1519.1430 −202.4872

A plot of the data and regression line are given in Figure 1.

Example 2 – Estimating Cost Function of a Hosiery Mill

113095.80 − 97575.53 15520.27

Figure 1: Regression of math score on LSD concentration (Wagner, et al, 1968)

cost Y=Math Score X=Tissue LSD Concentration

Figure 2: Estimated cost function for hosiery mill (Dean, 1941)

1.3 Analysis of Variance

We start with the total (uncorrected) sum of squares for Y:

Combining equations (11) and (12), we get:

Now summing ei Ŷi over all observations:

Now making use of the following two facts:

We obtain the desired result:

SS(TOTAL UNCORRECTED) = SS(MODEL) + SS(RESIDUAL)

The computational formulas are obtained as follows:

Source of Degrees of Sum of

Table 5: The Analysis of Variance

Example 1 – Pharmacodynamics of LSD

Source of Degrees of Sum of

Table 6: The Analysis of Variance for the LSD/math score data

The mean squares for regression and residuals are as follow:

Source of Degrees of Sum of

The mean squares for regression and residuals are as follow:

1.4 Precision and Distribution of Estimates

Under the simple linear regression model:

E[Yi ] = µi = β 0 + β 1 Xi V ar[Yi ] = σi2 = σ2

From equation (13):

From Equation (14):

where t∗i = t Pn X(X

By independence of the Yi , we get:

Expanding the first term in the exponent in equation (16), we get:

Expanding the second term in the exponent in equation (16), we get:

1.4.2 Distribution of β̂0

Then, we can write β̂0 as two linear functions of the {Yi }:

Putting these together in equation (19):

E[β̂0 ] = E[U − W ] = E[U ] − E[W ] = β0 + β1 X − β1 X = β0

V ar[U − W ] = V ar[U ] + V ar[W ] − 2Cov[U, W ]