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

2 Specification

Fortran Interface

Subroutine g02jhf (

hlmm, nvpr, gamma, effn, rnkx, ncov, lnlike, lb, b, se, czz, ldczz, cxx, ldcxx, cxz, ldcxz, rcomm, icomm, ifail)

Integer, Intent (In)	::	nvpr, lb, ldczz, ldcxx, ldcxz
Integer, Intent (Inout)	::	icomm(*), ifail
Integer, Intent (Out)	::	effn, rnkx, ncov
Real (Kind=nag_wp), Intent (Inout)	::	gamma(nvpr+1), czz(ldczz,rnlsvnrf), cxx(ldcxx,fnlsvnff), cxz(ldcxz,fnlsvnff), rcomm()
Real (Kind=nag_wp), Intent (Out)	::	lnlike, b(lb), se(lb)
Type (c_ptr), Intent (In)	::	hlmm

C Header Interface

#include <nag.h>

void

g02jhf_ (void **hlmm, const Integer *nvpr, double gamma[], Integer *effn, Integer *rnkx, Integer *ncov, double *lnlike, const Integer *lb, double b[], double se[], double czz[], const Integer *ldczz, double cxx[], const Integer *ldcxx, double cxz[], const Integer *ldcxz, double rcomm[], Integer icomm[], Integer *ifail)

The routine may be called by the names g02jhf or nagf_correg_lmm_fit.

3 Description

g02jhf 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

γ

, g02jhf 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. g02jhf 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).

g02jhf 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)))

or the log-likelihood function:

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

where

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

By default the restricted log-likelihood function is used, the log-likelihood function can be chosen through the optional parameter Solver as detailed in the documentation for g02jff.

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 Arguments

Note: prior to calling g02jhf the initialization routine g02jff must be called, therefore this documentation should be read in conjunction with the document for g02jff. In particular some argument names and conventions described in that document are also relevant here, but their definition has not been repeated. Specifically, hlmm should be interpreted identically in both routines.

1: $hlmm$ – Type (c_ptr) Input

On entry: a G22 handle to the internal data structure containing a description of the required model as returned in hlmm by g02jff.

2: $nvpr$ – Integer Input

On entry:

g

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

σ_{R}^{2}

This should be set to the value of nvpr returned by g02jff.

3: $gamma (nvpr + 1)$ – Real (Kind=nag_wp) array Input/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}

Labels for the variance components can be obtained using g22ydf as demonstrated in Section 10.

Constraint:

gamma (1) = - 1.0

gamma (i) \geq 0.0

, for

i = 1, 2, \dots, g

4: $effn$ – Integer Output

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

weight ='U'

), or all weights are nonzero,

effn = n

5: $rnkx$ – Integer Output

On exit: the rank of the design matrix,

X

, for the fixed effects.

6: $ncov$ – Integer Output

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

ncov = nvpr

7: $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.

8: $lb$ – Integer Input

On entry: the dimension of the arrays b and se as declared in the (sub)program from which g02jhf is called.

Constraint:

lb \geq fnlsv \times nff + rnlsv \times nrf

9: $b (lb)$ – Real (Kind=nag_wp) array Output

On exit: the parameter estimates, with the first

nrf \times rnlsv

elements of

b

containing the parameter estimates for the random effects,

ν

, and the remaining

nff \times fnlsv

elements containing the parameter estimates for the fixed effects,

β

Labels for the parameter estimates can be obtained using g22ydf as demonstrated in Section 10.

10: $se (lb)$ – Real (Kind=nag_wp) array Output

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

11: $czz (ldczz, rnlsv \times nrf)$ – Real (Kind=nag_wp) array Output

On exit: if

rnlsv = 1

, 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.

rnlsv > 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 random subject variable, the next nrf columns holding the part of the matrix corresponding to the second level of the overall random subject variable etc.

ldczz = 0

, czz is not referenced.

12: $ldczz$ – Integer Input

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

Constraints:

$ldczz = 0$ or $ldczz \geq nrf$ ;
if $ldczz = 0$ , $ldcxx = ldcxz = 0$ .

13: $cxx (ldcxx, fnlsv \times nff)$ – Real (Kind=nag_wp) array Output

On exit: if

fnlsv = 1

, 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

fnlsv > 1

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

ldcxx = 0

, cxx is not referenced.

14: $ldcxx$ – Integer Input

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

Constraint:

ldcxx = 0

ldcxx \geq nff

15: $cxz (ldcxz, fnlsv \times nff)$ – Real (Kind=nag_wp) array Output

On exit: 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.

ldcxz = 0

, cxz is not referenced.

16: $ldcxz$ – Integer Input

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

Constraint:

ldcxz = 0

ldcxz \geq fnlsv \times nff

17: $rcomm (*)$ – Real (Kind=nag_wp) array Communication Array

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

18: $icomm (*)$ – Integer array Communication Array

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

19: $ifail$ – Integer Input/Output

On entry: ifail must be set to

0

- 1

1

to set behaviour on detection of an error; these values have no effect when no error is detected.

A value of

0

causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of

- 1

means that an error message is printed while a value of

1

means that it is not.

If halting is not appropriate, the value

- 1

1

is recommended. If message printing is undesirable, then the value

1

is recommended. Otherwise, the value

0

is recommended. 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 = 11$: hlmm has not been initialized or is corrupt.

$ifail = 12$: hlmm is not a G22 handle as generated by g02jff.

$ifail = 21$: On entry, $nvpr = 〈value〉$ .
Constraint: $nvpr \geq 〈value〉$ .

$ifail = 31$: On entry, $gamma (〈value〉) = 〈value〉$ .
Constraint: $gamma (1) = - 1.0$ or $gamma (i) \geq 0.0$ .

$ifail = 81$: On entry, $lb = 〈value〉$ .
Constraint: $lb \geq 〈value〉$ .

$ifail = 121$: On entry, $ldczz = 〈value〉$ .
Constraint: $ldczz \geq 〈value〉$ .

$ifail = 141$: On entry, $ldcxx = 〈value〉$ .
Constraint: $ldcxx \geq 〈value〉$ .

$ifail = 161$: On entry, $ldcxz = 〈value〉$ .
Constraint: $ldcxz \geq 〈value〉$ .

$ifail = 171$: On entry, the communication arrays, icomm and rcomm, have not been initialized correctly.

$ifail = 1001$: Optimal solution found, but requested accuracy not achieved.

$ifail = 1002$: Too many major iterations.

$ifail = 1003$: Current point cannot be improved upon.

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

$ifail = - 99$: An unexpected error has been triggered by this routine. Please contact NAG.
See Section 7 in the Introduction to the NAG Library FL Interface for further information.

$ifail = - 399$: Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library FL Interface for further information.

$ifail = - 999$: Dynamic memory allocation failed.
See Section 9 in the Introduction to the NAG Library FL Interface for further information.

7 Accuracy

Not applicable.

8 Parallelism and Performance

g02jhf is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.

g02jhf makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.

Please consult the X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this routine. Please also consult the Users' Note for your implementation for any additional implementation-specific information.

9 Further Comments

If g02jgf has been called, then rather than use the values of fnlsv, nff, rnlsv, nrf and nvpr as returned by g02jff they should be obtained by calling g02zlf. See Section 9 in g02jgf for more details.

10 Example

This example fits a random effects model with three random submodels and two fixed effects to a simulated dataset with

90

observations and

12

variables. The model is fit using restricted maximum likelihood (REML). Custom labels for the parameter estimates and variance components are then constructed from information returned by g22ydf. See Section 10 in g02jff for an example that uses the standard parameter labels directly.

g02jh: FL CL CPP AD

NAG FL Interfaceg02jhf (lmm_​fit)

▸▿ Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

NAG FL Interface
g02jhf (lmm_fit)