overall yield prediction process.

13


where now b, = c6?/(G? + 63), c is a suitably chosen constant, and 6? and
6} are appropriately chosen estimators.

The extension from combining two estimators to combining k estimators
from k normal populations N(8,0?), i = 1,...,k, is discussed by Norwood
and Hinkelmann (1977). Here

(1.8) TS) = ay Ty +... + OTe

with ®; = 67°/ Sf 6;°. They show that Var (T®) <min {Var T;} if each
sample size is greater than 9, or if some sample size is equal to 9, and the
others need be greater than 17.

For the case k = 2 Nair (1980) computes the variance of T* as an infinite
series, and as a function of two parameters, a? and a = no?/n,03. Of
course, it is symmetric and can be restated as a function of 03 and 1/a.

Following the formulation of Cohen and Sackrowitz (1974), Kubokawa
(1987) provides a family of minimax estimators under normalized quadratic
loss functions. Green and Strawderman (1991) also consider quadratic loss
and provide a James-Stein shrinkage estimator. The use of a quadratic loss
function is extended to the multivariate case by Loh (1991), where now we
have normal populations N(#, 24) and N(0, Xz). As in the univariate case,
there are mean estimators 6), #2 and independent covariance matrix estima-
tors 51,55, each having a Wishart distribution, the loss function is

(1.9) L(, 0,51, D2) = (6 — 0)'(Dz + D57)(0 — 8)

The estimator

(1.10) 6 = (Sy! + Sz1)"1(Dy48, + EZ")

is shown to be best linear unbiased.

The model that we here consider is that there are k normal populations
N(6,0?),i=1,...,. The data available are k unbiased estimators T;,..., Ty
of 6. However, the vector T = (Tj,...,7,) has covariance matrix ¥, for

4


3. The Relative Efficiency of @ for k =2 andn=N-1

The result that Var (0) = ao) Var (6) coincides withwhat intuition sug-
gests: when k = 1, Var (8) = Var (9); when k > 1, Var (9) > Var (6), and for
all k, limy-> co Var (8) = Var (0). But the result gives more precise infor-

mation that helps one to assess the efficiency of the Graybill-Deal estimator
for a given sample size.

Consider the case k = 2, N = n—1. If, without loss of generality, we take
ou = min {011,022}, then Var (0) < min (011, 022) when

(3.1) ai < (n= an)?

s a
N—-3 ~ 031022 — 07g

In the special case for which cov (T,,72) = 0, (3.1) is 1/(N — 3) <
o11/022 < 1, which implies that Var (8) < min (011,092) for all N > 5.
Note that this does not contradict the previously quoted result of Graybill
and Deal (1959): their hypothesis allows N, and No, the sample sizes for
the respective constituent estimators, to be unequal; whereas the current
theorem was derived under the assumption that N; = No = N.

Writing a1, = 07?¢22,0 < a < 1, and denoting the correlation between T;
and T» by p, (3.1) can be written as

1g (ene?
N-37 1-9

(3.2)

i

Then it is apparent that if one varies the parameters a and p so that a— p >
0, the sample size N necessary for (3.2) to hold increases without bound.
But this also is intuitive: a — p + 0 is equivalent to 6 + T,. Given a rough
initial estimate for the parameters a and p, one may use (3.2) to obtain
some idea whether the Graybill-Deal estimator dominates the better of the
two constituent estimators for a given sample size.

Taking the special case 04; = 99, (3.2) becomes

1
N-3

i

<i f
T+

>

we
>

©


The estimator @ is unbiased and has variance

(1.17)
~ e'S1[(T — be)'(T — de)|S~*e
Var(@) = EsEp (es-tep
— , &ST1ES te
~ Ss (e'S~te)? *

In the next section we provide a proof of the basic result:

(1.18) Var() = (2= 1

) Var(6).

2. Proof of the Theorem

The Wishart density of S is

(2.1) f(S)=C(kn) |B") 9 |? eap(—5 tr B48), 3 > 0,

where

k .
Clbyn) = 2a rey
i=l

and X= > 0 (that is, is positive definite).
Let Y = E-3SD-2, so that the density of Y is

(2.2) f(Y) =C(kn) |Y [= exp(—3 trY),¥ >0.
With b = D-2e

x vy-2b
(2.3) Var(@) =€ laa .

6


which there is a sample covariance matrix S having a Wishart distribution
W(Z, k,n). Furthermore, 5 and (T;,...7,) are independent.

When © is known, the linear estimator

(1.11) O= wil +...+ welk, wit... twe=1,

with w;,i = 1,...,k, fixed is unbiased. Let w = (wy,...,w,)’ and e =
(i,...,1)!. For the choice

(1.12) w = (e'D"1)/(e'D7e),

@ is also minimum variance. Furthermore,

e'SME(T — Ge) (T ~ be)|Z-te 1

(1.13) Var(6) = Cane = Fyre"

That Var (0) is minimum variance follows from the Cauchy-Schwartz inequal-
ity:

(1.14) (w'Zw)(e'Z7e) > (w'e) = 1

with equality if and only if (1.12) holds. Also,

1.15 eS te)! < minfo?,..., a2},
1 i

which follows from (1.14) with w =e; = (0,...,0,1,0,..., 0).

When © is unknown it is estimated by S, and we consider the estimator

(1.16) 6 = (e'ST)/(e'S7e).

5


(2.6) Var(6) = (1 + E(zi u! Zu) |Var (6).

The density of Z has the form (2.2), which can be written as

(2.7)
f (Zo, 211, 4) = C(k,n) | Zo2 [oe exp(—5tr Zn) temp(—211)(1—w'u)

Again, using orthogonal invariance, the expectation in (2.6) is

(2.8) Elz! Zp! u] = C(k,n) hh,

where

co nee 1 M+ 2, nee
q, | ap *eap(—52u1)dan =T( )2 ean
0

2

bh = | ulu(1 — uu)" du = (k — ptt 2S 44 pr(2*),
wucr

ne 1
I = | (e, Zz3e) | Zoo [pe exp(—strZ22)dZa2.
Z22>0 .

The integral J) can be evaluated using polar coordinates; it is also a Dirichlet
Integral of Type-I, (see Sobel, Uppuluri and Frankowski, 1977). To simplify
notation in Iz let Q = Z22, so that Q is a (k — 1) x (k — 1) matrix having a
Wishart distributon W(,k —1,n). Then Iz = €(Q7")n/C(k —1,n). But
this expectation is known (see e.g. Kshirsagar, 1978, p. 72) so that

I= [(n— Bo a DERE ys

Combining these results we obtain

Var(@) = (1+ Lilels)Var(@)
n-1 4
= an gv

8


COMBINING CORRELATED UNBIASED ESTIMATORS
OF THE MEAN OF A NORMAL DISTRIBUTION

by

T. Keller
National Agricultural Statistics Service
U. S. Department of Agriculture

Ingram Olkin
Department of Statistics
Stanford University

Technical Report No. 2002-5
March 2002

Department of Statistics
Sequoia Hall
STANFORD UNIVERSITY
Stanford, California 94305-4065

http://www-stat.stanford.edu


time the covariance structure and the biases for a set of these key estimates of
yield (for a given combination of commodity, geographic region, and month
of the prediction) are relatively stable, and hence the covariance structure
and biases can be estimated from historical data.

There are (at least) two possible ways to deal with the biases in the
estimators. If one arbitrarily enforces the constraint that w; + ...+w, =1,
then among all such linear combinations, w)T,+...+-w,7}, of the constituent

estimators T,,..., 7%, the linear combination with the smallest mean square
error is:
, eM T
al = ———
(4.1) e'M—e’
where T = (Tj,...,Z,)/ and M is the mean square error matrix. Then, in
analogy with the development in the case for which T;,..., 7, were unbiased,

one actually computes

iV ets
g-eMot

4.2 —,
(42) e'M-e

where M is an estimate of M based on the historical data. The proper-
ties of the estimator in (4.2) are not well understood. The condition that
wi +...+w, = 1 has a very natural statistical interpretation in the unbiased
setting: the condition insures that the composite estimate is also unbiased.
In the more general setting for which the constituent indications may be
biased, the condition that w; +...+ w, = 1 has no such natural statisti-
cal interpretation, and seems to be merely a mathematical convenience. In
particular, the composite estimator in (4.2) is, in general, biased. An al-
ternative procedure is to use the historical data to estimate the biases of
the constituent estimators, and then use these bias estimates to calculate
‘corrected’ estimators which are unbiased, and proceed as discussed in the
first section of the paper. There are several statistical issues that arise in
considering this latter procedure that need to be addressed. However, in the
majority of applications to yield prediction that we have examined thus far,

1



This form of equation (3.1) implies that the sample size for (3.1) to hold
increases without bound as p > 1. Once again, this is intuitive: to say p is
close to 1 means the estimator T) provides essentially the same information
about @ as the estimator T,, and hence the composite estimator cannot be
expected to provide much more information than that provided by T; alone.

4, An Agricultural Application: Predicting Yield

The National Agricultural Statistics Service (NASS), an agency of the
USDA, makes predictions of yield (defined as production per unit area) for
the major US agricultural commodities, both at the level of the primary
region of production for the commodity, and also at the level of the indi-
vidual states that comprise that region. A diverse body of information is
employed in making a yield prediction, and the process of incorporating that
information into a single estimate of yield is a complex task that relies on
the experience and judgment of commodity experts, as well as statistical
methodology. The goal of research into composite estimation at NASS is not
to replace that experience and judgment, but to apply appropriate statisti-
cal methodology to obtain a consistent means of combining key estimators
of yield into a single yield prediction, and thereby provide a useful starting
point for the process of making a monthly yield prediction.

There are two major types of surveys that form the basis for estimates
of final yield. One is a survey that makes various measurements and counts
for plants from randomly located plots in randomly selected fields having the
commodity in question. The same plots are visited several times throughout
the growing season, and the particular counts and measurements which are
made depend on the point in the growing season at which the data is collected
and the development of the crop, as well as the particular commodity in
question. These survey results are used to produce a yield prediction using
biological models of yield. The other type of survey, directly asks randomly
selected producers of the commodity for their best estimate of the yield they
anticipate for their own fields. This survey is also conducted at intervals
throughout the growing season. These survey data are the basis for more
than one estimate of yield, the estimates depending on how the data are
summarized.

At some point after the crop has been harvested, the actual yield is known
with some degree of accuracy. Although yield certainly varies widely from
year to year, it is reasonable to assume that over some ‘moving window’ of

10


Because the density (2.2) is orthogonally invariant, that is, £(G'YG) =
L(Y) for any orthogonal matrix G, a judicious choice of G allows one to put
(2.3) in a more convenient form. Let e; = (1,0,...,0), and choose G so that.
the first row of G is b!/./b/b and the remaining k — 1 rows of G complete an
orthonormal basis for G. Then, by construction, Gb = Vb‘be;. Consequently,
with Z = G'YG, (2.3) becomes

~ dZe 1
Var(@) = [laze ye

Note that b'b = e’/D~'e, and recall that Var(8) = e/D7'e, so that
(2.4) Var(6) = €[

Remark: For any vector a of unit length, and a positive definite matrix
B,a'B’a > (a'Ba)’. Hence (2.4) demonstrates that Var(9) > Var(@) under
the hypothesis that S and T = (T,...,7,)! are independent, but with no
distributional assumptions on S or T.

Now the task of proving the theorem is reduced to computing the expec-
tation on the right side of equation (2.4). Towards that end, partition the
k x k matrix Z and its inverse as

ALL zy 1 Zn ca
ZL= Zr=( 0 5 ,
( 2 Log ) > ( ZH Zo

where Z. and Zog are both (k-—1) x (k-1).

In what follows we make use of well-known relationships between the
blocks of Z and Z~!. (See, for instance, Anderson, 1984.) Employing these
relationships, and that (I — wu’)-! = I+ ce the expression inside the
expectation brackets in (2.4) can be written as:

Ig? 52. SIS.
Ze By tha

°°) Eee)

=l+a,u/Zz)u,

where u = IP a | (Fa and (2.4) becomes:

7


condition, they show that T* is uniformly better than either 7; or T>, where
better means smaller variance.

Seshadri (1974), motivated by balanced incomplete block (BIB) design
considerations, assumes that there is an unbiased estimator 6 of the ratio
b = 0?/(0? + 03), independent of T; and T, . Then the estimator

(1.4) TO = (1—5)T, +8

is unbiased, and var T® < min (varT;, varT2) provided Var 6 < b? and Var
(1-6) < (1—b)?. The key point is that in certain BIB designs there is an
intra-block and inter-block estimator, and also an estimator b.

Zacks (1966) starts with the assumption that the ratio p = 02/0? is
unknown, and creates an estimator

(1.5) TO = (61, +T)/(O+,

where f is independent of T, and Ty. Then T®) is unbiased. The efficiency
of T®) cannot be given in closed form, and Zacks provides graphs of the
efficiency relative to the estimator T@ with p replacing f.

When the sample sizes of the two samples are equal to n, Cohen and
Sackrowitz (1974) discuss estimators of the form

(1.6) T = GT, + dTr,

where a; are functions of sample variances and are chosen with respect to
a squared error loss function normalized by o?. They determine the sample
size n for which T®) is superior to either T, or To.

Because the estimators T; of the mean and s? of the variances are location
and scale estimators, Cohen (1974) considers a location-scale family as a more
general construct than the normal family. Again, the combined estimator is

(1.7) T =hT,+hT,, b+ =1,


Table 1

Predicted yields (weight per area) of commodity Z for state X in month Y.

farmer

year reported yield yield model

1 88.0
2 82.5
3 83.0
4 73.5
5 79.0
6 82.0
7 83.0
8 80.8
9 81.0
10 79.0
i 64.0
12 80.5
13 83.0
14 81.5

biological

Root Mean Square Error:

Farmer Reported Yield
Biological Yield Model
Composite Estimator
Panel of Experts

3.06
3.92
2.68
2.58

87.5
80.0
86.5
79.0
84.5
83.5
79.8
84.0
83.0
79.0
76.0
83.8
87.0
78.5

14

composite
estimate

87.8
81.5
84.2
75.3
81.3
82.5
81.8
81.8
81.7
79.0
68.3
81.6
84.4
80.4

panel
of exports

89.5
82.5
85.8
76.3
83.3
83.8
85.0
81.3
81.8
81.0
67.5
83.0
85.0
82.0

true
yield

87.8
87.3
85.3
76.8
78.3
89.0
82.5
84.0
82.3
80.8
68.3
83.0
85.0
81.8


Seshadri, V. (1963) Constructing uniformly better estimators. Journal of the
American Statistical Association 58 172-175.

Sobel, M., Uppuluri, R., and Frankowski, K. (1977) Selected tables in math-
ematical statistics, Vol. 4: Dirichlet Distribution - Type 1, American Math-
ematical Society.

Tippett, L.H.C. (1931). The method of statistics Williams and Norgate:
London.

Zacks, S. (1966) Unbiased estimation of the common mean of two normal
distributions based on small samples of equal size, Journal of the American
Statistical Association 61 467-476

16


making a prediction requires the combination of estimators. The present
analysis was motivated by a model to predict agricultural yield. However,
the model is generic, and occurs in a variety of contexts. The specifics of the
application are discussed in Section 4.

It is perhaps surprising that the earliest methods for combining estimators
were nonparametric. Fisher (1930) and Tippett(1931) proposed methods for
combining p -values obtained from independent studies. Fisher was moti-
vated by agriculture and Tippett by industrial engineering.

The parametric problem was first posed by Cochran (1937), who was
also motivated by an agricultural problem. For simplicity suppose that we
have two estimators T; and T> of @ from a N (0,07) and (8,03) population,
respectively. The combined estimator

(1.1) T = wil, + wel

with w, = o7?/(o7? + 03°), w. = o9°/(o7" + og”) is unbiased and has
variance

252
102

(1.2) VT)= Bee

S min(oj, 03).

Consequently, the combined estimator dominates either single estimator.

In practice the variances are unknown, and estimates 67,63 independent
of T,,T2, are substituted in w; and we, that is,

(1.3) T* = thTy + tT.

Of course, now the variance of T* is no longer minimum variance, but it
is unbiased.

Cochran’s paper was the genesis for a sequence of papers. We briefly
describe these in chronological order. Graybill and Deal (1959) started with
the Cochran model and assumed that the estimators 6? and 63 are indepen-
dent and that each arises from a sample of size larger than 9. Under this


1.

Combining Correlated Unbiased Estimators of the Mean of a
Normal distribution

T. Keller and I. Olkin

National Agricultural Statistics Service and Stanford University

Abstract

There are many applications in which one seeks to combine multiple
estimators of the same parameter. If the constituent estimators are
unbiased, then the fixed linear combination which is minimum variance
unbiased is well-known, and may be written in terms of the covariance
matrix of the constituent estimators. In general, the covariance

matrix is unknown, and one computes a composite estimate of the
unknown parameter with the covariance matrix replaced by its maximum
likelihood estimator. The efficiency of this composite estimator

relative to the constituent estimators has been investigated in the
special case for which the constituent estimators are uncorrelated.

For the general case in which the estimators are normally distributed
and correlated, we give an explicit expression relating the variance

of the composite estimator computed using the covariance matrix, and
the variance of the composite estimator computed using the maximum
likelihood estimate of the covariance matrix. This result suggests

that the latter composite estimator may be usefulin applications in
which only a moderate sample size is available. Details of one such
application are presented: combining estimates of agricultural yield
obtained from multiple surveys into a single yield prediction.

Introduction

The need to combine estimators from different sources arises in many

fields of application. In agriculture estimates may come from different exper-
imental stations; in the medical sciences there may be multi-sites or multiple
studies; sample surveys may contain subsurveys at different locations. Often


the latter procedure has resulted in a somewhat smaller mean square error
of prediction than that resulting from the former procedure.

Security issues preclude a presentation of an actual NASS data set. Nev-
ertheless, in order to give the reader some idea how the methodology works
in practice, a masked data set has been prepared. In Table 1, the columns
headed ‘farmer reported yield’ and ‘biological yield model’ correspond to 14
years of yield predictions based on surveys of the anticipated yield reported
by producers, and of the crop itself, respectively. These data have already
been corrected for bias, using bias estimates based on historic data. To mimic
the operation of a moving window of time, each of the 14 estimates in the col-
umn labeled ‘composite estimate’ is calculated using the sample covariance
matrix based on the other 13 years. The column headed ‘panel of experts’
gives the yield prediction produced by a group of commodity experts, who
are privy to a wide variety of information. That information includes not
only the two constituent estimators of yield listed in the table, but other
yield estimators, the raw survey data, yield estimates from previous months,
data on weather, growing conditions, and cultural practices, reports on ge-
netic improvements, information on economic inputs to the crop, prices, and
so forth. The column headed ‘true yield’ is the last value of yield for the
combination of year, geographic area and commodity in question which has
been published by the USDA.

Note that the root mean square error for the composite estimator was less
than that for either of the individual estimators of yield, and only slightly
more than the root mean square error for the yield estimate produced by
the panel of commodity experts. This latter observation is particularly re-
markable in view of the bulk of additional information relevant to predicting
yield that is available to the panel. The results displayed in Table 1 general-
ize to all the combinations of crop, month, and geographic region which we
have tried so far: the composite estimator, on average, does about as well,
or a little better, as a predictor of yield than the estimate of yield produced
by the panel of experts. In some instances, the composite estimator does
quite a bit better than any of the competing estimators, and seldom is it the
case that the mean square error for the composite estimator is very much
larger than the mean square error for the prediction based on the estimate of
yield produced by the expert panel. Thus there is evidence that a significant
reduction in the mean square error in predicting yield can be achieved by
computing a Graybill-Deal type composite estimator as an initial step in the

12


References

Anderson, T.W. (1984) An Introduction to Multivariate Statistical Analysis,
John Wiley Sons, New York.

Cochran, W.G. (1937) Problems arising in the analysis of a series of similar
experiments, Supplement to the Journal of the Royal Statistical Society 4
102-118.

Cohen, A. (1976) Combining estimates of location, Journal of the American
Statistical Association 71 172-175.

Cohen, A. and Sackrowitz, H.B. (1974) On estimating the common mean of
two normal distributions. Annals of Statistics 2 1274-1282.

Fisher, R. A. (1932). Statistical methods for research workers (4th ed.) Oliver
and Boyd; London.

Graybill, F.A. and Deal, R.B. (1959) Combining unbiased estimators, Bio-
metrics 15 543-550.

Green, E.J. and Strawderman, W.E. (1991) A James-Stein type estimator for
combining unbiased and possibly biased estimators, Jornal of the American
Statistical Association 86 1001-1006.

Kshirsagar, A. (1978) Multivariate Analysis, Marcel Dekker, Inc.

Kubokawa, T. (1987) Admissible minimax estimation of a common mean of
two normal populations. Annals of Statistics, 15 1245-1256.

Loh, W-L. (1991) Estimating the common mean of two multivariate normal
distributions, Annals of Statistics, 19 297-313.

Nair, K.A. (1980) Variance and distribution of the Graybill-Deal estimator of
the common mean of two normal populations, Annals of Statistics 8 212-216.

Norwood, T.E. and Hinkelmann, K. Jr. (1977) Estimating the common mean
of several normal populations, Annals of Statistics 5 1047-1050.

Raj, D. (1998) Sampling Theory, McGraw-Hill, 16-17.

15