C-1

Appendix C

Generalized least-squares model description and assumptions
C-2

Consider a region with n gaging stations as follows.

At each gaged site, a streamflow characteristic is estimated, such as the logarithm of the 50-year 

peak flow,

, (1)

where ψi is the true (but unknown) log of the 50-year peak and ηi is a random error. If yi is an 

unbiased estimate of ψi, then ηi (sometimes called time sampling error) has a mean of zero and a 

variance that is a function of how many years of data are available for the site and the standard 

deviation of water-year peaks. In addition, there are k basin characteristics, such as log of drainage 

area, that are measured with negligible error.

Assuming that (within the region defined by the basin characteristics at the n stations) ψ is 

approximately linearly related to the basin characteristics (x’s), then the model formulation can be 

written as:

, (2)

where εi is a model error assumed uncorrelated from observation to observation, with mean zero and 

constant variance, γ2. Substituting into equation 1,

. (3)

In matrix notation:

, (4)

where

, (5)

where E[υ]=0, and E[υυT]=Λ. Now the GLS estimator of β is:

. (6)

yi Ψi ηi = +

Ψi β0 β1x1i β2x2i

... βkxki εi = ++ ++ (i=1,2,...,n; n>k)

yi β0 β1x1i β2x2i ... βkxki ηi εi = + + + + ++

Y X = β υ +

Y

y1

y2

...

yn

X=

1 x11 x21 ... xk1

1 x12 x22 ... xk2

...... .........

1 x1n x2n ... xkn

β

β0

β1

...

βk

υ

ε1 η1 +

ε2 η2 +

...

εn ηn +

= ==

b X = ( ) TΛ–1X –1XTΛ–1Y
C-3

The problem with this estimator is that Λ is unknown and must be estimated from the data. In 

OLS, Λ is estimated as σ2

I, which would be a good estimate if all stations in that region had 

approximately the same lengths of record, or if the variance of ηi is small relative to the variance of 

εi at every station in the region.

Because this assumption may be hard to justify, a better estimate of Λ is attempted. Denote this 

estimated covariance matrix Â, and the GLS estimator, b, will be referred to as an Estimated 

Generalized Least Squares (EGLS) estimator.

EGLS Regression

An example illustrates how Â is estimated. Suppose that yi is the log of the 50-year peak 

estimated from mi

 years of record and that the water-year peaks follow a log-Pearson Type III 

(LPIII) distribution at all sites. Further, to minimize notation, assume that the skew coefficient at all 

sites is zero. The elements of Â would be given by:

. (7)

In this equation, K (LPIII standard deviate for zero skewness and 50-year recurrence interval), 

mi

 (record length at station i), mj (record length at station j), and mij (concurrent record length for 

stations i and j) are known, but σi (standard deviation of water-year peaks at station i), ρij (cross 

correlation of water-year peaks at stations i and j), and γ2 (variance of model error) must be 

estimated from the data. Furthermore, we cannot use si

 (the sample estimate of σi) as an estimate of 

σi without introducing bias, and the use of rij (sample cross correlations) for ρij often causes 

numerical problems. Therefore, we estimate σi and ρij as follows.

The standard deviation of water-year peaks, σi

, is estimated from a regional regression of the form:

(8)

By estimating the standard deviations, si, that enter into equation 7 with equation 8, we are assured 

that the rows of the Λ matrix are not correlated with the observed dependent variable Y. This quality 

is necessary for the estimates of β to be unbiased.

The cross correlation coefficient, ρij, is estimated by developing an empirical relation between 

sample cross correlations, rij, and distance between stations of the form:

. (9)

λij

γ2 σ2

i 1 0.5K2 ( ) +

mi

+ ------------------------------------for i j ( ) =

or

ρijσiσjmij 1 0.5K2 ( ) +

mimj

-----------------------------------------------------for i j ( ) ≠



















=

ln si ( ) b0 b1x1i b2x2i ... bkxki = + + ++

rij Θ

dij

αdij + 1 ------------------

=
C-4

Estimating the cross correlations in this manner assures us that the matrix Λ will be positive 

definite. Figure 1 below shows a smooth curve with Θ=.9812 and α=.00412 based on data from 

Illinois. This curve was developed by running the GLSNET program that will be described later.

Figure 1. Relation between cross correlation and distance.

Now the only parameters left to find in the EGLS model are the regression coefficients, b, and 

variance of the model error, γ2

. The model error variance, γ2, and regression coefficients, b, are 

found by iteratively searching for the best non-negative solution to the equation:

. (10)

The GLSNET/AIDE package leads one through the development of equations 8 and 9 in 

preparation for the estimation of the GLS regression coefficients.

Reporting results and errors

The predicted response at ungaged site k with basin characteristics xk =(1, xk,1, xk,2, ..., xk,p) is:

=xkb. (11)

The standard error of the prediction in OLS regression is:

E ( ) Y X– β { } TΛ–1( ) Y X– β = n k – 1 –

yˆ

k
C-5

S( )= {σ2[1+ xk(X’X)-1x’k]}0.5 . (12)

In GLS regression, the standard error of prediction is:

. (13)

The S( ) is a function of x and the computed standard error of a prediction in percent will also be a 

function of x.

Standard Errors in Percent

When a standard error or average prediction error in log units follows a normal distribution, the 

error may be expressed in percent of the predicted value in cubic feet per second (ft3

/s). Denote σ as 

the standard error in log (base 10) units, Scfs as the standard error in ft3/s, and E(q|xk) as the 

predicted value of q, in ft3/s, given xk, and xk =(1, xk,1, xk,2, ..., xk,p) is a vector of basin 

characteristics. The standard error in percent, Spercent is given by:

 (14)

(Aitcheson and Brown, 1957).

Sometimes it is said in OLS that two-thirds of the points lie within one standard error of 

estimate of the regression function. This is true for the log unit standard error of estimate, σ, but it 

generally is not correct for Spercent. This is true because the errors in log space are symmetrically 

distributed under the assumption of normality of the log errors, but the errors in ft3/s are skewed. 

You can, however, calculate +percent and -percent errors with the following formulas:

; and (15)

. (16)

The three formulas above apply not only to the standard error of estimate for an regression, but 

they also apply to the standard error of the model, , in GLS regression, the average prediction 

error, and standard error of a prediction in both OLS and GLS.

Average prediction error (APE)

One overall measure of how good the regression model is for prediction is the average 

prediction error (Hardison, 1971), where the average is taken over prediction sites with X variables 

identical to the observed data. This measure assumes the observed data have been collected at a 

representative set of sites in the region. It is computed as:

yˆ

k

S yˆ

k () γˆ2 xkX'Λ

ˆ –1

X–1x' = + k

yˆ

k

Spercent 100

Scfs

Eqxk ( )

------------------ 100 e5.302σ2

= = ( ) – 1

Splus 100 10σ = ( ) – 1

Sminus 100 10–σ = ( ) – 1

γ

ˆ
C-6

. (17)

The first term in the brackets on the right side of equation 17 represents an estimate of the average 

squared model error for the n sites and the second term inside the brackets is an estimate of the 

average squared error due to estimating true model parameters from a sample of data.

Prediction interval

Users of the regression model are probably more interested in a measure of error in a particular 

prediction rather than an average prediction. A good measure of the error of a particular prediction is 

the confidence interval of a prediction, or prediction interval. Let x0 represent the usual row vector 

of basin characteristics at a prediction site. As usual x0 is augmented by a 1 as the first element. The 

predicted value is . A 100(1-α) prediction interval would be:

 , (18)

where

, (19)

where tα/2, n-p′

 is the critical value from a t-distribution for n-p′ degrees of freedom.

If a log transform had been made so that y0 = log10(q0), then the prediction interval would be:

. (20)

SELECTED REFERENCES

Acreman, M.C., and Wiltshire, S.E., 1987, Identification of regions for regional flood frequency 

analysis [abs,], EOS, v. 68, no. 44.

Aitchison, J., and Brown, J.A.C., 1957, The Log-Normal Distribution: London, Cambridge 

University Press.

Hardison, C.H., 1971, Prediction error of regression estimates of streamflow characteristics at 

ungaged sites: U.S. Geological Survey Professional Paper 750-C, p. 228-236.

APE γ

ˆ2

i

n

------ xi X'Λ

ˆ –1

( ) X –1x'

i

n -------------------------------------

i = 1

n

+ ∑

i = 1

n

∑

   

  1 2/

=

yˆ

0 x0 = b

yˆ

0 – T y0 yˆ0 ≤ ≤ + T

T tα

2

---, n p – ' γ

ˆ2

0 x0 X'Λ

ˆ –1

( ) X –1x'

0 = ( ) +

10yˆ 0 – T q0 10yˆ 0 + T ≤ ≤