NAG Library Routine Document
G02JCF
1 Purpose
G02JCF preprocesses a dataset prior to fitting a linear mixed effects regression model of the following form via either
G02JDF or
G02JEF.
2 Specification
SUBROUTINE G02JCF ( |
WEIGHT, N, NCOL, DAT, LDDAT, LEVELS, Y, WT, FIXED, LFIXED, NRNDM, RNDM, LDRNDM, NFF, NLSV, NRF, RCOMM, LRCOMM, ICOMM, LICOMM, IFAIL) |
INTEGER |
N, NCOL, LDDAT, LEVELS(NCOL), FIXED(LFIXED), LFIXED, NRNDM, RNDM(LDRNDM,NRNDM), LDRNDM, NFF, NLSV, NRF, LRCOMM, ICOMM(LICOMM), LICOMM, IFAIL |
REAL (KIND=nag_wp) |
DAT(LDDAT,NCOL), Y(N), WT(*), RCOMM(LRCOMM) |
CHARACTER(1) |
WEIGHT |
|
3 Description
G02JCF must be called prior to fitting a linear mixed effects regression model with either
G02JDF or
G02JEF.
The model fitting routines
G02JDF and
G02JEF fit a model of the following form:
where |
is a vector of observations on the dependent variable,
|
|
is an by design matrix of fixed independent variables, |
|
is a vector of unknown fixed effects, |
|
is an by design matrix of random independent variables, |
|
is a vector of length of unknown random effects, |
|
is a vector of length of unknown random errors,
|
and
and
are Normally distributed with expectation zero and variance/covariance matrix defined by
where
,
is the
identity matrix and
is a diagonal matrix.
Case weights can be incorporated into the model by replacing and with and respectively where is a diagonal weight matrix.
4 References
None.
5 Parameters
- 1: WEIGHT – CHARACTER(1)Input
On entry: indicates if weights are to be used.
- No weights are used.
- Case weights are used and must be supplied in array WT.
Constraint:
or .
- 2: N – INTEGERInput
On entry:
, the number of observations.
The effective number of observations, that is the number of observations with nonzero weight (see
WT for more detail), must be greater than the number of fixed effects in the model (as returned in
NFF).
Constraint:
.
- 3: NCOL – INTEGERInput
On entry: the number of columns in the data matrix,
DAT.
Constraint:
.
- 4: DAT(LDDAT,NCOL) – REAL (KIND=nag_wp) arrayInput
On entry: a matrix of data, with
holding the
th observation on the
th variable. The two design matrices
and
are constructed from
DAT and the information given in
FIXED (for
) and
RNDM (for
).
Constraint:
if .
- 5: LDDAT – INTEGERInput
On entry: the first dimension of the array
DAT as declared in the (sub)program from which G02JCF is called.
Constraint:
.
- 6: LEVELS(NCOL) – INTEGER arrayInput
On entry:
contains the number of levels associated with the
th variable held in
DAT.
If the th variable is continuous or binary (i.e., only takes the values zero or one) then must be set to . Otherwise the th variable is assumed to take an integer value between and , (i.e., the th variable is discrete with levels).
Constraint:
, for .
- 7: Y(N) – REAL (KIND=nag_wp) arrayInput
On entry: , the vector of observations on the dependent variable.
- 8: WT() – REAL (KIND=nag_wp) arrayInput
-
Note: the dimension of the array
WT
must be at least
if
.
On entry: if
,
WT must contain the diagonal elements of the weight matrix
.
If , the th observation is not included in the model and the effective number of observations is the number of observations with nonzero weights.
If
,
WT is not referenced and the effective number of observations is
.
Constraint:
if , , for .
- 9: FIXED(LFIXED) – INTEGER arrayInput
On entry: defines the structure of the fixed effects design matrix,
.
- The number of variables, , to include as fixed effects (not including the intercept if present).
- The fixed intercept flag which must contain if a fixed intercept is to be included and otherwise.
- The column of DAT holding the
th fixed variable, for .
See
Section 8.1 for more details on the construction of
.
Constraints:
- ;
- ;
- , for .
- 10: LFIXED – INTEGERInput
On entry: length of the vector
FIXED.
Constraint:
.
- 11: NRNDM – INTEGERInput
On entry: the second dimension of the array
RNDM as declared in the (sub)program from which G02JCF is called.
Constraint:
.
- 12: RNDM(LDRNDM,NRNDM) – INTEGER arrayInput
On entry:
defines the structure of the
random effects design matrix,
. The
th column of
RNDM defines a block of columns in the design matrix
.
- The number of variables, , to include as random effects in the th block (not including the random intercept if present).
- The random intercept flag which must contain if block includes a random intercept and otherwise.
- The column of DAT holding the
th random variable in the th block, for .
- The number of subject variables, , for the th block. The subject variables define the nesting structure for this block.
- The column of DAT holding the
th subject variable in the th block, for .
See
Section 8.2 for more details on the construction of
.
Constraints:
- ;
- ;
- at least one random variable or random intercept must be specified in each block, i.e., ;
- the column identifiers associated with the random variables must be in the range to NCOL, i.e., , for ;
- ;
- the column identifiers associated with the subject variables must be in the range to NCOL, i.e., , for .
- 13: LDRNDM – INTEGERInput
On entry: the first dimension of the array
RNDM as declared in the (sub)program from which G02JCF is called.
Constraint:
.
- 14: NFF – INTEGEROutput
On exit: , the number of fixed effects estimated, i.e., the number of columns in the design matrix .
- 15: NLSV – INTEGEROutput
On exit: the number of levels for the overall subject variable (see
Section 8.2 for a description of what this means). If there is no overall subject variable,
.
- 16: NRF – INTEGEROutput
On exit: the number of random effects estimated in each of the overall subject blocks. The number of columns in the design matrix is given by .
- 17: RCOMM(LRCOMM) – REAL (KIND=nag_wp) arrayCommunication Array
On exit: communication array as required by the analysis routines
G02JDF and
G02JEF.
- 18: LRCOMM – INTEGERInput
On entry: the dimension of the array
RCOMM as declared in the (sub)program from which G02JCF is called.
Constraint:
.
- 19: ICOMM(LICOMM) – INTEGER arrayCommunication Array
On exit: if
,
holds the minimum required value for
LICOMM and
holds the minimum required value for
LRCOMM, otherwise
ICOMM is a communication array as required by the analysis routines
G02JDF and
G02JEF.
- 20: LICOMM – INTEGERInput
On entry: the dimension of the array
ICOMM as declared in the (sub)program from which G02JCF is called.
Constraint:
or
.
where
- ,
- ,
- ,
- ,
- , and
-
- 21: IFAIL – INTEGERInput/Output
On entry:
IFAIL must be set to
,
. 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
is recommended. If the output of error messages is undesirable, then the value
is recommended. Otherwise, if you are not familiar with this parameter, the recommended value is
.
When the value is used it is essential to test the value of IFAIL on exit.
On exit:
unless the routine detects an error or a warning has been flagged (see
Section 6).
6 Error Indicators and Warnings
If on entry
or
, explanatory error messages are output on the current error message unit (as defined by
X04AAF).
Errors or warnings detected by the routine:
On entry, or .
On entry, .
On entry, .
On entry, at least one data point for variable has data taking a value less than one or greater than .
On entry, .
On entry, is less than one for at least one .
On entry, there is at least one negative weight.
On entry, number of fixed parameters is less than zero.
On entry,
LFIXED is too small.
On entry, .
On entry, number of random parameters for at least one random statement is less than zero.
On entry,
LDRNDM is too small.
On entry,
LRCOMM is too small.
On entry,
LICOMM is too small.
On entry, more fixed factors than observations.
On entry, no observations due to zero weights.
On entry, invalid value for fixed intercept flag.
On entry, invalid value for random intercept flag in at least one random statement.
On entry, the identifier associated with at least one fixed variable, as supplied in
FIXED is less than one or greater than
NCOL.
On entry, must be at least one parameter, or an intercept in each random statement.
On entry, the identifier associated with random variable
in random statement
, is less than one or greater than
NCOL.
On entry, number of subject parameters for random statement .
On entry, nesting variable in random statement has less than two levels.
7 Accuracy
Not applicable.
8.1 Construction of the fixed effects design matrix,
Let
- denote the number of fixed variables, that is ;
- denote the th fixed variable, that is the vector of values held in the th column of DAT when ;
- denote the th element of ;
- denote the number of levels for , that is ;
- denoted an indicator function that returns a vector of values whose th element is if and otherwise.
The design matrix for the
fixed effects,
, is constructed as follows:
- set to
one
and the flag to false;
- if a fixed intercept is included, that is ,
- set the first column of to a vector of s;
- set ;
- set to true;
- loop over each fixed variable, so for each ,
- if ,
- set the th column of to be ;
- set ;
- else
- if is false then
- set the columns, to , of to , for ;
- set ;
- set to true;
-
else
- set the columns, to , of to , for ;
- set .
The number of columns in the design matrix,
, is therefore given by
This quantity is returned in
NFF.
In summary, G02JCF converts all non-binary categorical variables (i.e., where
) to dummy variables. If a fixed intercept is included in the model then the first level of all such variables is dropped. If a fixed intercept is not included in the model then the first level of all such variables, other than the first, is dropped. The variables are added into the model in the order they are specified in
FIXED.
8.2 Construction of random effects design matrix,
Let
- denote the number of random variables in the th random statement, that is ;
- denote the th random variable from the th random statement, that is the vector of values held in the th column of DAT when ;
- denote the th element of ;
- denote the number of levels for , that is ;
- denoted an indicator function that returns a vector of values whose th element is if and otherwise;
- denote the number of subject variables in the th random statement, that is ;
- denote the th subject variable from the th random statement, that is the vector of values held in the th column of DAT when ;
- denote the th element of ;
- denote the number of levels for , that is ;
- denoted an indicator function that returns a vector of values whose th element is if for all and otherwise.
The design matrix for the
random effects,
, is constructed as follows:
- set to one;
- loop over each random statement, so for each ,
- loop over each level of the last subject variable, so for each ,
-
- loop over each level of the second subject variable, so for each ,
- loop over each level of the first subject variable, so for each ,
- if a random intercept is included, that is ,
- set the th column of to ;
- set ;
- loop over each random variable in the th random statement, so for each ,
- if ,
- set the th column of to where indicates an element-wise multiplication between the two vectors, and ;
- set ;
- else
- set the columns, to , of to , for . As before, indicates an element-wise multiplication between the two vectors, and ;
- set .
In summary, each column of
RNDM defines a block of consecutive columns in
. G02JCF converts all non-binary categorical variables (i.e., where
or
) to dummy variables. All random variables defined within a column of
RNDM are nested within all subject variables defined in the same column of
RNDM. In addition each of the subject variables are nested within each other, starting with the first (i.e., each of the
are nested within
which in turn is nested within
, which in turn is nested within
, etc.).
If the last subject variable in each column of
RNDM are the same (i.e.,
) then all random effects in the model are nested within this variable. In such instances the last subject variable (
) is called the overall subject variable. The fact that all of the random effects in the model are nested within the overall subject variable means that
is block diagonal in structure. This fact can be utilised to improve the efficiency of the underlying computation and reduce the amount of internal storage required. The number of levels in the overall subject variable is returned in
.
If the last
subject variables in each column of
RNDM are the same, for
then the overall subject variable is defined as the interaction of these
variables and
If there is no overall subject variable then .
The number of columns in the design matrix is given by .
8.3 The RNDM parameter
To illustrate some additional points about the
RNDM parameter, we assume that we have a dataset with three discrete variables,
,
and
, with
and
levels respectively, and that
is in the first column of
DAT,
in the second and
the third. Also assume that we wish to fit a model containing
along with
nested within
, as random effects. In order to do this the
RNDM matrix requires two columns:
The first column,
, indicates one random variable (
), no intercept (
), the random variable is in the first column of
DAT (
), there are no subject variables; as no nesting is required for
(
). The last element in this column is ignored.
The second column,
, indicates one random variable (
), no intercept (
), the random variable is in the second column of
DAT , there is one subject variable (
), and the subject variable is in the third column of
DAT .
The corresponding
matrix would have
columns, with
coming from
and
(
) from
nested within
. The, symmetric,
matrix has the form
where
indicates a structural zero, i.e., it always takes the value
, irrespective of the data, and
a value that is not a structural zero. The first two rows and columns of
correspond to
. The block diagonal matrix in the 12 rows and columns in the bottom right correspond to
nested within
. With the
blocks corresponding to the levels of
. There are three blocks as the subject variable (
) has three levels.
The model fitting routines,
G02JDF and
G02JEF, use the sweep algorithm to calculate the log likelihood function for a given set of variance components. This algorithm consists of moving down the diagonal elements (called pivots) of a matrix which is similar in structure to
, and updating each element in that matrix. When using the
diagonal element of a matrix
, an element
, is adjusted by an amount equal to
. This process can be referred to as sweeping on the
th pivot. As there are no structural zeros in the first row or column of the above
, sweeping on the first pivot of
would alter each element of the matrix and therefore destroy the structural zeros, i.e., we could no longer guarantee they would be zero.
Reordering the
RNDM matrix to
i.e., the swapping the two columns, results in a
matrix of the form
This matrix is identical to the previous one, except the first two rows and columns have become the last two rows and columns. Sweeping a matrix, , of this form on the first pivot will only affect those elements , where , which is only the th and th row and columns, and the top left hand block of rows and columns. The block diagonal nature of the first rows and columns therefore greatly reduces the amount of work the algorithm needs to perform.
G02JCF constructs the
as specified by the
RNDM matrix, and does not attempt to reorder it to improve performance. Therefore for best performance some thought is required on what ordering to use. In general it is more efficient to structure
RNDM in such a way that the first row relates to the deepest level of nesting, the second to the next level, etc..
9 Example
See Section 9 in
G02JDF and
G02JEF.