NAG FL Interface
g02jaf (mixeff_reml)

Note: this routine is deprecated. Replaced by g02jff followed by g02jhf.

The new suite of routines consists of:

g22ybf used to describe the dataset of interest. Calling this routine allows labels to be assign to variables, which can then be used when specifying the model.
g22yaf multiple calls to this routine are used to specify the fixed and random part of the model. The model is specified using strings and a modelling language, for example the string: $V_{1} + V_{2} + V_{1} \dots V_{1}$ would specify a model with the main effects of variables $V_{1}$ and $V_{2}$ and the interaction term between them. The modelling language is explained in detail in Section 3 in g22yaf.
g02jff pre-processes the dataset prior to calling the model fitting routine.
g02jhf fits the model and returns the parameter estimates etc.

In addition to the routines listed above, the following can also be used:

g02jgf combines information returned by multiple calls to g02jff. This is useful for large problems as it allows the dataset to be split up into smaller subsets of data, pre-processing each one separately before combining them into a single set of information as if g02jff had been called on the full dataset.
g22ydf can be used to obtain labels for the parameter estimates returned by g02jhf.
g22zmf can be used to set any optional arguments.
g22znf can be used to return the value of any optional arguments.

By default, the model fitting routine, g02jhf, fits the linear mixed effects model using restricted maximum likelihood (REML). In order to fit the model using maximum likelihood (ML) you need to call the optional argument setting routine, g22zmf with optstr set to

Likelihood = ML

, between the call to g02jff and the call to g02jhf.

Keyword Search:

NAG Library Manual, Mark 28.6

Interfaces: FL CL CPP AD PY MB

NAG FL Interface Introduction

G02 (Correg) Chapter Contents

G02 (Correg) Chapter Introduction

g02ja: FL CL CPP AD PY MB

▸▿ 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

1 Purpose

g02jaf fits a linear mixed effects regression model using restricted maximum likelihood (REML).

2 Specification

Fortran Interface

Subroutine g02jaf (

n, ncol, lddat, dat, levels, yvid, cwid, nfv, fvid, fint, nrv, rvid, nvpr, vpr, rint, svid, gamma, nff, nrf, df, reml, lb, b, se, maxit, tol, warn, ifail)

Integer, Intent (In)	::	n, ncol, lddat, levels(ncol), yvid, cwid, nfv, fvid(nfv), fint, nrv, rvid(nrv), nvpr, vpr(nrv), rint, svid, lb, maxit
Integer, Intent (Inout)	::	ifail
Integer, Intent (Out)	::	nff, nrf, df, warn
Real (Kind=nag_wp), Intent (In)	::	dat(lddat,ncol), tol
Real (Kind=nag_wp), Intent (Inout)	::	gamma(nvpr+2)
Real (Kind=nag_wp), Intent (Out)	::	reml, b(lb), se(lb)

C Header Interface

#include <nag.h>

void

g02jaf_ (const Integer *n, const Integer *ncol, const Integer *lddat, const double dat[], const Integer levels[], const Integer *yvid, const Integer *cwid, const Integer *nfv, const Integer fvid[], const Integer *fint, const Integer *nrv, const Integer rvid[], const Integer *nvpr, const Integer vpr[], const Integer *rint, const Integer *svid, double gamma[], Integer *nff, Integer *nrf, Integer *df, double *reml, const Integer *lb, double b[], double se[], const Integer *maxit, const double *tol, Integer *warn, Integer *ifail)

The routine may be called by the names g02jaf or nagf_correg_mixeff_reml.

3 Description

g02jaf 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 \times p$ design matrix for the fixed independent variables,
$β$ is a vector of length $p$ of unknown fixed effects,
$Z$ is a known $n \times 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

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 expectations 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

γ

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

g02jaf fits the model using a quasi-Newton algorithm to maximize the restricted log-likelihood function:

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

where

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

Once the final estimates for

γ^{*}

have been obtained, the value of

σ_{R}^{2}

is given by:

σ_{R}^{2} = (r^{'} V^{- 1} r) / (n - p) .

Case weights,

W_{c}

, can be incorporated into the model by replacing

X^{'} X

and

Z^{'} Z

with

X^{'} W_{c} X

and

Z^{'} 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

1: $n$ – Integer Input

On entry:

n

, the number of observations.

Constraint:

n \geq 1

2: $ncol$ – Integer Input

On entry: the number of columns in the data matrix, dat.

Constraint:

ncol \geq 1

3: $lddat$ – Integer Input

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

Constraint:

lddat \geq n

4: $dat (lddat, ncol)$ – Real (Kind=nag_wp) array Input

On entry: array containing all of the data. For the

i

th observation:

$dat (i, yvid)$ holds the dependent variable, $y$ ;
if $cwid \neq 0$ , $dat (i, cwid)$ holds the case weights;
if $svid \neq 0$ , $dat (i, svid)$ holds the subject variable.

The remaining columns hold the values of the independent variables.

Constraints:

if $cwid \neq 0$ , $dat (i, cwid) \geq 0.0$ ;
if $levels (j) \neq 1$ , $1 \leq dat (i, j) \leq levels (j)$ .

5: $levels (ncol)$ – Integer array Input

On entry:

levels (i)

contains the number of levels associated with the

i

th variable of the data matrix dat. If this variable is continuous or binary (i.e., only takes the values zero or one) then

levels (i)

should be

1

; if the variable is discrete then

levels (i)

is the number of levels associated with it and

dat (j, i)

is assumed to take the values

1

levels (i)

, for

j = 1, 2, \dots, n

Constraint:

levels (i) \geq 1

, for

i = 1, 2, \dots, ncol

6: $yvid$ – Integer Input

On entry: the column of dat holding the dependent,

y

, variable.

Constraint:

1 \leq yvid \leq ncol

7: $cwid$ – Integer Input

On entry: the column of dat holding the case weights.

cwid = 0

, no weights are used.

Constraint:

0 \leq cwid \leq ncol

8: $nfv$ – Integer Input

On entry: the number of independent variables in the model which are to be treated as being fixed.

Constraint:

0 \leq nfv < ncol

9: $fvid (nfv)$ – Integer array Input

On entry: the columns of the data matrix dat holding the fixed independent variables with

fvid (i)

holding the column number corresponding to the

i

th fixed variable.

Constraint:

1 \leq fvid (i) \leq ncol

, for

i = 1, 2, \dots, nfv

10: $fint$ – Integer Input

On entry: flag indicating whether a fixed intercept is included (

fint = 1

Constraint:

fint = 0

1

11: $nrv$ – Integer Input

On entry: the number of independent variables in the model which are to be treated as being random.

Constraints:

$0 \leq nrv < ncol$ ;
$nrv + rint > 0$ .

12: $rvid (nrv)$ – Integer array Input

On entry: the columns of the data matrix

dat

holding the random independent variables with

rvid (i)

holding the column number corresponding to the

i

th random variable.

Constraint:

1 \leq rvid (i) \leq ncol

, for

i = 1, 2, \dots, nrv

13: $nvpr$ – Integer Input

On entry: if

rint = 1

and

svid \neq 0

, nvpr is the number of variance components being

estimated - 2

, (

g - 1

), else

nvpr = g

nrv = 0

nvpr

is not referenced.

Constraint: if

nrv \neq 0

1 \leq nvpr \leq nrv

14: $vpr (nrv)$ – Integer array Input

On entry:

vpr (i)

holds a flag indicating the variance of the

i

th random variable. The variance of the

i

th random variable is

σ_{j}^{2}

, where

j = vpr (i) + 1

rint = 1

and

svid \neq 0

and

j = vpr (i)

otherwise. Random variables with the same value of

j

are assumed to be taken from the same distribution.

Constraint:

1 \leq vpr (i) \leq nvpr

, for

i = 1, 2, \dots, nrv

15: $rint$ – Integer Input

On entry: flag indicating whether a random intercept is included (

rint = 1

svid = 0

, rint is not referenced.

Constraint:

rint = 0

1

16: $svid$ – Integer Input

On entry: the column of dat holding the subject variable.

svid = 0

, 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

Z_{1} \times Z_{S}

denote the interaction between variables

Z_{1}

and

Z_{S}

, fitting a model with

rint = 0

, random-effects

Z_{1} + Z_{2}

and subject variable

Z_{S}

is equivalent to fitting a model with random-effects

Z_{1} \times Z_{S} + Z_{2} \times Z_{S}

and no subject variable. If

rint = 1

the model is equivalent to fitting

Z_{S} + Z_{1} \times Z_{S} + Z_{2} \times Z_{S}

and no subject variable.

Constraint:

0 \leq svid \leq ncol

17: $gamma (nvpr + 2)$ – 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, g

. If

rint = 1

and

svid \neq 0

g = nvpr + 1

, else

g = 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, g

, holds the final estimate of

σ_{i}^{2}

and

gamma (g + 1)

holds the final estimate for

σ_{R}^{2}

Constraint:

gamma (1) = - 1.0

gamma (i) \geq 0.0

, for

i = 1, 2, \dots, g

18: $nff$ – Integer Output

On exit: the number of fixed effects estimated (i.e., the number of columns,

p

, in the design matrix

X

19: $nrf$ – Integer Output

On exit: the number of random effects estimated (i.e., the number of columns,

q

, in the design matrix

Z

20: $df$ – Integer Output

On exit: the degrees of freedom.

21: $reml$ – 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.

22: $lb$ – Integer Input

On entry: the size of the array b.

Constraint:

lb \geq fint + \sum_{i = 1}^{nfv} \max (levels (fvid (i)) - 1, 1) + L_{S} \times (rint + \sum_{i = 1}^{nrv} levels (rvid (i)))

where

L_{S} = levels (svid)

svid \neq 0

and

1

otherwise.

23: $b (lb)$ – Real (Kind=nag_wp) array 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

fint = 1

b (1)

contains the estimate of the fixed intercept. Let

L_{i}

denote the number of levels associated with the

i

th fixed variable, that is

L_{i} = levels (fvid (i))

. Define

if $fint = 1$ , $F_{1} = 2$ else if $fint = 0$ , $F_{1} = 1$ ;
$F_{i + 1} = F_{i} + \max (L_{i} - 1, 1)$ , $i \geq 1$ .

Then for

i = 1, 2, \dots, nfv

if $L_{i} > 1$ , $b (F_{i} + j - 2)$ contains the parameter estimate for the $j$ th level of the $i$ th fixed variable, for $j = 2, 3, \dots, L_{i}$ ;
if $L_{i} \leq 1$ , $b (F_{i})$ contains the parameter estimate for the $i$ th fixed variable.

Random effects

Redefining

L_{i}

to denote the number of levels associated with the

i

th random variable, that is

L_{i} = levels (rvid (i))

. Define

if $rint = 1$ , $R_{1} = 2$ else if $rint = 0$ , $R_{1} = 1$ ;
$R_{i + 1} = R_{i} + L_{i}$ , $i \geq 1$ .

Then for

i = 1, 2, \dots, nrv

if $svid = 0$ ,
- if $L_{i} > 1$ , $b (nff + R_{i} + j - 1)$ contains the parameter estimate for the $j$ th level of the $i$ th random variable, for $j = 1, 2, \dots, L_{i}$ ;
- if $L_{i} \leq 1$ , $b (nff + R_{i})$ contains the parameter estimate for the $i$ th random variable;
if $svid \neq 0$ ,
- let $L_{S}$ denote the number of levels associated with the subject variable, that is $L_{S} = levels (svid)$ ;
- if $L_{i} > 1$ , $b (nff + (s - 1) L_{S} + R_{i} + j - 1)$ contains the parameter estimate for the interaction between the $s$ th level of the subject variable and the $j$ th level of the $i$ th random variable, for $s = 1, 2, \dots, L_{S}$ and $j = 1, 2, \dots, L_{i}$ ;
- if $L_{i} \leq 1$ , $b (nff + (s - 1) L_{S} + R_{i})$ contains the parameter estimate for the interaction between the $s$ th level of the subject variable and the $i$ th random variable, for $s = 1, 2, \dots, L_{S}$ ;
- if $rint = 1$ , $b (nff + 1)$ contains the estimate of the random intercept.

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

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

25: $maxit$ – Integer Input

On entry: the maximum number of iterations.

maxit < 0

, the default value of

100

is used.

maxit = 0

, the parameter estimates

(β, ν)

and corresponding standard errors are calculated based on the value of

γ_{0}

supplied in gamma.

26: $tol$ – Real (Kind=nag_wp) Input

On entry: the tolerance used to assess convergence.

tol \leq 0.0

, the default value of

ε^{0.7}

is used, where

ε

is the machine precision.

27: $warn$ – Integer Output

On exit: is set to

1

if a variance component was estimated to be a negative value during the fitting process. Otherwise warn is set to

0

warn = 1

, the negative estimate is set to zero and the estimation process allowed to continue.

28: $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 = 1$: On entry, $cwid = ⟨ value ⟩$ and $ncol = ⟨ value ⟩$ .
Constraint: $0 \leq cwid \leq ncol$ and any supplied weights must be $\geq 0.0$ .

On entry, $fint = ⟨ value ⟩$ .
Constraint: $fint = 0$ or $1$ .

On entry, lb too small: $lb = ⟨ value ⟩$ .

On entry, $lddat = ⟨ value ⟩$ and $n = ⟨ value ⟩$ .
Constraint: $lddat \geq n$ .

On entry, $n = ⟨ value ⟩$ .
Constraint: $n \geq 1$ .

On entry, $n = ⟨ value ⟩$ .
Constraint: number of observations with nonzero weights must be greater than one.

On entry, $ncol = ⟨ value ⟩$ .
Constraint: $ncol \geq 1$ .

On entry, $nfv = ⟨ value ⟩$ and $ncol = ⟨ value ⟩$ .
Constraint: $0 \leq nfv < ncol$ .

On entry, $nrv = ⟨ value ⟩$ and $ncol = ⟨ value ⟩$ .
Constraint: $0 \leq nrv < ncol$ and $nrv + rint > 0$ .

On entry, $nvpr = ⟨ value ⟩$ and $nrv = ⟨ value ⟩$ .
Constraint: $0 \leq nvpr \leq nrv$ and ( $nrv \neq 0$ or $nvpr \geq 1$ ).

On entry, $rint = ⟨ value ⟩$ .
Constraint: $rint = 0$ or $1$ .

On entry, $svid = ⟨ value ⟩$ and $ncol = ⟨ value ⟩$ .
Constraint: $0 \leq svid \leq ncol$ .

On entry, $yvid = ⟨ value ⟩$ and $ncol = ⟨ value ⟩$ .
Constraint: $1 \leq yvid \leq ncol$ .

$ifail = 2$: On entry, $gamma (i) < 0.0$ , for at least one $i$ .

On entry, invalid data: categorical variable with value greater than that specified in levels.

On entry, $levels (i) < 1$ , for at least one $i$ .

On entry, $ncol = ⟨ value ⟩$ .
Constraint: $1 \leq fvid (i) \leq ncol$ , for all $i$ .

On entry, $ncol = ⟨ value ⟩$ .
Constraint: $1 \leq rvid (i) \leq ncol$ , for all $i$ .

On entry, $nvpr = ⟨ value ⟩$ .
Constraint: $1 \leq vpr (i) \leq nvpr$ , for all $i$ .

$ifail = 3$: Degrees of freedom $< 1$ : $df = ⟨ value ⟩$ .
This is due to the number of parameters exceeding the effective number of observations.

$ifail = 4$: Routine failed to converge in maxit iterations: $maxit = ⟨ value ⟩$ .
See Section 10 for advice.

Routine failed to converge to specified tolerance: $tol = ⟨ value ⟩$ .
See Section 10 for advice.

$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

The accuracy of the results can be adjusted through the use of the tol argument.

8 Parallelism and Performance

Background information to multithreading can be found in the Multithreading documentation.

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

g02jaf 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

Wherever possible any block structure present in the design matrix

Z

should be modelled through a subject variable, specified via svid, rather than being explicitly entered into dat.

g02jaf uses an iterative process to fit the specified model and for some problems this process may fail to converge (see

ifail = 4

). 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 maximum likelihood (see g02jbf) or using the noniterative MIVQUE0.

To fit the model just using MIVQUE0, the first element of gamma should be set to

- 1.0

and maxit should be set to zero.

Although the quasi-Newton algorithm used in g02jaf 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,

Z

, is given by:

Z = (\begin{matrix} 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 \\ 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 \end{matrix})

= (\begin{matrix} A & 0 & 0 & 0 \\ 0 & A & 0 & 0 \\ 0 & 0 & A & 0 \\ 0 & 0 & 0 & A \\ A & 0 & 0 & 0 \\ 0 & A & 0 & 0 \\ 0 & 0 & A & 0 \\ 0 & 0 & 0 & A \end{matrix}),

(1)

where

A = (\begin{matrix} 1 & 1 & 0 & 0 \\ 1 & 0 & 1 & 0 \\ 1 & 0 & 0 & 1 \end{matrix}) .

The block structure evident in (1) is modelled by specifying a four-level subject variable, taking the values

{1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4}

. The first column of

1 s

is added to

A

by setting

rint = 1

. The remaining columns of

A

are specified by a three level factor, taking the values,

{1, 2, 3, 1, 2, 3, 1, \dots}

g02ja: FL CL CPP AD PY MB

NAG FL Interfaceg02jaf (mixeff_​reml)

▸▿ 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
g02jaf (mixeff_reml)