NAG Library Routine Document

G02JDF

+− Contents

1 Purpose

2 Specification

3 Description

4 References

5 Parameters

6 Error Indicators and Warnings

7 Accuracy

8 Further Comments

+− 9 Example

9.1 Program Text

9.2 Program Data

9.3 Program Results

1 Purpose

G02JDF fits a multi-level linear mixed effects regression model using restricted maximum likelihood (REML). Prior to calling G02JDF the initialization routine G02JCF must be called.

2 Specification

SUBROUTINE G02JDF (

LVPR, VPR, NVPR, GAMMA, EFFN, RNKX, NCOV, LNLIKE, LB, ID, LDID, B, SE, CZZ, LDCZZ, CXX, LDCXX, CXZ, LDCXZ, RCOMM, ICOMM, IOPT, LIOPT, ROPT, LROPT, IFAIL)

INTEGER	LVPR, VPR(LVPR), NVPR, EFFN, RNKX, NCOV, LB, ID(LDID,LB), LDID, LDCZZ, LDCXX, LDCXZ, ICOMM(*), IOPT(LIOPT), LIOPT, LROPT, IFAIL
REAL (KIND=nag_wp)	GAMMA(NVPR+1), LNLIKE, B(LB), SE(LB), CZZ(LDCZZ,), CXX(LDCXX,), CXZ(LDCXZ,), RCOMM(), ROPT(LROPT)

3 Description

G02JDF fits a model of the form:

y = X β + Z ν + ε

where	$y$ is a vector of $n$ observations on the dependent variable,
	$X$ is a known $n$ by $p$ design matrix for the fixed independent variables,
	$β$ is a vector of length $p$ of unknown fixed effects,
	$Z$ is a known $n$ by $q$ design matrix for the random independent variables,
	$ν$ is a vector of length $q$ of unknown random effects,
and	$ε$ is a vector of length $n$ of unknown random errors.

Both

ν

and

ε

are assumed to have a Gaussian distribution with expectation zero and variance/covariance matrix defined by

Var [\begin{matrix} ν \\ ε \end{matrix}] = [\begin{matrix} G & 0 \\ 0 & R \end{matrix}]

where

R = σ_{R}^{2} I

I

is the

n \times n

identity matrix and

G

is a diagonal matrix. It is assumed that the random variables,

Z

, can be subdivided into

g \leq q

groups with each group being identically distributed with expectation zero and variance

σ_{i}^{2}

. The diagonal elements of matrix

G

therefore take one of the values

\{σ_{i}^{2} : i = 1, 2, \dots, g\}

, depending on which group the associated random variable belongs to.

The model therefore contains three sets of unknowns: the fixed effects

β

, the random effects

ν

and a vector of

g + 1

variance components

γ

, where

γ = \{σ_{1}^{2}, σ_{2}^{2}, \dots, σ_{g - 1}^{2}, σ_{g}^{2}, σ_{R}^{2}\}

. Rather than working directly with

γ

, G02JDF uses an iterative process to estimate

γ^{*} = \{σ_{1}^{2} / σ_{R}^{2}, σ_{2}^{2} / σ_{R}^{2}, \dots, σ_{g - 1}^{2} / σ_{R}^{2}, σ_{g}^{2} / σ_{R}^{2}, 1\}

. Due to the iterative nature of the estimation a set of initial values,

γ_{0}

, for

γ^{*}

is required. G02JDF allows these initial values either to be supplied by you or calculated from the data using the minimum variance quadratic unbiased estimators (MIVQUE0) suggested by Rao (1972).

G02JDF fits the model by maximizing the restricted log-likelihood function:

- 2 l_{R} = \log (|V|) + (n - p) \log (r^{T} V^{- 1} r) + \log |X^{T} V^{- 1} X| + (n - p) (1 + \log (2 π / (n - p)))

where

V = Z G Z^{T} + R, r = y - X b and b = {(X^{T} V^{- 1} X)}^{- 1} X^{T} V^{- 1} y .

Once the final estimates for

γ^{*}

have been obtained, the value of

σ_{R}^{2}

is given by

σ_{R}^{2} = (r^{T} V^{- 1} r) / (n - p) .

Case weights,

W_{c}

, can be incorporated into the model by replacing

X^{T} X

and

Z^{T} Z

with

X^{T} W_{c} X

and

Z^{T} W_{c} Z

respectively, for a diagonal weight matrix

W_{c}

The log-likelihood,

l_{R}

, is calculated using the sweep algorithm detailed in Wolfinger et al. (1994).

4 References

Goodnight J H (1979) A tutorial on the SWEEP operator The American Statistician 33(3) 149–158

Harville D A (1977) Maximum likelihood approaches to variance component estimation and to related problems JASA 72 320–340

Rao C R (1972) Estimation of variance and covariance components in a linear model J. Am. Stat. Assoc. 67 112–115

Stroup W W (1989) Predictable functions and prediction space in the mixed model procedure Applications of Mixed Models in Agriculture and Related Disciplines Southern Cooperative Series Bulletin No. 343 39–48

Wolfinger R, Tobias R and Sall J (1994) Computing Gaussian likelihoods and their derivatives for general linear mixed models SIAM Sci. Statist. Comput. 15 1294–1310

5 Parameters

Note: prior to calling G02JDF the initialization routine G02JCF must be called, therefore this documention should be read in conjunction with the document for G02JCF.

In particular some parameter names and conventions described in that document are also relevant here, but their definition has not been repeated. Specifically, RNDM, WEIGHT, N, NFF, NRF, NLSV, LEVELS, FIXED, DAT, LICOMM and LRCOMM should be interpreted identically in both routines.

1: LVPR – INTEGERInput

On entry: the sum of the number of random parameters and the random intercept flags specified in the call to G02JCF.

Constraint:

LVPR = \sum_{i} RNDM (1, i) + RNDM (2, i)

2: VPR(LVPR) – INTEGER arrayInput

On entry: a vector of flags indicating the mapping between the random variables specified in RNDM and the variance components,

σ_{i}^{2}

. See Section 8 for more details.

Constraint:

1 \leq VPR (i) \leq NVPR

, for

i = 1, 2, \dots, LVPR

3: NVPR – INTEGERInput

On entry:

g

, the number of variance components being estimated (excluding the overall variance,

σ_{R}^{2}

Constraint:

1 \leq NVPR \leq LVPR

4: GAMMA( $NVPR + 1$ ) – REAL (KIND=nag_wp) arrayInput/Output

On entry: holds the initial values of the variance components,

γ_{0}

, with

GAMMA (i)

the initial value for

σ_{i}^{2} / σ_{R}^{2}

, for

i = 1, 2, \dots, NVPR

GAMMA (1) = - 1.0

, the remaining elements of GAMMA are ignored and the initial values for the variance components are estimated from the data using MIVQUE0.

On exit:

GAMMA (i)

, for

i = 1, 2, \dots, NVPR

, holds the final estimate of

σ_{i}^{2}

and

GAMMA (NVPR + 1)

holds the final estimate for

σ_{R}^{2}

Constraint:

GAMMA (1) = - 1.0 ​ or ​ GAMMA (i) \geq 0.0

, for

i = 1, 2, \dots, g

5: EFFN – INTEGEROutput

On exit: effective number of observations. If there are no weights (i.e.,

WEIGHT ='U'

), or all weights are nonzero, then

EFFN = N

6: RNKX – INTEGEROutput

On exit: the rank of the design matrix,

X

, for the fixed effects.

7: NCOV – INTEGEROutput

On exit: number of variance components not estimated to be zero. If none of the variance components are estimated to be zero, then

NCOV = NVPR

8: LNLIKE – REAL (KIND=nag_wp)Output

On exit:

- 2 l_{R} (\hat{γ})

where

l_{R}

is the log of the restricted maximum likelihood calculated at

\hat{γ}

, the estimated variance components returned in GAMMA.

9: LB – INTEGERInput

On entry: the dimension of the arrays B and SE and the second dimension of the array ID as declared in the (sub)program from which G02JDF is called.

Constraint:

LB \geq NFF + NRF \times NLSV

10: ID(LDID,LB) – INTEGER arrayOutput

On exit: an array describing the parameter estimates returned in B. The first

NLSV \times NRF

columns of ID describe the parameter estimates for the random effects and the last NFF columns the parameter estimates for the fixed effects.

A print routine for decoding the parameter estimates given in B using information from ID is supplied with the example program for this routine.

For fixed effects:

for l=NRF×NLSV+1 ,…, NRF×NLSV+NFF
- if $B (l)$ contains the parameter estimate for the intercept then
  $ID (1, l) = ID (2, l) = ID (3, l) = 0;$
- if $B (l)$ contains the parameter estimate for the $i$ th level of the $j$ th fixed variable, that is the vector of values held in the $k$ th column of DAT when $FIXED (j + 2) = k$ then
  $\begin{array}{l} ID (1, l) = 0, \\ ID (2, l) = j, \\ ID (3, l) = i; \end{array}$
- if the $j$ th variable is continuous or binary, that is $LEVELS (FIXED (j + 2)) = 1$ , then $ID (3, l) = 0$ ;
- any remaining rows of the $l$ th column of ID are set to $0$ .

For random effects:

let
- $N_{R_{b}}$ denote the number of random variables in the $b$ th random statement, that is $N_{R_{b}} = RNDM (1, b)$ ;
- $R_{j b}$ denote the $j$ th random variable from the $b$ th random statement, that is the vector of values held in the $k$ th column of DAT when $RNDM (2 + j, b) = k$ ;
- $N_{S_{b}}$ denote the number of subject variables in the $b$ th random statement, that is $N_{S_{b}} = RNDM (3 + N_{R_{b}}, b)$ ;
- $S_{j b}$ denote the $j$ th subject variable from the $b$ th random statement, that is the vector of values held in the $k$ th column of DAT when $RNDM (3 + N_{R_{b}} + j, b) = k$ ;
- $L (S_{j b})$ denote the number of levels for $S_{j b}$ , that is $L (S_{j b}) = LEVELS (RNDM (3 + N_{R_{b}} + j, b))$ ;
then
- for $l = 1, 2, \dots NRF \times NLSV$ , if $B (l)$ contains the parameter estimate for the $i$ th level of $R_{j b}$ when $S_{k b} = s_{k}$ , for $k = 1, 2, \dots, N_{S_{b}}$ and $1 \leq s_{k} \leq L (S_{j b})$ , i.e., $s_{k}$ is a valid value for the $k$ th subject variable, then
  $\begin{array}{l} ID (1, l) = b, \\ ID (2, l) = j, \\ ID (3, l) = i, \\ ID (3 + k, l) = s_{k}, k = 1, 2, \dots, N_{S_{b}}; \end{array}$
- if the parameter being estimated is for the intercept then $ID (2, l) = ID (3, l) = 0$ ;
- if the $j$ th variable is continuous, or binary, that is $L (S_{j b}) = 1$ , then $ID (3, l) = 0$ ;
- the remaining rows of the $l$ th column of ID are set to $0$ .

In some situations, certain combinations of variables are never observed. In such circumstances all elements of the

l

th row of ID are set to

- 999

11: LDID – INTEGERInput

On entry: the first dimension of the array ID as declared in the (sub)program from which G02JDF is called.

Constraint:

LDID \geq 3 + \max_{j} (RNDM (3 + RNDM (1, j), j))

, i.e.,

3 +

maximum number of subject variables (see G02JCF).

12: B(LB) – REAL (KIND=nag_wp) arrayOutput

On exit: the parameter estimates, with the first

NRF \times NLSV

elements of

B

containing the parameter estimates for the random effects,

ν

, and the remaining NFF elements containing the parameter estimates for the fixed effects,

β

. The order of these estimates are described by the ID parameter.

13: SE(LB) – REAL (KIND=nag_wp) arrayOutput

On exit: the standard errors of the parameter estimates given in B.

14: CZZ(LDCZZ, $*$ ) – REAL (KIND=nag_wp) arrayOutput

Note: the second dimension of the array CZZ must be at least

NRF \times NLSV

(see G02JCF).

On exit: if

NLSV = 1

, then CZZ holds the lower triangular portion of the matrix

(1 / σ^{2}) (Z^{T} {\hat{R}}^{- 1} Z + {\hat{G}}^{- 1})

, where

\hat{R}

and

\hat{G}

are the estimates of

R

and

G

respectively. If

NLSV > 1

then CZZ holds this matrix in compressed form, with the first NRF columns holding the part of the matrix corresponding to the first level of the overall subject variable, the next NRF columns the part corresponding to the second level of the overall subject variable etc.

15: LDCZZ – INTEGERInput

On entry: the first dimension of the array CZZ as declared in the (sub)program from which G02JDF is called.

Constraint:

LDCZZ \geq NRF

16: CXX(LDCXX, $*$ ) – REAL (KIND=nag_wp) arrayOutput

Note: the second dimension of the array CXX must be at least

NFF

(see G02JCF).

On exit: CXX holds the lower triangular portion of the matrix

(1 / σ^{2}) X^{T} {\hat{V}}^{- 1} X

, where

\hat{V}

is the estimated value of

V

17: LDCXX – INTEGERInput

On entry: the first dimension of the array CXX as declared in the (sub)program from which G02JDF is called.

Constraint:

LDCXX \geq NFF

18: CXZ(LDCXZ, $*$ ) – REAL (KIND=nag_wp) arrayOutput

Note: the second dimension of the array CXZ must be at least

NLSV \times NRF

(see G02JCF).

On exit: if

NLSV = 1

, then CXZ holds the matrix

(1 / σ^{2}) (X^{T} {\hat{V}}^{- 1} Z) \hat{G}

, where

\hat{V}

and

\hat{G}

are the estimates of

V

and

G

respectively. If

NLSV > 1

then CXZ holds this matrix in compressed form, with the first NRF columns holding the part of the matrix corresponding to the first level of the overall subject variable, the next NRF columns the part corresponding to the second level of the overall subject variable etc.

19: LDCXZ – INTEGERInput

On entry: the first dimension of the array CXZ as declared in the (sub)program from which G02JDF is called.

Constraint:

LDCXZ \geq NFF

20: RCOMM( $*$ ) – REAL (KIND=nag_wp) arrayCommunication Array

Note: the dimension of the array RCOMM must be at least

LRCOMM

(see G02JCF).

On entry: communication array initialized by a call to G02JCF.

21: ICOMM( $*$ ) – INTEGER arrayCommunication Array

Note: the dimension of the array ICOMM must be at least

LICOMM

(see G02JCF).

On entry: communication array initialized by a call to G02JCF.

22: IOPT(LIOPT) – INTEGER arrayInput

On entry: optional parameters passed to the optimization routine.

By default G02JDF fits the specified model using a modified Newton optimization algorithm as implemented in E04LBF. In some cases, where the calculation of the derivatives is computationally expensive it may be more efficient to use a sequential QP algorithm. The sequential QP algorithm as implemented in E04UCA can be chosen by setting

IOPT (5) = 1

. If

LIOPT < 4

IOPT (5) \neq 1

then E04LBF will be used.

Different optional parameters are available depending on the optimization routine used. In all cases, using a value of

- 1

will cause the default value to be used. In addition only the first LIOPT values of IOPT are used, so for example, if only the first element of IOPT needs changing and default values for all other optional parameters are sufficient LIOPT can be set to

1

E04LBF is being used

$i$	Description	Equivalent E04LBF parameter	Default Value
$1$	Number of iterations	MAXCAL	$1000$
$2$	Unit number for monitoring information	n/a	As returned by X04ABF
$3$	Print optional parameters ( $1 =$ print)	n/a	$- 1$ (no printing performed)
$4$	Frequency that monitoring information is printed	IPRINT	$- 1$
$5$	Optimizer used	n/a	n/a

If requested, monitoring information is displayed in a similar format to that given by E04LBF.

E04UCA is being used

$i$	Description	Equivalent E04UCA parameter	Default Value
$1$	Number of iterations	Major Iteration Limit	$\max (50, 3 \times NVPR)$
$2$	Unit number for monitoring information	n/a	As returned by X04ABF
$3$	Print optional parameters ( $1 =$ print, otherwise no print)	List/Nolist	$- 1$ (no printing performed)
$4$	Frequency that monitoring information is printed	Major Print Level	$0$
$5$	Optimizer used	n/a	n/a
$6$	Number of minor iterations	Minor Iteration Limit	$\max (50, 3 \times NVPR)$
$7$	Frequency that additional monitoring information is printed	Minor Print Level	$0$

23: LIOPT – INTEGERInput

On entry: length of the options array IOPT. If

LIOPT \leq 0

then IOPT is not referenced and default values are used for all optional parameters.

24: ROPT(LROPT) – REAL (KIND=nag_wp) arrayInput

On entry: optional parameters passed to the optimization routine.

Different optional parameters are available depending on the optimization routine used. In all cases, using a value of

- 1.0

will cause the default value to be used. In addition only the first LROPT values of ROPT are used, so for example, if only the first element of ROPT needs changing and default values for all other optional parameters are sufficient LROPT can be set to

1

E04LBF is being used

$i$	Description	Equivalent E04LBF parameter	Default Value
$1$	Sweep tolerance	n/a	$\max (\sqrt{eps}, \sqrt{eps} \times \max_{i} ({zz}_{i i}))$
$2$	Accuracy of linear minimizations	ETA	$0.9$
$3$	Accuracy to which solution is required	XTOL	$0.0$
$4$	Initial distance from solution	STEPMX	$100000.0$

E04UCA is being used

$i$	Description	Equivalent E04UCA parameter	Default Value
$1$	Sweep tolerance	n/a	$\max (\sqrt{eps}, \sqrt{eps} \times \max_{i} ({zz}_{i i}))$
$2$	Lower bound for $γ^{*}$	n/a	$eps / 100$
$3$	Upper bound for $γ^{*}$	n/a	$10^{20}$
$4$	Line search tolerance	Line Search Tolerance	$0.9$
$5$	Optimality tolerance	Optimality Tolerance	${eps}^{0.72}$

where

eps

is the machine precision returned by X02AJF and

{zz}_{i i}

denotes the

i

diagonal element of

Z^{T} Z

25: LROPT – INTEGERInput

On entry: length of the options array ROPT. If

LROPT \leq 0

then ROPT is not referenced and default values are used for all optional parameters.

26: IFAIL – INTEGERInput/Output

On entry: IFAIL must be set to

0

- 1 ​ or ​ 1

. If you are unfamiliar with this parameter you should refer to Section 3.3 in the Essential Introduction for details.

For environments where it might be inappropriate to halt program execution when an error is detected, the value

- 1 ​ or ​ 1

is recommended. If the output of error messages is undesirable, then the value

1

is recommended. Otherwise, if you are not familiar with this parameter, the recommended value is

0

. When the value $- 1 or 1$ is used it is essential to test the value of IFAIL on exit.

On exit:

IFAIL = 0

unless the routine detects an error or a warning has been flagged (see Section 6).

6 Error Indicators and Warnings

If on entry

IFAIL = 0

- 1

, explanatory error messages are output on the current error message unit (as defined by X04AAF).

Errors or warnings detected by the routine:

$IFAIL = 1$: On entry, LVPR is too small.

$IFAIL = 2$: On entry, $VPR (i) < 1$ or $> NVPR$ .

$IFAIL = 3$

On entry,	$NVPR < 1$ ,
or	$NVPR > LVPR$ .

$IFAIL = 4$: On entry, $GAMMA (1) \neq - 1.0$ and $GAMMA (i) < 0$ .

$IFAIL = 9$: On entry, LB is too small.

$IFAIL = 11$: On entry, LDID is too small.

$IFAIL = 15$: On entry, LDCZZ is too small.

$IFAIL = 17$: On entry, LDCXX is too small.

$IFAIL = 19$: On entry, LDCXZ is too small.

$IFAIL = 21$: On entry, ICOMM has not been initialized.

$IFAIL = 101$: Optimal solution found, but requested accuracy not achieved.

$IFAIL = 102$: Too many major iterations.

$IFAIL = 103$: Current point cannot be improved upon.

$IFAIL = 104$: At least one negative estimate for gamma was obtained. All negative estimates have been set to zero.

7 Accuracy

Not applicable.

8 Further Comments

The parameter VPR gives the mapping between the random variables and the variance components. In most cases

VPR (i) = i

, for

i = 1, 2, \dots, \sum_{i} RNDM (1, i) + RNDM (2, i)

. However, in some cases it might be necessary to associate more than one random variable with a single variance component, for example, when the columns of DAT hold dummy variables.

Consider a dataset with three variables:

DAT = (\begin{matrix} 1 & 1 & 3.6 \\ 2 & 1 & 4.5 \\ 3 & 1 & 1.1 \\ 1 & 2 & 8.3 \\ 2 & 2 & 7.2 \\ 3 & 2 & 6.1 \end{matrix})

where the first column corresponds to a categorical variable with three levels, the next to a categorical variable with two levels and the last column to a continuous variable. So in a call to G02JCF

LEVELS = (\begin{matrix} 3 & 2 & 1 \end{matrix})

also assume a model with no fixed effects, no random intercept, no nesting and all three variables being included as random effects, then

\begin{matrix} FIXED = (\begin{matrix} 0 & 0 \end{matrix}); \\ RNDM = {(\begin{matrix} 3 & 0 & 1 & 2 & 3 \end{matrix})}^{T} . \end{matrix}

Each of the three columns in DAT therefore correspond to a single variable and hence there are three variance components, one for each random variable included in the model, so

VPR = (\begin{matrix} 1 & 2 & 3 \end{matrix}) .

This is the recommended way of supplying the data to G02JDF, however it is possible to reformat the above dataset by replacing each of the categorical variables with a series of dummy variables, one for each level. The dataset then becomes

DAT = (\begin{matrix} 1 & 0 & 0 & 1 & 0 & 3.6 \\ 0 & 1 & 0 & 1 & 0 & 4.5 \\ 0 & 0 & 1 & 1 & 0 & 1.1 \\ 1 & 0 & 0 & 0 & 1 & 8.3 \\ 0 & 1 & 0 & 0 & 1 & 7.2 \\ 0 & 0 & 1 & 0 & 1 & 6.1 \end{matrix})

where each column only has one level

LEVELS = (\begin{matrix} 1 & 1 & 1 & 1 & 1 & 1 \end{matrix}) .

Again a model with no fixed effects, no random intercept, no nesting and all variables being included as random effects is required, so

\begin{matrix} FIXED = (\begin{matrix} 0 & 0 \end{matrix}); \\ RNDM = {(\begin{matrix} 6 & 0 & 1 & 2 & 3 & 4 & 5 & 6 \end{matrix})}^{T} . \end{matrix}

With the data entered in this manner, the first three columns of DAT correspond to a single variable (the first column of the original dataset) as do the next two columns (the second column of the original dataset). Therefore VPR must reflect this

VPR = (\begin{matrix} 1 & 1 & 1 & 2 & 2 & 3 \end{matrix}) .

In most situations it is more efficient to supply the data to G02JCF in terms of categorical variables rather than transform them into dummy variables.

9 Example

This example fits a random effects model with three levels of nesting to a simulated dataset with

90

observations and

12

variables.

NAG Library Routine DocumentG02JDF

+− Contents

1 Purpose

2 Specification

3 Description

4 References

5 Parameters

6 Error Indicators and Warnings

7 Accuracy

8 Further Comments

9 Example

9.1 Program Text

9.2 Program Data

9.3 Program Results

NAG Library Routine Document

G02JDF