NAG CL Interface
g02jhc (lmm_fit)
1
Purpose
g02jhc fits a multi-level linear mixed effects regression model using restricted maximum likelihood (REML) or maximum likelihood (ML). Prior to calling
g02jhc the initialization function
g02jfc must be called.
2
Specification
void |
g02jhc (void *hlmm,
Integer nvpr,
double gamma[],
Integer *effn,
Integer *rnkx,
Integer *ncov,
double *lnlike,
Integer lb,
double b[],
double se[],
double czz[],
Integer pdczz,
double cxx[],
Integer pdcxx,
double cxz[],
Integer pdcxz,
double rcomm[],
Integer icomm[],
NagError *fail) |
|
The function may be called by the names: g02jhc or nag_correg_lmm_fit.
3
Description
g02jhc fits a model of the form:
where |
is a vector of observations on the dependent variable, |
|
is a known by design matrix for the fixed independent variables, |
|
is a vector of length of unknown fixed effects, |
|
is a known by design matrix for the random independent variables, |
|
is a vector of length of unknown random effects, |
and |
is a vector of length of unknown random errors. |
Both
and
are assumed to have a Gaussian distribution with expectation zero and variance/covariance matrix defined by
where
,
is the
identity matrix and
is a diagonal matrix. It is assumed that the random variables,
, can be subdivided into
groups with each group being identically distributed with expectation zero and variance
. The diagonal elements of matrix
therefore take one of the values
, 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
variance components
, where
. Rather than working directly with
,
g02jhc uses an iterative process to estimate
. Due to the iterative nature of the estimation a set of initial values,
, for
is required.
g02jhc 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).
g02jhc fits the model by maximizing the restricted log-likelihood function:
or the log-likelihood function:
where
By default the restricted log-likelihood function is used, the log-likelihood function can be chosen through the optional parameter
as detailed in the documentation for
g02jfc.
Once the final estimates for
have been obtained, the value of
is given by
Case weights, , can be incorporated into the model by replacing and with and respectively, for a diagonal weight matrix .
The log-likelihood,
, 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
g02jhc the initialization function
g02jfc must be called, therefore this documentation should be read in conjunction with the document for
g02jfc. 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 functions.
-
1:
– void *
Input
-
On entry: a G22 handle to the internal data structure containing a description of the required model as returned in
hlmm by
g02jfc.
-
2:
– Integer
Input
-
On entry:
, the number of variance components being estimated (excluding the overall variance,
).
This should be set to the value of
nvpr returned by
g02jfc.
-
3:
– double
Input/Output
-
On entry: holds the initial values of the variance components,
, with
the initial value for
, for
.
If
, the remaining elements of
gamma are ignored and the initial values for the variance components are estimated from the data using MIVQUE0.
On exit:
, for
, holds the final estimate of
and
holds the final estimate for
.
Labels for the variance components can be obtained using
g22ydc as demonstrated in
Section 10.
Constraint:
, for .
-
4:
– Integer *
Output
-
On exit: effective number of observations. If there are no weights (i.e., ), or all weights are nonzero, .
-
5:
– Integer *
Output
-
On exit: the rank of the design matrix, , for the fixed effects.
-
6:
– Integer *
Output
-
On exit: number of variance components not estimated to be zero. If none of the variance components are estimated to be zero, .
-
7:
– double *
Output
-
On exit:
where
is the log of the restricted maximum likelihood calculated at
, the estimated variance components returned in
gamma.
-
8:
– Integer
Input
-
On entry: the dimension of the arrays
b and
se.
Constraint:
.
-
9:
– double
Output
-
On exit: the parameter estimates, with the first
elements of
containing the parameter estimates for the random effects,
, and the remaining
elements containing the parameter estimates for the fixed effects,
.
Labels for the parameter estimates can be obtained using
g22ydc as demonstrated in
Section 10.
-
10:
– double
Output
-
On exit: the standard errors of the parameter estimates given in
b.
-
11:
– double
Output
-
Note: where appears in this document, it refers to the array element
.
On exit: if
,
CZZ holds the lower triangular portion of the matrix
, where
and
are the estimates of
and
respectively.
If
, 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.
If
,
czz is not referenced and may be
NULL.
-
12:
– Integer
Input
-
On entry: the stride separating matrix row elements in the array
czz.
Constraints:
- or ;
- if , .
(see
g02jfc)
-
13:
– double
Output
-
Note: where appears in this document, it refers to the array element
.
On exit: if
,
CXX holds the lower triangular portion of the matrix
, where
is the estimated value of
.
If
, 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.
If
,
cxx is not referenced and may be
NULL.
-
14:
– Integer
Input
-
On entry: the stride separating matrix row elements in the array
cxx.
Constraint:
or
.
(see
g02jfc)
-
15:
– double
Output
-
Note: where appears in this document, it refers to the array element
.
On exit:
CXZ holds the matrix
, where
and
are the estimates of
and
respectively.
If
,
cxz is not referenced and may be
NULL.
-
16:
– Integer
Input
-
On entry: the stride separating matrix row elements in the array
cxz.
Constraint:
or .
(see )
-
17:
– double
Communication Array
-
On entry: communication array initialized by a call to
g02jfc.
-
18:
– Integer
Communication Array
-
On entry: communication array initialized by a call to
g02jfc.
-
19:
– NagError *
Input/Output
-
The NAG error argument (see
Section 7 in the Introduction to the NAG Library CL Interface).
6
Error Indicators and Warnings
- NE_ALLOC_FAIL
-
Dynamic memory allocation failed.
See
Section 3.1.2 in the Introduction to the NAG Library CL Interface for further information.
- NE_ARRAY_SIZE
-
On entry, .
Constraint: .
On entry, .
Constraint: .
On entry, .
Constraint: .
On entry, .
Constraint: .
- NE_BAD_PARAM
-
On entry, argument had an illegal value.
- NE_HANDLE
-
hlmm has not been initialized or is corrupt.
hlmm is not a G22 handle as generated by
g02jfc.
- NE_ILLEGAL_COMM
-
On entry, the communication arrays,
icomm and
rcomm, have not been initialized correctly.
- NE_INT
-
On entry, .
Constraint: .
- NE_INTERNAL_ERROR
-
An internal error has occurred in this function. Check the function call and any array sizes. If the call is correct then please contact
NAG for assistance.
See
Section 7.5 in the Introduction to the NAG Library CL Interface for further information.
- NE_NEG_ELEMENT
-
At least one negative estimate for
gamma was obtained. All negative estimates have been set to zero.
- NE_NO_LICENCE
-
Your licence key may have expired or may not have been installed correctly.
See
Section 8 in the Introduction to the NAG Library CL Interface for further information.
- NE_REAL_ARRAY
-
On entry, .
Constraint: or .
- NW_KT_CONDITIONS
-
Current point cannot be improved upon.
- NW_NOT_CONVERGED
-
Optimal solution found, but requested accuracy not achieved.
- NW_TOO_MANY_ITER
-
Too many major iterations.
7
Accuracy
Not applicable.
8
Parallelism and Performance
g02jhc is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
g02jhc 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 function. Please also consult the
Users' Note for your implementation for any additional implementation-specific information.
If
g02jgc has been called, then rather than use the values of
fnlsv,
nff,
rnlsv,
nrf and
nvpr as returned by
g02jfc they should be obtained by calling
g02zlc. See
Section 9 in
g02jgc 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
observations and
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
g22ydc. See
Section 10 in
g02jfc for an example that uses the standard parameter labels directly.
10.1
Program Text
10.2
Program Data
10.3
Program Results