SoM [Xy.¢ Klay — Xe KKB]

éL ,

= xo (Ky: — XX) + =
op Rol ¥ tu>t
aon 1 — 2p"

Sat gen g (vine — 8

aL,

So KK

t
S xo" (XX, XL)
tt<u

"B78

(4)

= NCaystry F esoen)
bata

ae Nat
rattatt
Now
n
coe aay
atta! 2

= Ky. + Kye — (XX, + XUKB
ift-fu

pe

x’

—XiX.p ift=u

CALCULATING
THE MAXIMUM
LIKELIHOOD ESTIMATES
The likelihood equations are nonlinear and
have to be solved by iterative procedures. The
initial estimates of the model parameters that
we use are the ordinary least squares esti-
mates 8, = (X'X)“X’y, and the estimated
variances and covariances based on A.

Bun = (Vu — BEXyn — BLK. +
8.X!.X,8,) /N.

The program uses three procedures for
solving the likelihood equations; (1) Newton-

2

(3)

y. +p’/X'.X.B)

Raphson; (2) Fisher scoring; and (3) the
iterated Aitken procedure. Different iterative
procedures are needed because one of the pro-
cedures may not converge to a solution; and
the different procedures converge at different
rates. The method of determining which pro-
cedure to use at each iteration is: if the change
in L is less than a specified value at the pre-
ceeding iteration, the Newton-Raphson pro-
cedure is used; if the change in L is greater
than the specified value, or if the Newton-
Raphson procedure does not work, the Fisher
scoring method is used. Finally, the iterated
Aitken method is used if both Newton-Raph-
son and Fisher scoring procedures fail to
work. These two procedures may not work at
a particular iteration because of an ill-
conditioned matrix of second derivatives, an
ill-conditioned estimate of %, or a decrease in
the value of the likelihood function when
evaluated at the new estimates.

The three methods of solving the likelihood
equations are based on calculating a vector
of increments (A) to the vector of parameter
estimates (6). The iterated parameter esti-
mates are 6, Ot Aw

The vector of increments A for the Newton-
Raphson procedure is

aL ebay pa
BB BaD
A> ® (5)
eh aL a
emp anee,,., Lad.

The first matrix on the right side of (5) is
the inverse of the matrix of second deriva-
tives of the likelihood function and the sec-
ond matrix is the vector of first derivatives.
The notation [ ],., indicates that the ele-
ments of the vector of first derivatives (3)


BERR Re Rok RR Ee gee EE EE EAE EER EE EE REE EE EE EE EERE AE SEE RE
*#* INSTRUCTIONS FOR CHANGING ARRAY DIMENSIUNS. THIS IS REQUIRED ONLY
IN

MAIN ROUTINE.

1. DETERMINE P AND Q.
2. FROM THESE VALUES CALCULATE PPsPPP,QQ,MOL,MD2 USING ARRAY SIZE

DEFINITIONS BELOW.
3. CHANGE DIMENSION SIZE OF ARRAYS TO APPKUPRIATE VALUES.

4. CHANGE MXP TO Py MXQ TO Qy AND MXMD1 TO MDL.
SRR Sa tot a AOI RO HOR SOR ORR I RE SOE TOO A AR RE EE

RRR BRI dea HR EE RAO RE RR LOR SR OO ARE OR

eee ARRAY SIZE DEFINITIONS
PP = P#(P+1)I/2
PPP= PPH(PP4LIs/2
QQ = QelQeL)/2
MDL= Q+PP MAX = MXMDL
MD2= MDL*®(MO1L4+1)/2
HRI RR HER EE EEE RE EEE EERE ER EERE ROA EEE EEE EE EER EEE REESE
*##* ARRAY DIMENSIONS, MUST BE AT LEAST THE SIZE BELOW.
REAL*8& DEREV2(MD2),SIGMACPP),YPY(PP)sSIGINVEPP) PLS(PP),
DER TV 1(0MO1L} p PLBB(QQ), PLSS( PPP), DELTA(MDL) ,Y(P) yX(PyQ)y
REAL*8& XPX(PP4Q,Q) ,XPYIP,QeP) BETALQ) pPLBSLQy PP), PLB(Q)
ERE EER EE SHEERS CER ERER EEE EERE ES EEK EEE EERE E EEK EES EERE EE REE RER ORK EE EEE

ida dhaddladiadladlasisindintidindndindiatndnininideihinshthdedniednhinthinthinindhiniindaininbibinte ett eee
*** THERE ARE 4 PLACES WHERE FATAL MATRIX INVERSIONS MIGHT OCCUR.
THESE ARE NUMBERED WITH MIER. A MESSAGE IS PRINTED WITH THIS
VALUE TO TELL THE USER THE MATRIX IN WHICH THE INVERSION FAILED.
IF MEER = Ly THE INITIAL X*X COULD NOT BE INVERTED
25 THE INITIAL ESTIMATE OF SIGMA COULD NOT BE INVERTED
3e X*S XyTHE MATRIX OF 2ND DERIVS. WITH RESPECT TO BETA
4y X*S Xy AT LAST ITERATION
REREE EE ER EEREE ERE EE EE EERE ESS REREEE RENEE EEE EREEHEEREER EERE REREEERE EE EEE


to B and &. The elements of these matrices
can be used as estimates of the variances and
covariances of B and Z.

othe value wk can be used to set up likeli-
hood ratio tests to determine whether some
of the independent variables can be dropped
from the model. If Lu is the value of the log
likelihood function under the full model, and
L,, is the value obtained with the restricted
model (deletion of a portion of 8 of length
qi), then ~2 (L,, — Lu) is asymptotically dis-
tributed as y° with q, degrees of freedom.

‘The estimates of B and & are both printed
and punched. The punched cards can be read
as starting values for further iterations (see
program control information).

PROGRAM APPLICATION

‘The program was initially developed to ob-
tain solutions to the forest growth and yield
functions as described by Sullivan and Clut-
ter (1972) when there are more than two
repeated measurements per plot. The model
and procedure presented here are more gen-
eral, however, With proper construction of
the data, the program can obtain solutions to
the usual multivariate linear regression model
as well as solutions to some systems of equa-
tions, where the dependent variables are
combinations of repeated measurements and
measurements on several correlated variates.
Note that the X, can be combinations of ex.
perimental design variables and covariates
and need not be the same for the different y,.

A sample problem is presented in Appendix
II for 3 repeated measurements on each of 22
fixed plots. We use Clutter’s (1963) growth
and yield model to describe plot. volume and
basal area. The equations can be written as

E(InV.,) = Bo + & Sit Br An? + B, IB,
E(InB,,) = AyAut MB, + a.(1~ A) Aut)

+ a(1~AvAn)S,

where V,, = cubic foot volume per acre on
the i" plot at time t,
S. = site index of i plot (in feet),
A. = stand age of i" plot at time t,
and —_B,, = basal area per acre of the i plot

at time t (in square feet).

We consider stand ages (A,) and initial
basal area (B,) as design variables. Note that
the basal area equation for time 1 is an iden-
tity, while the coefficient for A\,/A,, InB,, in
the basal area equation for the second and
third remeasurement periods is always 1. To
obtain the computing form for the data, we
delete the first basal area equation and sub-
tract A.A,‘ InB, from both sides of the tt
hasal area equation for t = 2, 3.

Appendix IT contains a partial listing of the
sample data and the control cards required
to run the program with the data. The first
four fields of the sample data cards contain
identification information only, and are not
used by the program. The values are plot
number, measurement period, equation type,
and variate number.

The data for the five equations for each
plot are together. The first three cards have
the data for the three volume equations, the
last two cards have the data for the basal
area equations for time 2 and 3. The data on
the data cards are the values of the variables
defined as shown in Table 1.

Table 1.—Data layout

Y

xX, XX, xX, X: Xs
inV, A:* InB, 0.0 0.0
Inv, A>? InB, 0.0 0.0
nV, A;?— InB, 0.0 0.0
nB.* 0.05 0.0 I-AA Gg ~AAS)S
InB, 80° 00 I-AAS (HAAS
where InB, — AA. InB,.
natin Tn

4


and the elements of the matrix of second
derivatives (4) are evaluated at 6... 0L/0Z"
denotes a vector whose tu" element is ¢L/do'",
&L/cBe=" is a matrix whose tu" column is
the vector L/P 0", and ?L/2E-0E" is a
matrix whose element in row (r,s) and column
(tu) is PL/Po"te", for lc r<s<i,...,
pRilstsu=1,...,p

The vector of increments for the Fisher
scoring procedure is

eh 1 aL
apoB op
A=- (6)
eL aL
aE, LEB ds

The first step of the iterated Aitken proce-
duce is to calculate the increment to the vec-
tor of regression coefficients (8) which is

. eL ft éL
a6.) = [zal GL.

The iterated estimate of the vector of regres-
sion coefficients is 6, = &. + A(B). It can
be shown that , is the generalized least.
squares estimator of @ evaluated at &;.,. The
iterated variance and covariance estimates are

Jo = (yye — BY (Xyn +
Xd) + BX! XB) /N.

The iterations are terminated when L, —
Li. < « where « is a small predetermined
value, and L, is the logarithm of the likeli-
hood function evaluated at 9;.

The major computational problem is the
matrix inversion required for the calculation
of A,, For the Newton-Raphson procedure,
the matrix of second derivatives evaluated at
6,, may be singular. We also found that,
when the residuals are highly correlated, the
matrix of second derivatives with respect to
the elements of E* is ill-conditioned. This
second problem was alleviated when we no-
ticed that the inverse of #°L/eZ*0Z* could be
calculated directly from the estimates of o'!
calculated at the previous iteration.

Denote @*L/d E~dE" as A. We found that
the element in row (7,8) and column (t,u) of
the inverse of A is

ACEC) oe

rat

eget ot oh),

n

Thus, the Fisher scoring procedure requires
only one matrix inversion, the qxq matrix of
second derivatives with respect to p. If the
matrix of second derivatives needed for the
Newton-Raphson procedure is negative defi-
nite, it’s inverse can be calculated with one
more inversion of a matrix of order q by using
definitions of the inverse of a partitioned
matrix.

The point at which an iteration should be
started with the Newton-Raphson procedure
has not been determined. In their procedure
for variance component estimation, Jennrich
and Sampson (1976) use the scoring algo-
rithm whenever the change in L is greater
than 1.0. On the data sets we have worked
with, we have found that the Newton-
Raphson procedure converges more rapidly
than the other two procedures when the full
matrix of second derivatives can be inverted.
The point. at which this occurs seems a char-
acteristic of each data set.

When a solution is obtained by a particular
algorithm, but the change in L is negative,
the program switches to a partial stepping
procedure similar to that suggested by Jenn-
rich and Sampson (1976). At each step, A is
divided by 4 and new estimates are obtained.
The partial stepping is terminated at the first
step which increases L.

PROGRAM INFORMATION

The program arrays are dimensioned for a
maximum p of 10 and a maximum q of 12. If
larger problems are encountered, the size of
the appropriate arrays should be changed. In-
structions for doing this are found in the
comments section at the beginning of the
program. The control card instructions are
also located here, as well as in Appendix I.

When a matrix cannot be inverted, the pro-
gram terminates with an error message. Ap-
pendix I contains an explanation of these
messages.

The program output includes the estimates
of B, &, |&|, the value of L evaluated at these
estimates, and the inverse of the negative of
the matrix of second derivatives with respect


APPENDIX |

EKG ERK EEE EE ERE EE EER EER REE RE ERAS EE EE REE EEE ERE EEE EERE SERRE EE EER
#ex CONTROL CARDS

le
2e

be
Te

TITLE CARO FORMAT = (10A8)

JOB CONTROL CARD

COLS. 1-10 DELT = END ITERATIONS INDICATOR. WHEN THE
PROPORTIONATE CHANGE OF THE LUG LIKELIHOOD
FUNCTION IS LESS THAN DELTy ITERATIONS
TERMINATE. IF 0 IS READ, DELT IS SET TO 1.E-1O

11-20 SNW = START NEWTGN-RAPHSON PRUCEOURE WHEN THE

PROPORTIONATE CHANGE OF THE LOG LIKELTHUOD
FUNCTION IS LESS THAN *SNW*, DEFAULT = .1 .
IF THE PROPURTIONATE CHANGE IS GREATER THAN
*SNW*, THE FISHER SCORING PROCEOURE IS USED.

21-25 N = NO. OF EXPERIMENTAL UNITS (NO. UF PLOTS).
26-30 P = NO. OF REPEATED MEASUREMENTS OR DEPENDENT
VARTATES MEASURED ON EACH PLUT.
31-35 Q = NO. OF INOEPENDENT VARIABLES
36-40 NEITR = MAX. NO. ITERATIONS. (IF 0, NITR=20).
41-45 IPR = PRINT OPTIUN
= 0 PRINT RESULTS AT LAST ITERATION ONLY.
= 1 PRINT LOG LIKELIHOOD FUNCTION, OIFFERENCE,
AND VALUE OF QUADRATIC FURM EACH ITERATION.
= 2 PRINT EVERYTHING BUT NEG. OF INVERSE OF
MATRIX OF SECOND DERIVATIVES EACH ITERATION.
46-50 F = FILE ON WHICH DATA IS LUCATED (DEFAULT=5).
51-55 NF = NU. FORMAT CARDS (MAX = 3).
56-60 INT = INTERCEPT INDICATOR: IF INT NUT EQUAL Oy THE
NO. OF INDEPENDENT VARS. IS INCREASED BY L AND
A COL. OF L'S IS ADDED TO INDEP. VARS.
61-65 NC = NO.UF ADDITIONAL CUMMENTS CARDS (CARD TYPE 3)
66-70 INP) = INPUT INDICATOR. IF INP eee

O INITIAL ESTIMATES OF 3ETA AND SIGMA ARE
CALCULATED FROM DATA
1 READ STARTING VALUES UF BETA AND SIGMA
{CARD TYPE 6 AND 7) PUNCHED ON PREVIGUS RUN.
COMMENTS CARDS - OPTIONAL, READ UNLY IF ¥*NC# UN CARD 2 IS «GT. 0%
FURMAT = (10A8) - VARIABLE NAMES CAN BE READ HERES ILeEe
YL = LN({VOL) Y2 = LN(BA) XL = L/AGE eee ETC.
FURMAT OF DATA DESCRIBING THE Pe(Utl) VARIABLES FOR A PLOT.
(ALL VARIABLES ARE REAL)
DATA (1¥ FOLLOWED BY Q X*S3 THIS IS REPEATED P TIMES FOR
THE P CORRELATED OBSERVATIONS ON 1 PLOT. THERE ARE N
OF THESE OBSERVATIONS)

READ CARDS 6 AND 7 GNLY IF INP ON CARD 2 IS EQUAL TO lL.
INITIAL ESTIMATES OF BETAy FORMAT = (5016.8)

INITIAL ESTIMATES OF SIGMA (UPPER TRIANGULAR PORTIUN STORED
COLUMNWISE - [.E. Chel) Cle2) (242) (193) (293) 393d aoe
FORMAT = (5016.8).

SOR RR ORR RRO EMO TO RE ER RE RR EERE EH EA ERE RE RE EEE

7


ODC [853:525.1]--051.5

A COMPUTER PROGRAM

for the Maximum Likelihood Estimator
of the General Multivariate Linear Model
with Correlated Errors

FOREST SERVICE GENERAL TECHNICAL REPORT NE-51
1979

FOREST SERVICE, U.S. DEPARTMENT OF AGRICULTURE
NORTHEASTERN FOREST EXPERIMENT STATION
370 REED ROAD, BROOMALL, PA. 19008


1 0.7220395ne 36
De1999 2678 55,

2599699069

3 3 2.29975 5290 va errcuerers
& dlsterecrysg Lalitibeciongs rth 3oe" papartya7cer4
. ease 9.22772224 etisiteaces
F  GrERIK ees “aang ss ALL Zevuase de <M DeRWHIE ME TTTSATOGLD
Tae vatue ur rar auapsatic Fok» bas, vagvo

Leo of tae Ce Cine re

aN ON

oe eam Draseetes

a

ITeeAT ICN =

bats sie

apayit

PROCEDUSE USED At

ugncut LeMATE JiFEE

sien tue aids sc = meENce

U6 oF He Seley

ues ALISL EMC Co functyeM «

ance

TSS ASLUEOK IME auacoaric foes

Sete guocene tare beeteeenes

REGRESSED Curt F

CEENIS (CTA)
O.297S397IDE DE DLTITATST ANNO

e471 958620001 Je97a EL LED OE

deta etr yo2eat

ores7
SSENMTE UE coy
eC erereteer)
2 G.lehgaza0-99
3 dlasseurar1>

ware ce

2038 OMI?
Deasngecona a>
= TE O24 150-94
URS RAIA

devs 72a9i9-a2
asenreroe
aleeaansnoen2

BLAND eo 74
cr ereereery

dal Tey nec- 32

SINGS 20-95 94249567 aem92

LTS (PARAM NES Cusaht ty
2e2 Da O8
Clie Wem o1
dlegtreune de

CF guar teers sts et
=O, 294912 0-U5 “C,991 1400-07
salle wigs 51 2el9z4aa79-37
Berea cers

Ww. SETA) WEL GD ¢ SUSHAC TAI) y (Lede? EseL 2)
Dln2e43d-97 2.6 1584ge-LU -O.54Ts6v0- IL
3ldusse9-92 24 ane le-22 -916940999-92

0.260085 01

1
Tplasteer 33

THE NeRW

tar CAT EIHCH +

etna beseacar

SETACedy eta a stomata Gh
POrreySterrrrrnereererrisrrirtttt

re cryatayey,

vec cr eresy,

TESaranoa Oana, fey acnsay suena niensys sesemitceds, etscesageor 3.41478 10-97
DEtw reat aT BLT bred Sipe tancan o 3M 91 IF 9,54 RZEC-IT OF1159120-05 ~9.5137230-97

OI6OACHONAIS <BIPEANOT OLagHey Fie

We pereaaruany er se sernace wart

UVESSC OL SEL SAE DS OF Feet viurs ob ocean tk var bars inthe SesefoT Te wera
Lost rmae a te
4 2 deracar rye
“ AIVANPDAMn= ROTHER ewer? ess sen vocmae
. Mas2go Tens ye OTLMIMAL IO va ENeizetn oe gees res reeae
BORO RIED SHAT TNE EMO IN MNS REP TIONG SM buNnsee de <I TMEsHTErD2 965022020096
tovense nh Card # tatu ah earcery AE AWNINGS jut oespecr to stowact yu) Crags sat ey
selecdesreesena cease cetgness eet coames eseranepeeecerenanee pereereecceressei te
earreserenys
2 Ciawatenee sy arta syste

yee ode ede
6 close ear cae
CCCP aes
© give Tal 6t= de
+ SORES SSSI?
no -byrsigur nese

9 css aam eee
Die 92 en

to giau ant see
ipo T9y Pa 8
ay anlaaz2 ara
Dir keeg cm Jet
12 OL2TdMow aL =e
Olea 2a 6
ta Diet Wve ee
alasteaae oo?
le sere ror 2o4c- 27
Tolease aw ae = 37
1s aetterveD~3F

CaMEvcrrmnns

Covarrance mate TAT
1 o.tyap0agceat
+ RNR ayaa

3 ory reaun oreo

5

2019 ITP 12D
trcanene

ORISSA EO

deAP tau eton de
DEN 42e SEI I>
AH TaD ECs
Sag Tes e998
Wate sere
Nigaateg ys
237082 220-97
361 10192 D>
solepzorato= 27
SUVs 1-08
DELL 81 390= 9%
Wo4227 819+ Ie

(ORATLATICR FOR

a pvaycone dt
Perrone)
wale roe 10809
De rean tere C

ian trates
SOL@aeANISUn Oeartaaceor

aewamdyeje Senne se
Manette le Wiesnigaua a?
MnetaI98 haw alten oe
TIADEWS Bae Me Naame pede aT

Veh kay Wage

oer Lote

detaeevory

Tah eeusce oe

sredteges r=a6
less 2e12>:0F

cree
Nugsae sou,
Dees P92 19-97
=e getas5- 2c
Pol iseane ate
Veescstp-ar
DIE 94d~ 6
Wie Ne

Del PpAI EOLA Ie

deverecstce ma
Siitesisic36
wlsgaset gray
eas restteceo7
SP enget at
seep dent tain de
Perera
Dn65107 IHS

detageeacpear
= 2eu2eat prego

Yetean0c8
dSsorive sete Poneae

M2 M92 Oe 30

10

D.1802263E-96
Druk 19 798-96
D132 Lane nee
32820483 76-99

928 36356096
2.5E 5 TSLREIE
Dentasro2t-3¢
3.4408 1L-96
3.276009 16-97
124 tee at-26
0193 78395E-37
9257622500-97

25251901 30-96,
922661673E-07

3.40 90000691

9230 495980-96
Dossyas 9-16
0213728 590-96

942595478

0.1L 85992

70,23542390036
Bea S4eNGb0-IE
2.8 786806097

~0.92602080-07
O21 72a2630-36

071 7835019496
9114736580-96

DebireseTo-us
30124229 8erd0

Ose aren ten 26
O21320240-05
-0.43083820-94
0.81516520-96
0.9931687G-97
0.16783 640-96

Peveerertays
9281590390-36

Tepoeect

94295 4480-95
Disotai1y-o2 -vizsesstOmae

0.897 Tuod-97
91105 2220-vs

3413902 935-90
Bee 978 S4u-ve
= 244 3446 5-96
202 S9ae 300~D6
964884) 30-95
-9.77339T9u-9T
9436829900-07

9690998 Tu-aT


APPENDIX II

SAMPLE PROBLEM, GROwTH AN YLELO ustng eLurea's EOvA TOMS,

a

STMULTANEGUS EQUATIONS

io

: Porras

i 2 = nivour at Fo

Ye 2 TN Gaar ! Cwreal 1 #aGe1/AGERD =

XBT SOMCGA Jcttaalieacenjaccai ye nt = [Mrencerts 32 = S106, 1006
DAUEED ach RRR LITRES 22" GALATI est

YU G2yV3 2 FOL XD OB KG) ord

BO vaayS = FUK5 9 XOP
y RL owe ay * KS xb

Fan ceseeasoesers sone sesen: ferrtrererrrotc

Ueegr 105 02n5 01481008) reeeeeres

PLT TSUts iio 6e.9 ovorite 4.26818 960 2.0
L212 719505 120 6620 Slows siseatt | OF0 20
13.13 Tooter 11a Gerd gldite  Sleeall 920 0:3
122% oretads af3 “60a Ola on Diiasae | 7261539
1325 ogee of0 oro Ole oro Drives? 2Lc79572
211i 617689 110 0r9 Oloe00 5:3 ato
2212 toess ilo eos Slores2 oro 020
2313 730038 Ilo coro Ololnes dro 029
2224 or61085 919 Or9 3l9 OLiMtt Breese?
25.25 1590391 910 oso oso Shaves 10.3448
aria 2.02864 9.96575 96 o.0
a2ie Diowe2  Aligo1s 909 323
aaa bledder Ubebss ole oro
B24 ato 033 Olizs3s  Treeee7
goes 200 ono Slaceoa 12204082
geiti Dloass Sheate 359 ato
Brie ¢ orosize se5317 950 do
22313 Fikes 1155725 oloiset 3¥ere2 30 or9
22224 or6s1y 919 G19 O10 on Diere 722680
22425 1132957 019 013 oto a Slits wali reer

SAMPLE PROBLEM, GROMTM AND YIELO USING CLUTTER ES CQUATIONS
PROGRAM CONTROL INF LPMATION

NO. PLOTS =
NOs CORRELATED VARIATES =
NO, TNOEDENDENT VARIABLES =
MAKe NOs ITERATIONS ©
ORLRT OPTION =

INPUT FILE =

NO. FORMAT CaRUS
INTERCERT INDICATNE =
DELTA (END TTERATIONS INDICATCRE » 2.1926-05
SNK INEWTON RAPHSOW CUTOFF POENT? # 0.1904 60C

15 WINER PROBLEM LDENTIELCSTICN

porppeerrrrrests eocrecnesecencasenseese
SEMULTANEOUS EOLATICYS tariance RAPS
wh = UN(yOL) AT TEM ty V2 = UNUVCUD AT TiME

V3 CMEVOLY AT TEAC 32 Fae MLULGAZI = UNCSATT SAGER ZAG E 2S
Mis URtEBCEPTE x2 0° SITE INDEK:
AS 8 CAISRUTI NS A6 e ETAT SATII
VaevS = FUSy ROD see

y moo as ” xe

AOK FLD. 5 2F5e1 AFL 909T
aogusten 6o6 4e am

ctissisicees).
; DIL29LIDEMEII Oe LZISEHYIO

2 RUSE) Shadeniateo, sasase2seo-v1

2 gears eheaasie i? <203laeseeorse — astuer2ur-o2

$ rapzpenesiody pr Uagsarieat caabstesiocss  STHASSSER2 y spscconeeoe
tue nereaniawwr CE ChISS PROBED ATIC LY = YORE) vgtsssgarvent

un veag.CSTTSATES, 36 MESSER CLEC Len eer)

Or3edsssoo IL HULIVELZILD-Ol ~P.7S1STIICVIZ II s¢aracvede reagzesr tae y4151S4296U-91
Nos SC MysTeCeTS 2 9

INL THAL ESTIAATE UE oe SAS

area)

ive enenarnsnt G6 SUMED,



The first control card is a general title card.
The second is the job control card. On this
card, the first two fields are blank, indicating
the program assumes the default values for
DELT and SNW. The third field contains the
number of plots (22) and the fourth field
contains the number of variates (or equa-
tions). There are six independent variables;
the maximum number of iterations is set at
30; the print option is set to zero to print the
results at the last iteration only; the data is
located on file 5; there is one format card; no
column of 1’s is added to the data; there are
11 additional comments cards; and the blank
in the last field indicates that the initial esti-
mates are calculated from the data.

Following the job control card are the 11
additional comments cards. Next is the data
format card followed by the data itself.

The data are a portion of those collected
for a growth and yield study of managed hard-
wood stands being made by the Northeastern
Forest Experiment Station. The program was
written in Fortran for an IBM 370/168.

A deck and listing of the computer program
described in this publication is available on
request with the understanding that the U.S.
Department of Agriculture cannot assure its

2 U.S. GOVERNMENT PRINTING OFFICE: 1979—603-0116

accuracy, completeness, reliability, or suit-
ability for any other purpose than that re-
ported. The recipient may not assert any
proprietary rights thereto nor represent it to
anyone as other than a Government-produced
computer program. For cost information,
please write the authors at the: Northeastern
Forest Experiment Station, 370 Reed Road,
Broomall, Pa. 19008.

LITERATURE CITED

Clutter, J. L.
1963. Compatible growth and yield models for
loblolly pine. For. Sci. 9(3): 355-371.

Jennrich, R. L, and P. F, Sampson.
1976. Newton-Raphson and related algorithms
for maximum likelihood variance component
estimation. Technometrics 18(1): 11-17.

Seegrist, D. W., and S. L. Arner.
1978. Statistical analysis of linear growth and
yield models with correlated observations from
permanent plots measured at fixed intervals. In
Growth models for long term forecasting of
timber yields. Publ. FWS-1-78. (Jéran Fries,
Harold E. Burkhart, and Timothy A. Max,
eds.) Sch. For. Wildl. Res. Va. Polytech. Inst.
State Univ., Blacksburg. p. 209-223.

Sullivan, A. D., and J. L. Clutter.
1972. A simultaneous growth and yield model
for loblolly pine. For. Sci. 18: 76-86.


Ix THIS REPORT we describe a computer
program that obtains maximum likelihood
estimates of the parameters for a general mul-
tivariate linear model with two or more corre-
lated observations. The variates can be a set
of different types of variates, repeated meas-
urements on the same variate, or combina-
tions of repeated measurements and different
variates.

We present, as an example, the growth and
yield model described by Clutter (1963). This
model uses recursive simultaneous equations
to describe plot volume and basal area where
the values are measured repeatedly on the
same plot.

The general multivariate linear model for p
correlated measurements on each of n plots
can be written

y= Xpte q)
where y’ = (y'1, ¥'y+ ++ Wo)
X= (X4, Xy,..-, Xb)
and oy, = (Yiy-+ +) Yin)
Xi = (Rye Rin)
My = (Xiu. Kiya)

The element y,,; is the observed value of
variate Y, on plot j and x;;, is the observed
value of the k" independent variable associ-
ated with Y, on plot j.

We assume that the vectors of observations
for the n plots are mutually independent and
that the p observations for each plot have
the covariance matrix £. Therefore, the co-
variance matrix for y can be written as Y =
= X I, where X represents the Kronecker
product. The order of £ is p, and I is the
identity matrix of order n.

Due to correlation among the observations,
the ordinary least squares procedure cannot
be used to solve the system (1) for estimates
of B. Therefore, we developed a computer pro-
gram to obtain maximum likelihood estimates
for Band Z.

The program maximizes the logarithm of
the likelihood function of y;

L = —np/2 In(an) + n/2/8"|-1/2Q (2)

where Q = (y— XB)’ ¥ (y- XB)

and |£+| denotes the determinant of matrix
=,

For the maximization of L, we require the
first and second derivatives of L with respect
to B and E+. Denote the tu element of E-* as
o', The required derivatives are presented in
a form convenient for computation. Their
derivation can be found in Seegrist and Arner
(1978).

The first derivatives can be written


A COMPUTER PROGRAM

for the Maximum Likelihood Estimator
of the General Multivariate Linear Model
with Correlated Errors

Abstract

This report describes a computer program that obtains maximum
likelihood estimates of the parameters for a general linear model
with correlated observations. The variates can be a set of different
types of variates, repeated measurement on the same variate, or com-
binations of repeated measurements and different variates. A forest
growth and yield model describing repeated measurement of plot
volume and basal area is presented as an example.



The Authors

STANFORD 1. ARNER joined the Northeastern Forest Experi-
ment Station in 1966 and has served as a computer programmer and
biological statistician with the Station's Biometrics Group, located
in Broomall, PA. He received a B. S. degree in forest science in 1966
from the Pennsylvania State University, and a M. S. degree in 1971
from the School of Forest Resources at Pennsylvania State Univer-
sity.

DONALD W. SEEGRIST joined the Northeastern Forest Experi-
ment Station in 1968 and currently is Leader of the Station’s Bio-
metrics Group. He received a B.S, degree in zoology in 1953 from
George Washington University, and M.S. and Ph, D. degrees in zo-
ology in 1959 and 1965, respectively, from the University of Califor-
nia, Berkeley. He joined the Forest Service in 1963 as a biological
statistician with the Pacific Southwest Forest and Range Experiment
Station at Berkeley.

MANUSCRIPT RECEIVED FOR PUBLICATION 30 NOVEMBER 1978