The routine may be called by the names g02jbf or nagf_correg_mixeff_ml.
g02jbf fits a model of the form:
is a vector of observations on the dependent variable,
is a known design matrix for the fixed independent variables,
is a vector of length of unknown fixed effects,
is a known design matrix for the random independent variables,
is a vector of length of unknown random effects;
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 , g02jbf uses an iterative process to estimate . Due to the iterative nature of the estimation a set of initial values, , for is required. g02jbf 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).
g02jbf fits the model using a quasi-Newton algorithm to maximize the log-likelihood function:
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 .
Goodnight J H (1979) A tutorial on the SWEEP operator The American Statistician33(3) 149–158
Harville D A (1977) Maximum likelihood approaches to variance component estimation and to related problems JASA72 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 DisciplinesSouthern 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
1: – IntegerInput
On entry: , the number of observations.
2: – IntegerInput
On entry: the number of columns in the data matrix, dat.
3: – IntegerInput
On entry: the first dimension of the array dat as declared in the (sub)program from which g02jbf is called.
4: – Real (Kind=nag_wp) arrayInput
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.
if , ;
if , .
5: – Integer arrayInput
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 .
, for .
6: – IntegerInput
On entry: the column of dat holding the dependent, , variable.
7: – IntegerInput
On entry: the column of dat holding the case weights.
If , no weights are used.
8: – IntegerInput
On entry: the number of independent variables in the model which are to be treated as being fixed.
9: – Integer arrayInput
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.
, for .
10: – IntegerInput
On entry: flag indicating whether a fixed intercept is included ().
11: – IntegerInput
On entry: the number of independent variables in the model which are to be treated as being random.
12: – Integer arrayInput
On entry: the columns of the data matrix holding the random independent variables with holding the column number corresponding to the th random variable.
, for .
13: – IntegerInput
On entry: if and , nvpr is the number of variance components being , (), else .
If , is not referenced.
if , .
14: – Integer arrayInput
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.
, for .
15: – IntegerInput
On entry: flag indicating whether a random intercept is included ().
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.
17: – Real (Kind=nag_wp) arrayInput/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 .
or , for .
18: – IntegerOutput
On exit: the number of fixed effects estimated (i.e., the number of columns, , in the design matrix ).
19: – IntegerOutput
On exit: the number of random effects estimated (i.e., the number of columns, , in the design matrix ).
20: – IntegerOutput
On exit: the degrees of freedom.
21: – Real (Kind=nag_wp)Output
On exit: where is the log of the maximum likelihood calculated at , the estimated variance components returned in gamma.
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, .
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 :
contains the parameter estimate for the th level of the th fixed variable, for ;
if , contains the parameter estimate for the th fixed variable.
Redefining to denote the number of levels associated with the th random variable, that is . Define
if , else if , ; , .
Then for :
contains the parameter estimate for the th level of the th random variable, for ;
if , contains the parameter estimate for the th random variable;
let denote the number of levels associated with the subject variable, that is ;
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 ;
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: – Real (Kind=nag_wp) arrayOutput
On exit: the standard errors of the parameter estimates given in b.
25: – IntegerInput
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: – Real (Kind=nag_wp)Input
On entry: the tolerance used to assess convergence.
If , the default value of is used, where is the machine precision.
27: – IntegerOutput
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: – IntegerInput/Output
On entry: ifail must be set to , or to set behaviour on detection of an error; these values have no effect when no error is detected.
A value of causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of means that an error message is printed while a value of means that it is not.
If halting is not appropriate, the value or is recommended. If message printing is undesirable, then the value is recommended. Otherwise, the value is recommended. When the value or 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).
6Error 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, and .
Constraint: and any supplied weights must be .
Routine failed to converge to specified tolerance: .
An unexpected error has been triggered by this routine. Please
See Section 7 in the Introduction to the NAG Library FL Interface for further information.
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.
Dynamic memory allocation failed.
See Section 9 in the Introduction to the NAG Library FL Interface for further information.
The accuracy of the results can be adjusted through the use of the tol argument.
8Parallelism and Performance
g02jbf is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
g02jbf 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.
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.
g02jbf uses an iterative process to fit the specified model and for some problems this process may fail to converge (see ). If the routine 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 restricted maximum likelihood (see g02jaf) 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 g02jbf 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.
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:
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, .