NAG CL Interface
g02jac (mixeff_reml)
1
Purpose
g02jac fits a linear mixed effects regression model using restricted maximum likelihood (REML).
2
Specification
void |
g02jac (Integer n,
Integer ncol,
const double dat[],
Integer tddat,
const Integer levels[],
Integer yvid,
Integer cwid,
Integer nfv,
const Integer fvid[],
Integer fint,
Integer nrv,
const Integer rvid[],
Integer nvpr,
const Integer vpr[],
Integer rint,
Integer svid,
double gamma[],
Integer *nff,
Integer *nrf,
Integer *df,
double *reml,
Integer lb,
double b[],
double se[],
Integer maxit,
double tol,
Integer *warn,
NagError *fail) |
|
The function may be called by the names: g02jac, nag_correg_mixeff_reml or nag_reml_mixed_regsn.
3
Description
g02jac 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
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 expectations 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
,
g02jac uses an iterative process to estimate
. Due to the iterative nature of the estimation a set of initial values,
, for
is required.
g02jac 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).
g02jac fits the model using a quasi-Newton algorithm to maximize the restricted log-likelihood function:
where
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
-
1:
– Integer
Input
-
On entry: , the number of observations.
Constraint:
.
-
2:
– Integer
Input
-
On entry: the number of columns in the data matrix,
DAT.
Constraint:
.
-
3:
– const double
Input
-
Note: where appears in this document, it refers to the array element
.
On entry: array containing all of the data. For the
th observation:
- holds the dependent variable, ;
- if , holds the case weights;
- if , holds the subject variable.
The remaining columns hold the values of the independent variables.
Constraints:
- if , ;
- if , .
-
4:
– Integer
Input
-
On entry: the stride separating matrix column elements in the array
dat.
Constraint:
.
-
5:
– const Integer
Input
-
On entry:
contains the number of levels associated with the
th variable of the data matrix
DAT. If this variable is continuous or binary (i.e., only takes the values zero or one) then
should be
; if the variable is discrete then
is the number of levels associated with it and
is assumed to take the values
to
, for
.
Constraint:
, for .
-
6:
– Integer
Input
-
On entry: the column of
DAT holding the dependent,
, variable.
Constraint:
.
-
7:
– Integer
Input
-
On entry: the column of
DAT holding the case weights.
If , no weights are used.
Constraint:
.
-
8:
– Integer
Input
-
On entry: the number of independent variables in the model which are to be treated as being fixed.
Constraint:
.
-
9:
– const Integer
Input
-
On entry: the columns of the data matrix
DAT holding the fixed independent variables with
holding the column number corresponding to the
th fixed variable.
Constraint:
, for .
-
10:
– Integer
Input
-
On entry: flag indicating whether a fixed intercept is included ().
Constraint:
or .
-
11:
– Integer
Input
-
On entry: the number of independent variables in the model which are to be treated as being random.
-
12:
– const Integer
Input
-
On entry: the columns of the data matrix holding the random independent variables with holding the column number corresponding to the th random variable.
Constraint:
, for .
-
13:
– Integer
Input
-
On entry: if
and
,
nvpr is the number of variance components being
, (
), else
.
If , is not referenced.
Constraint:
if , .
-
14:
– const Integer
Input
-
On entry: holds a flag indicating the variance of the th random variable. The variance of the th random variable is , where if and and otherwise. Random variables with the same value of are assumed to be taken from the same distribution.
Constraint:
, for .
-
15:
– Integer
Input
-
On entry: flag indicating whether a random intercept is included (
).
If
,
rint is not referenced.
Constraint:
or .
-
16:
– Integer
Input
-
On entry: the column of
DAT holding the subject variable.
If , no subject variable is used.
Specifying a subject variable is equivalent to specifying the interaction between that variable and all of the random-effects. Letting the notation denote the interaction between variables and , fitting a model with , random-effects and subject variable is equivalent to fitting a model with random-effects and no subject variable. If the model is equivalent to fitting and no subject variable.
Constraint:
.
-
17:
– double
Input/Output
-
On entry: holds the initial values of the variance components,
, with
the initial value for
, for
. If
and
,
, else
.
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 .
Constraint:
, for .
-
18:
– Integer *
Output
-
On exit: the number of fixed effects estimated (i.e., the number of columns, , in the design matrix ).
-
19:
– Integer *
Output
-
On exit: the number of random effects estimated (i.e., the number of columns, , in the design matrix ).
-
20:
– Integer *
Output
-
On exit: the degrees of freedom.
-
21:
– double *
Output
-
On exit:
where
is the log of the restricted maximum likelihood calculated at
, the estimated variance components returned in
gamma.
-
22:
– Integer
Input
-
On entry: the size of the array
b.
Constraint:
where if and otherwise.
-
23:
– double
Output
-
On exit: the parameter estimates,
, with the first
nff elements of
b containing the fixed effect parameter estimates,
and the next
nrf elements of
b containing the random effect parameter estimates,
.
Fixed effects
If
,
contains the estimate of the fixed intercept. Let
denote the number of levels associated with the
th fixed variable, that is
. Define
- if , else if , ;
- , .
Then for
:
- if ,
contains the parameter estimate for the th level of the th fixed variable, for ;
- if , contains the parameter estimate for the th fixed variable.
Random effects
Redefining
to denote the number of levels associated with the
th random variable, that is
. Define
- if , else if , ;
, .
Then for
:
- if ,
- if ,
contains the parameter estimate for the th level of the th random variable, for ;
- if , contains the parameter estimate for the th random variable;
- if ,
- let denote the number of levels associated with the subject variable, that is ;
- if ,
contains the parameter estimate for the interaction between the th level of the subject variable and the th level of the th random variable, for and ;
- if ,
contains the parameter estimate for the interaction between the th level of the subject variable and the th random variable, for ;
- if , contains the estimate of the random intercept.
-
24:
– double
Output
-
On exit: the standard errors of the parameter estimates given in
b.
-
25:
– Integer
Input
-
On entry: the maximum number of iterations.
If , the default value of is used.
If
, the parameter estimates
and corresponding standard errors are calculated based on the value of
supplied in
gamma.
-
26:
– double
Input
-
On entry: the tolerance used to assess convergence.
If , the default value of is used, where is the machine precision.
-
27:
– Integer *
Output
-
On exit: is set to
if a variance component was estimated to be a negative value during the fitting process. Otherwise
warn is set to
.
If , the negative estimate is set to zero and the estimation process allowed to continue.
-
28:
– 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_BAD_PARAM
-
On entry, argument had an illegal value.
On entry, invalid data: categorical variable with value greater than that specified in
levels.
- NE_CONV
-
Routine failed to converge in
maxit iterations:
.
See
Section 10 for advice.
- NE_FAIL_TOL
-
Routine failed to converge to specified tolerance:
.
See
Section 10 for advice.
- NE_INT
-
On entry, .
Constraint: or .
On entry,
lb too small:
.
On entry, , for at least one .
On entry, .
Constraint: number of observations with nonzero weights must be greater than one.
On entry, .
Constraint: .
On entry, .
Constraint: , for all .
On entry, .
Constraint: , for all .
On entry, .
Constraint: .
On entry, .
Constraint: , for all .
On entry, .
Constraint: or .
- NE_INT_2
-
On entry, and .
Constraint: and any supplied weights must be .
On entry, and .
Constraint: .
On entry, and .
Constraint: and .
On entry, and .
Constraint: and ( or ).
On entry, and .
Constraint: .
On entry, and .
Constraint: .
On entry, and .
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_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
-
On entry, , for at least one .
- NE_ZERO_DOF_ERROR
-
Degrees of freedom : .
This is due to the number of parameters exceeding the effective number of observations.
7
Accuracy
The accuracy of the results can be adjusted through the use of the
tol argument.
8
Parallelism and Performance
g02jac is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
g02jac 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.
Wherever possible any block structure present in the design matrix
should be modelled through a subject variable, specified via
svid, rather than being explicitly entered into
dat.
g02jac uses an iterative process to fit the specified model and for some problems this process may fail to converge (see
NE_CONV or
NE_FAIL_TOL). If the function fails to converge then the maximum number of iterations (see
maxit) or tolerance (see
tol) may require increasing; try a different starting estimate in
gamma. Alternatively, the model can be fit using maximum likelihood (see
g02jbc) or using the noniterative MIVQUE0.
To fit the model just using MIVQUE0, the first element of
gamma should be set to
and
maxit should be set to zero.
Although the quasi-Newton algorithm used in
g02jac tends to require more iterations before converging compared to the Newton–Raphson algorithm recommended by
Wolfinger et al. (1994), it does not require the second derivatives of the likelihood function to be calculated and consequentially takes significantly less time per iteration.
10
Example
The following dataset is taken from
Stroup (1989) and arises from a balanced split-plot design with the whole plots arranged in a randomized complete block-design.
In this example the full design matrix for the random independent variable,
, is given by:
where
The block structure evident in
(1) is modelled by specifying a four-level subject variable, taking the values
. The first column of
is added to
by setting
. The remaining columns of
are specified by a three level factor, taking the values,
.
10.1
Program Text
10.2
Program Data
10.3
Program Results