Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox: nag_correg_mixeff_hier_ml (g02je)

## Purpose

nag_correg_mixeff_hier_ml (g02je) fits a multi-level linear mixed effects regression model using maximum likelihood (ML). Prior to calling nag_correg_mixeff_hier_ml (g02je) the initialization function nag_correg_mixeff_hier_init (g02jc) must be called.

## Syntax

[gamma, effn, rnkx, ncov, lnlike, id, b, se, czz, cxx, cxz, ifail] = g02je(vpr, nvpr, gamma, lb, rcomm, icomm, iopt, ropt, 'lvpr', lvpr, 'liopt', liopt, 'lropt', lropt)
[gamma, effn, rnkx, ncov, lnlike, id, b, se, czz, cxx, cxz, ifail] = nag_correg_mixeff_hier_ml(vpr, nvpr, gamma, lb, rcomm, icomm, iopt, ropt, 'lvpr', lvpr, 'liopt', liopt, 'lropt', lropt)

## Description

nag_correg_mixeff_hier_ml (g02je) 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$ by $p$ design matrix for the fixed independent variables, $\beta$ is a vector of length $p$ of unknown fixed effects, $Z$ is a known $n$ by $q$ design matrix for the random independent variables, $\nu$ is a vector of length $q$ of unknown random effects, and $\epsilon$ is a vector of length $n$ of unknown random errors.
Both $\nu$ and $\epsilon$ are assumed to have a Gaussian distribution with expectation zero and variance/covariance matrix defined by
 $Var ν ε = G 0 0 R$
where $R={\sigma }_{R}^{2}I$, $I$ is the $n×n$ identity matrix and $G$ is a diagonal matrix. It is assumed that the random variables, $Z$, can be subdivided into $g\le q$ groups with each group being identically distributed with expectation zero and variance ${\sigma }_{i}^{2}$. The diagonal elements of matrix $G$ therefore take one of the values $\left\{{\sigma }_{i}^{2}:i=1,2,\dots ,g\right\}$, depending on which group the associated random variable belongs to.
The model therefore contains three sets of unknowns: the fixed effects $\beta$, the random effects $\nu$ and a vector of $g+1$ variance components $\gamma$, where $\gamma =\left\{{\sigma }_{1}^{2},{\sigma }_{2}^{2},\dots ,{\sigma }_{g-1}^{2},{\sigma }_{g}^{2},{\sigma }_{R}^{2}\right\}$. Rather than working directly with $\gamma$, nag_correg_mixeff_hier_ml (g02je) uses an iterative process to estimate ${\gamma }^{*}=\left\{{\sigma }_{1}^{2}/{\sigma }_{R}^{2},{\sigma }_{2}^{2}/{\sigma }_{R}^{2},\dots ,{\sigma }_{g-1}^{2}/{\sigma }_{R}^{2},{\sigma }_{g}^{2}/{\sigma }_{R}^{2},1\right\}$. Due to the iterative nature of the estimation a set of initial values, ${\gamma }_{0}$, for ${\gamma }^{*}$ is required. nag_correg_mixeff_hier_ml (g02je) 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).
nag_correg_mixeff_hier_ml (g02je) fits the model by maximizing the log-likelihood function:
 $-2 l R = log V + n log rT V-1 r + log 2 π / n$
where
 $V = ZG ZT + R, r=y-Xb and b = XT V-1 X -1 XT V-1 y .$
Once the final estimates for ${\gamma }^{*}$ have been obtained, the value of ${\sigma }_{R}^{2}$ is given by
 $σR2 = rT V-1 r / n - p .$
Case weights, ${W}_{c}$, can be incorporated into the model by replacing ${X}^{\mathrm{T}}X$ and ${Z}^{\mathrm{T}}Z$ with ${X}^{\mathrm{T}}{W}_{c}X$ and ${Z}^{\mathrm{T}}{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).

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

## Parameters

Note: prior to calling nag_correg_mixeff_hier_ml (g02je) the initialization function nag_correg_mixeff_hier_init (g02jc) must be called, therefore this documention should be read in conjunction with the document for nag_correg_mixeff_hier_init (g02jc).
In particular some argument names and conventions described in that document are also relevant here, but their definition has not been repeated. Specifically, rndm, weight, n, nff, nrf, nlsv, levels, fixed, dat, licomm and lrcomm should be interpreted identically in both functions.

### Compulsory Input Parameters

1:     $\mathrm{vpr}\left({\mathbf{lvpr}}\right)$int64int32nag_int array
A vector of flags indicating the mapping between the random variables specified in rndm and the variance components, ${\sigma }_{i}^{2}$. See Further Comments for more details.
Constraint: $1\le {\mathbf{vpr}}\left(\mathit{i}\right)\le {\mathbf{nvpr}}$, for $\mathit{i}=1,2,\dots ,{\mathbf{lvpr}}$.
2:     $\mathrm{nvpr}$int64int32nag_int scalar
$g$, the number of variance components being estimated (excluding the overall variance, ${\sigma }_{R}^{2}$).
Constraint: $1\le {\mathbf{nvpr}}\le {\mathbf{lvpr}}$.
3:     $\mathrm{gamma}\left({\mathbf{nvpr}}+1\right)$ – double array
Holds the initial values of the variance components, ${\gamma }_{0}$, with ${\mathbf{gamma}}\left(\mathit{i}\right)$ the initial value for ${\sigma }_{\mathit{i}}^{2}/{\sigma }_{R}^{2}$, for $\mathit{i}=1,2,\dots ,{\mathbf{nvpr}}$.
If ${\mathbf{gamma}}\left(1\right)=-1.0$, the remaining elements of gamma are ignored and the initial values for the variance components are estimated from the data using MIVQUE0.
Constraint: ${\mathbf{gamma}}\left(1\right)=-1.0\text{​ or ​}{\mathbf{gamma}}\left(\mathit{i}\right)\ge 0.0$, for $\mathit{i}=1,2,\dots ,g$.
4:     $\mathrm{lb}$int64int32nag_int scalar
The dimension of the arrays b and se and the second dimension of the array id.
Constraint: ${\mathbf{lb}}\ge {\mathbf{nff}}+{\mathbf{nrf}}×{\mathbf{nlsv}}$.
5:     $\mathrm{rcomm}\left(:\right)$ – double array
The dimension of the array rcomm must be at least ${\mathbf{lrcomm}}$ (see nag_correg_mixeff_hier_init (g02jc))
Communication array initialized by a call to nag_correg_mixeff_hier_init (g02jc).
6:     $\mathrm{icomm}\left(:\right)$int64int32nag_int array
The dimension of the array icomm must be at least ${\mathbf{licomm}}$ (see nag_correg_mixeff_hier_init (g02jc))
Communication array initialized by a call to nag_correg_mixeff_hier_init (g02jc).
7:     $\mathrm{iopt}\left({\mathbf{liopt}}\right)$int64int32nag_int array
Optional parameters passed to the optimization function.
By default nag_correg_mixeff_hier_ml (g02je) fits the specified model using a modified Newton optimization algorithm as implemented in nag_opt_bounds_mod_deriv2_comp (e04lb). In some cases, where the calculation of the derivatives is computationally expensive it may be more efficient to use a sequential QP algorithm. The sequential QP algorithm as implemented in nag_opt_nlp1_solve (e04uc) can be chosen by setting ${\mathbf{iopt}}\left(5\right)=1$. If ${\mathbf{liopt}}<5$ or ${\mathbf{iopt}}\left(5\right)\ne 1$ then nag_opt_bounds_mod_deriv2_comp (e04lb) will be used.
Different optional parameters are available depending on the optimization function used. In all cases, using a value of $-1$ will cause the default value to be used. In addition only the first liopt values of iopt are used, so for example, if only the first element of iopt needs changing and default values for all other optional parameters are sufficient liopt can be set to $1$.
nag_opt_bounds_mod_deriv2_comp (e04lb) is being used.
 $i$ Description Equivalent nag_opt_bounds_mod_deriv2_comp (e04lb) argument Default Value $1$ Number of iterations maxcal $1000$ $2$ Unit number for monitoring information n/a As returned by nag_file_set_unit_advisory (x04ab) $3$ Print optional parameters ($1=$ print) n/a $-1$ (no printing performed) $4$ Frequency that monitoring information is printed iprint $-1$ $5$ Optimizer used n/a n/a
If requested, monitoring information is displayed in a similar format to that given by nag_opt_bounds_mod_deriv2_comp (e04lb).
nag_opt_nlp1_solve (e04uc) is being used.
 $i$ Description Equivalent nag_opt_nlp1_solve (e04uc) argument Default Value $1$ Number of iterations Major Iteration Limit $\mathrm{max}\left(50,3×{\mathbf{nvpr}}\right)$ $2$ Unit number for monitoring information n/a As returned by nag_file_set_unit_advisory (x04ab) $3$ Print optional parameters ($1=\text{}$ print, otherwise no print) List/Nolist $-1$ (no printing performed) $4$ Frequency that monitoring information is printed Major Print Level $0$ $5$ Optimizer used n/a n/a $6$ Number of minor iterations Minor Iteration Limit $\mathrm{max}\left(50,3×{\mathbf{nvpr}}\right)$ $7$ Frequency that additional monitoring information is printed Minor Print Level $0$
If ${\mathbf{liopt}}\le 0$ then default values are used for all optional parameters and iopt is not referenced.
8:     $\mathrm{ropt}\left({\mathbf{lropt}}\right)$ – double array
Optional parameters passed to the optimization function.
Different optional parameters are available depending on the optimization function used. In all cases, using a value of $-1.0$ will cause the default value to be used. In addition only the first lropt values of ropt are used, so for example, if only the first element of ropt needs changing and default values for all other optional parameters are sufficient lropt can be set to $1$.
nag_opt_bounds_mod_deriv2_comp (e04lb) is being used.
 $i$ Description Equivalent nag_opt_bounds_mod_deriv2_comp (e04lb) argument Default Value $1$ Sweep tolerance n/a $\mathrm{max}\left(\sqrt{\mathrm{eps}},\sqrt{\mathrm{eps}}×\underset{i}{\mathrm{max}}\phantom{\rule{0.25em}{0ex}}\left({\mathrm{zz}}_{ii}\right)\right)$ $2$ Lower bound for ${\gamma }^{*}$ n/a $\mathrm{eps}/100$ $3$ Upper bound for ${\gamma }^{*}$ n/a ${10}^{20}$ $4$ Accuracy of linear minimizations eta $0.9$ $5$ Accuracy to which solution is required xtol $0.0$ $6$ Initial distance from solution stepmx $100000.0$
nag_opt_nlp1_solve (e04uc) is being used.
 $i$ Description Equivalent nag_opt_nlp1_solve (e04uc) argument Default Value $1$ Sweep tolerance n/a $\mathrm{max}\left(\sqrt{\mathrm{eps}},\sqrt{\mathrm{eps}}×\underset{i}{\mathrm{max}}\phantom{\rule{0.25em}{0ex}}\left({\mathrm{zz}}_{ii}\right)\right)$ $2$ Lower bound for ${\gamma }^{*}$ n/a $\mathrm{eps}/100$ $3$ Upper bound for ${\gamma }^{*}$ n/a ${10}^{20}$ $4$ Line search tolerance Line Search Tolerance $0.9$ $5$ Optimality tolerance Optimality Tolerance ${\mathrm{eps}}^{0.72}$
where $\mathrm{eps}$ is the machine precision returned by nag_machine_precision (x02aj) and ${\mathrm{zz}}_{ii}$ denotes the $i$ diagonal element of ${Z}^{\mathrm{T}}Z$.
If ${\mathbf{lropt}}\le 0$ then default values are used for all optional parameters and ropt is not referenced.

### Optional Input Parameters

1:     $\mathrm{lvpr}$int64int32nag_int scalar
Default: the dimension of the array vpr.
The sum of the number of random parameters and the random intercept flags specified in the call to nag_correg_mixeff_hier_init (g02jc).
Constraint: ${\mathbf{lvpr}}={\sum }_{i}{\mathbf{rndm}}\left(1,i\right)+{\mathbf{rndm}}\left(2,i\right)$.
2:     $\mathrm{liopt}$int64int32nag_int scalar
Default: the dimension of the array iopt.
Length of the options array iopt.
3:     $\mathrm{lropt}$int64int32nag_int scalar
Default: the dimension of the array ropt.
Length of the options array ropt.

### Output Parameters

1:     $\mathrm{gamma}\left({\mathbf{nvpr}}+1\right)$ – double array
${\mathbf{gamma}}\left(\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,{\mathbf{nvpr}}$, holds the final estimate of ${\sigma }_{\mathit{i}}^{2}$ and ${\mathbf{gamma}}\left({\mathbf{nvpr}}+1\right)$ holds the final estimate for ${\sigma }_{R}^{2}$.
2:     $\mathrm{effn}$int64int32nag_int scalar
Effective number of observations. If there are no weights (i.e., ${\mathbf{weight}}=\text{'U'}$), or all weights are nonzero, then ${\mathbf{effn}}={\mathbf{n}}$.
3:     $\mathrm{rnkx}$int64int32nag_int scalar
The rank of the design matrix, $X$, for the fixed effects.
4:     $\mathrm{ncov}$int64int32nag_int scalar
Number of variance components not estimated to be zero. If none of the variance components are estimated to be zero, then ${\mathbf{ncov}}={\mathbf{nvpr}}$.
5:     $\mathrm{lnlike}$ – double scalar
$-2{l}_{R}\left(\stackrel{^}{\gamma }\right)$ where ${l}_{R}$ is the log of the maximum likelihood calculated at $\stackrel{^}{\gamma }$, the estimated variance components returned in gamma.
6:     $\mathrm{id}\left(\mathit{ldid},{\mathbf{lb}}\right)$int64int32nag_int array
An array describing the parameter estimates returned in b. The first ${\mathbf{nlsv}}×{\mathbf{nrf}}$ columns of id describe the parameter estimates for the random effects and the last nff columns the parameter estimates for the fixed effects.
The example program for this function includes a demonstration of decoding the parameter estimates given in b using information from id.
For fixed effects:
• for $l={\mathbf{nrf}}×{\mathbf{nlsv}}+1,\dots ,{\mathbf{nrf}}×{\mathbf{nlsv}}+{\mathbf{nff}}$
• if ${\mathbf{b}}\left(l\right)$ contains the parameter estimate for the intercept then
 $id1l = id2l = id3l = 0 ;$
• if ${\mathbf{b}}\left(l\right)$ contains the parameter estimate for the $i$th level of the $j$th fixed variable, that is the vector of values held in the $k$th column of dat when ${\mathbf{fixed}}\left(j+2\right)=k$ then
 $id1l=0, id2l=j, id3l=i;$
• if the $j$th variable is continuous or binary, that is ${\mathbf{levels}}\left({\mathbf{fixed}}\left(j+2\right)\right)=1$, then ${\mathbf{id}}\left(3,l\right)=0$;
• any remaining rows of the $l$th column of id are set to $0$.
For random effects:
• let
• ${N}_{{R}_{b}}$ denote the number of random variables in the $b$th random statement, that is ${N}_{{R}_{b}}={\mathbf{rndm}}\left(1,b\right)$;
• ${R}_{jb}$ denote the $j$th random variable from the $b$th random statement, that is the vector of values held in the $k$th column of dat when ${\mathbf{rndm}}\left(2+j,b\right)=k$;
• ${N}_{{S}_{b}}$ denote the number of subject variables in the $b$th random statement, that is ${N}_{{S}_{b}}={\mathbf{rndm}}\left(3+{N}_{{R}_{b}},b\right)$;
• ${S}_{jb}$ denote the $j$th subject variable from the $b$th random statement, that is the vector of values held in the $k$th column of dat when ${\mathbf{rndm}}\left(3+{N}_{{R}_{b}}+j,b\right)=k$;
• $L\left({S}_{jb}\right)$ denote the number of levels for ${S}_{jb}$, that is $L\left({S}_{jb}\right)={\mathbf{levels}}\left({\mathbf{rndm}}\left(3+{N}_{{R}_{b}}+j,b\right)\right)$;
• then
• for $l=1,2,\dots {\mathbf{nrf}}×{\mathbf{nlsv}}$, if ${\mathbf{b}}\left(l\right)$ contains the parameter estimate for the $i$th level of ${R}_{jb}$ when ${S}_{\mathit{k}b}={s}_{\mathit{k}}$, for $\mathit{k}=1,2,\dots ,{N}_{{S}_{b}}$ and $1\le {s}_{\mathit{k}}\le L\left({S}_{jb}\right)$, i.e., ${s}_{k}$ is a valid value for the $k$th subject variable, then
 $id1l=b, id2l=j, id3l=i, id3+kl=sk, ​k=1,2,…,NSb;$
• if the parameter being estimated is for the intercept then ${\mathbf{id}}\left(2,l\right)={\mathbf{id}}\left(3,l\right)=0$;
• if the $j$th variable is continuous, or binary, that is $L\left({S}_{jb}\right)=1$, then ${\mathbf{id}}\left(3,l\right)=0$;
• the remaining rows of the $l$th column of id are set to $0$.
In some situations, certain combinations of variables are never observed. In such circumstances all elements of the $l$th row of id are set to $-999$.
7:     $\mathrm{b}\left({\mathbf{lb}}\right)$ – double array
The parameter estimates, with the first ${\mathbf{nrf}}×{\mathbf{nlsv}}$ elements of ${\mathbf{b}}$ containing the parameter estimates for the random effects, $\nu$, and the remaining nff elements containing the parameter estimates for the fixed effects, $\beta$. The order of these estimates are described by the id argument.
8:     $\mathrm{se}\left({\mathbf{lb}}\right)$ – double array
The standard errors of the parameter estimates given in b.
9:     $\mathrm{czz}\left(\mathit{ldczz},:\right)$ – double array
The first dimension of the array czz will be ${\mathbf{nrf}}$.
The second dimension of the array czz will be ${\mathbf{nrf}}×{\mathbf{nlsv}}$.
If ${\mathbf{nlsv}}=1$, then czz holds the lower triangular portion of the matrix $\left(1/{\sigma }^{2}\right)\left({Z}^{\mathrm{T}}{\stackrel{^}{R}}^{-1}Z+{\stackrel{^}{G}}^{-1}\right)$, where $\stackrel{^}{R}$ and $\stackrel{^}{G}$ are the estimates of $R$ and $G$ respectively. If ${\mathbf{nlsv}}>1$ 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 subject variable, the next nrf columns the part corresponding to the second level of the overall subject variable etc.
10:   $\mathrm{cxx}\left(\mathit{ldcxx},:\right)$ – double array
The first dimension of the array cxx will be ${\mathbf{nff}}$.
The second dimension of the array cxx will be ${\mathbf{nff}}$.
cxx holds the lower triangular portion of the matrix $\left(1/{\sigma }^{2}\right){X}^{\mathrm{T}}{\stackrel{^}{V}}^{-1}X$, where $\stackrel{^}{V}$ is the estimated value of $V$.
11:   $\mathrm{cxz}\left(\mathit{ldcxz},:\right)$ – double array
The first dimension of the array cxz will be ${\mathbf{nff}}$.
The second dimension of the array cxz will be ${\mathbf{nlsv}}×{\mathbf{nrf}}$.
If ${\mathbf{nlsv}}=1$, then cxz holds the matrix $\left(1/{\sigma }^{2}\right)\left({X}^{\mathrm{T}}{\stackrel{^}{V}}^{-1}Z\right)\stackrel{^}{G}$, where $\stackrel{^}{V}$ and $\stackrel{^}{G}$ are the estimates of $V$ and $G$ respectively. If ${\mathbf{nlsv}}>1$ then cxz 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 subject variable, the next nrf columns the part corresponding to the second level of the overall subject variable etc.
12:   $\mathrm{ifail}$int64int32nag_int scalar
${\mathbf{ifail}}={\mathbf{0}}$ unless the function detects an error (see Error Indicators and Warnings).

## Error Indicators and Warnings

Errors or warnings detected by the function:

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

${\mathbf{ifail}}=1$
lvpr is too small.
${\mathbf{ifail}}=2$
Constraint: $1\le {\mathbf{vpr}}\left(i\right)\le {\mathbf{nvpr}}$.
${\mathbf{ifail}}=3$
Constraint: $1\le {\mathbf{nvpr}}\le _$.
${\mathbf{ifail}}=4$
Constraint: ${\mathbf{gamma}}\left(1\right)=-1.0$ or ${\mathbf{gamma}}\left(i\right)\ge 0.0$.
${\mathbf{ifail}}=9$
lb is too small.
${\mathbf{ifail}}=11$
ldid is too small.
${\mathbf{ifail}}=15$
ldczz is too small.
${\mathbf{ifail}}=17$
ldcxx is too small.
${\mathbf{ifail}}=19$
ldcxz is too small.
${\mathbf{ifail}}=21$
On entry, icomm has not been initialized correctly.
${\mathbf{ifail}}=32$
On entry, at least one value of i, for $\mathit{i}=1,2,\dots ,{\mathbf{nvpr}}$, does not appear in vpr.
W  ${\mathbf{ifail}}=101$
Optimal solution found, but requested accuracy not achieved.
${\mathbf{ifail}}=102$
Too many major iterations.
W  ${\mathbf{ifail}}=103$
Current point cannot be improved upon.
W  ${\mathbf{ifail}}=104$
At least one negative estimate for gamma was obtained. All negative estimates have been set to zero.
${\mathbf{ifail}}=-99$
${\mathbf{ifail}}=-399$
Your licence key may have expired or may not have been installed correctly.
${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.

## Accuracy

Not applicable.

The argument vpr gives the mapping between the random variables and the variance components. In most cases ${\mathbf{vpr}}\left(\mathit{i}\right)=\mathit{i}$, for $\mathit{i}=1,2,\dots ,{\sum }_{\mathit{i}}{\mathbf{rndm}}\left(1,\mathit{i}\right)+{\mathbf{rndm}}\left(2,\mathit{i}\right)$. However, in some cases it might be necessary to associate more than one random variable with a single variance component, for example, when the columns of dat hold dummy variables.
Consider a dataset with three variables:
 $dat= 113.6 214.5 311.1 128.3 227.2 326.1$
where the first column corresponds to a categorical variable with three levels, the next to a categorical variable with two levels and the last column to a continuous variable. So in a call to nag_correg_mixeff_hier_init (g02jc)
 $levels=321$
also assume a model with no fixed effects, no random intercept, no nesting and all three variables being included as random effects, then
 $fixed=00; rndm=30123T.$
Each of the three columns in dat therefore correspond to a single variable and hence there are three variance components, one for each random variable included in the model, so
 $vpr=123.$
This is the recommended way of supplying the data to nag_correg_mixeff_hier_ml (g02je), however it is possible to reformat the above dataset by replacing each of the categorical variables with a series of dummy variables, one for each level. The dataset then becomes
 $dat= 100103.6 010104.5 001101.1 100018.3 010017.2 001016.1$
where each column only has one level
 $levels= 111111 .$
Again a model with no fixed effects, no random intercept, no nesting and all variables being included as random effects is required, so
 $fixed=00 ; rndm= 60123456T .$
With the data entered in this manner, the first three columns of dat correspond to a single variable (the first column of the original dataset) as do the next two columns (the second column of the original dataset). Therefore vpr must reflect this
 $vpr= 111223 .$
In most situations it is more efficient to supply the data to nag_correg_mixeff_hier_init (g02jc) in terms of categorical variables rather than transform them into dummy variables.

## Example

This example fits a random effects model with three levels of nesting to a simulated dataset with $90$ observations and $12$ variables.
```function g02je_example

fprintf('g02je example results\n\n');

n = 90;
ncol = 12;
nrndm = 3;
nvpr = int64(7);

% The number of levels associated with each of the independent variables
levels = [int64(2); 3; 2; 3; 2; 3; 1; 4; 5; 2; 3; 3];

% The Fixed part of the model
fixed = [int64(2); 1; 1; 2];

% The Random part of the model
rndm = [int64(2), 2,  3;
0,  0,  0;
3,  5,  7;
4,  6,  8;
3,  2,  9;
10, 11,  1;
11, 12, 12;
12,  0,  0];

% Dependant data
y = [ 3.1100;   2.8226;   7.4543;   4.4313;   6.1543;  -0.1783;
4.6748;   7.0667;   1.4262;   7.7290;  -2.1806;   6.8419;
1.2590;   8.8405;   6.1657;  -4.5605;  -1.2367; -12.2932;
-2.3374;   0.0716;   0.1895;   1.5608;  -0.8529;  -4.1169;
3.9977;  -8.1277;  -4.9656;  -0.6428;  -5.5152;  -5.5657;
14.8177;  16.9783;  13.8966;  14.8166;  19.3640;   9.5299;
12.0102;   6.1551;  -1.7048;   2.7640;   2.8065;   0.0974;
-7.8080; -18.0450;  -2.8199;   8.9893;   3.7978;  -6.3493;
8.1411;  -7.5483;  -0.4600;  -3.2135;  -6.6562;   5.1267;
3.5592;  -4.4420;  -8.5965;  -6.3187;  -7.8953; -10.1383;
-7.8850;  23.2001;   5.5829;  -4.3698;   2.1274;  -2.7184;
-17.9128;  -1.2708; -24.2735; -14.7374;   0.1713;   8.0006;
1.2100;   3.3307; -22.6713;   7.5562;  -7.0694;   3.7159;
-4.3135; -14.5577; -12.5107;   4.7708;  13.2797;  -6.3243;
-7.0549;  -9.2713; -18.7788;  -7.7230; -22.7230; -11.6609];

% Independent data
dat = [1, 3, 2, 1, 2, 2,-0.3160, 4, 2, 1, 1, 1;
1, 1, 1, 3, 1, 2,-1.3377, 1, 4, 1, 1, 1;
1, 3, 1, 3, 1, 3,-0.7610, 4, 2, 1, 1, 1;
2, 3, 2, 1, 1, 3,-2.2976, 4, 2, 1, 1, 1;
2, 2, 1, 3, 2, 3,-0.4263, 2, 1, 1, 1, 1;
2, 1, 2, 3, 1, 3, 1.4067, 3, 3, 2, 1, 1;
2, 3, 2, 1, 2, 1,-1.4669, 1, 2, 2, 1, 1;
1, 1, 1, 3, 2, 3, 0.4717, 2, 4, 2, 1, 1;
1, 3, 2, 3, 2, 1, 0.4436, 1, 3, 2, 1, 1;
1, 1, 1, 2, 2, 3,-0.5950, 3, 4, 2, 1, 1;
1, 3, 1, 3, 1, 1,-1.7981, 4, 2, 1, 2, 1;
2, 3, 1, 2, 1, 1, 0.2397, 1, 4, 1, 2, 1;
1, 2, 2, 1, 2, 3, 0.4742, 1, 1, 1, 2, 1;
2, 2, 2, 2, 2, 3, 0.6888, 3, 1, 1, 2, 1;
2, 1, 2, 3, 1, 3,-1.0616, 3, 5, 1, 2, 1;
1, 2, 2, 2, 2, 1,-0.5356, 1, 3, 2, 2, 1;
1, 3, 2, 2, 1, 1,-1.2963, 2, 5, 2, 2, 1;
1, 2, 2, 1, 2, 2,-1.5389, 3, 2, 2, 2, 1;
2, 3, 1, 1, 2, 2,-0.6408, 2, 1, 2, 2, 1;
1, 2, 2, 2, 1, 1, 0.6574, 1, 1, 2, 2, 1;
2, 1, 1, 1, 1, 3, 0.9259, 1, 2, 1, 3, 1;
2, 2, 2, 1, 2, 2, 1.5080, 3, 1, 1, 3, 1;
2, 3, 1, 1, 1, 3, 2.5821, 2, 3, 1, 3, 1;
1, 2, 2, 1, 2, 3, 0.4102, 1, 4, 1, 3, 1;
2, 1, 2, 3, 2, 2, 0.7839, 2, 5, 1, 3, 1;
1, 2, 2, 3, 2, 1,-1.8812, 4, 2, 2, 3, 1;
1, 2, 1, 3, 2, 3, 0.7770, 4, 1, 2, 3, 1;
2, 2, 1, 2, 1, 3, 0.2590, 3, 1, 2, 3, 1;
2, 3, 2, 2, 2, 3,-0.9250, 3, 3, 2, 3, 1;
2, 2, 1, 3, 2, 3,-0.4831, 1, 5, 2, 3, 1;
2, 2, 1, 3, 1, 3, 0.5046, 3, 3, 1, 1, 2;
2, 1, 1, 2, 2, 1,-0.6903, 2, 1, 1, 1, 2;
1, 3, 2, 2, 2, 1, 1.6166, 2, 5, 1, 1, 2;
2, 2, 2, 2, 1, 3, 0.2778, 2, 3, 1, 1, 2;
2, 3, 2, 2, 1, 2, 1.9586, 4, 2, 1, 1, 2;
1, 3, 1, 1, 1, 3, 1.0506, 2, 5, 2, 1, 2;
2, 1, 1, 3, 2, 3, 0.4871, 1, 1, 2, 1, 2;
2, 1, 2, 3, 2, 1, 2.0891, 4, 4, 2, 1, 2;
1, 2, 1, 1, 2, 2, 1.4338, 4, 3, 2, 1, 2;
1, 1, 2, 3, 1, 2,-1.1196, 3, 4, 2, 1, 2;
1, 3, 1, 1, 2, 1, 0.3367, 3, 2, 1, 2, 2;
2, 2, 1, 3, 1, 1, 0.1092, 2, 2, 1, 2, 2;
1, 1, 1, 2, 2, 2, 0.4007, 4, 1, 1, 2, 2;
2, 3, 1, 1, 1, 2, 0.1460, 3, 5, 1, 2, 2;
2, 1, 2, 3, 1, 3,-0.3877, 3, 4, 1, 2, 2;
1, 1, 1, 2, 2, 1, 0.6957, 4, 3, 2, 2, 2;
2, 1, 1, 1, 2, 1,-0.4664, 3, 3, 2, 2, 2;
1, 1, 1, 1, 2, 3, 0.2067, 2, 4, 2, 2, 2;
2, 1, 2, 1, 1, 2, 0.4112, 1, 4, 2, 2, 2;
2, 2, 1, 1, 1, 2,-1.3734, 3, 3, 2, 2, 2;
2, 1, 2, 3, 1, 3, 0.7065, 1, 3, 1, 3, 2;
1, 2, 2, 2, 1, 2, 1.3628, 4, 2, 1, 3, 2;
2, 1, 2, 2, 2, 3,-0.5052, 4, 5, 1, 3, 2;
2, 1, 1, 1, 2, 1,-1.3457, 2, 5, 1, 3, 2;
1, 1, 2, 1, 2, 3,-1.8022, 3, 4, 1, 3, 2;
2, 3, 1, 2, 1, 1, 0.0116, 2, 4, 2, 3, 2;
2, 2, 1, 3, 2, 3,-0.9075, 1, 3, 2, 3, 2;
2, 2, 2, 2, 2, 3,-1.4707, 1, 1, 2, 3, 2;
2, 2, 1, 1, 2, 1,-1.2938, 2, 3, 2, 3, 2;
1, 3, 1, 3, 2, 2,-1.1660, 4, 4, 2, 3, 2;
1, 2, 1, 1, 2, 3, 0.0397, 4, 4, 1, 1, 3;
1, 3, 1, 2, 1, 3,-0.5987, 3, 2, 1, 1, 3;
2, 3, 2, 2, 1, 1, 0.6683, 3, 3, 1, 1, 3;
2, 2, 1, 1, 2, 2,-0.0106, 1, 3, 1, 1, 3;
1, 2, 1, 3, 2, 2, 0.5885, 1, 3, 1, 1, 3;
1, 1, 1, 1, 1, 2, 0.4555, 1, 5, 2, 1, 3;
2, 2, 2, 1, 1, 2, 0.6502, 4, 3, 2, 1, 3;
1, 1, 1, 3, 1, 1,-0.1601, 1, 3, 2, 1, 3;
2, 2, 1, 3, 2, 3, 1.6910, 1, 1, 2, 1, 3;
2, 2, 2, 3, 1, 2, 0.1053, 4, 4, 2, 1, 3;
2, 1, 2, 3, 2, 2,-0.4037, 3, 4, 1, 2, 3;
1, 3, 2, 3, 1, 3,-0.5853, 3, 2, 1, 2, 3;
2, 3, 2, 1, 1, 1,-0.3037, 1, 3, 1, 2, 3;
1, 3, 1, 1, 2, 2,-0.0774, 1, 4, 1, 2, 3;
2, 3, 1, 2, 2, 1, 0.4733, 4, 5, 1, 2, 3;
1, 3, 2, 2, 1, 2,-0.0354, 4, 2, 2, 2, 3;
1, 3, 2, 2, 1, 1,-0.6640, 2, 1, 2, 2, 3;
2, 3, 1, 3, 1, 1, 0.0335, 4, 4, 2, 2, 3;
1, 2, 2, 2, 1, 3, 0.1351, 1, 1, 2, 2, 3;
1, 1, 2, 1, 2, 3,-0.5951, 3, 4, 2, 2, 3;
2, 2, 2, 3, 1, 3, 0.2735, 3, 2, 1, 3, 3;
2, 2, 1, 1, 1, 3, 0.3157, 1, 2, 1, 3, 3;
2, 2, 2, 1, 1, 1,-1.0843, 2, 3, 1, 3, 3;
1, 2, 2, 1, 2, 2,-0.0836, 4, 2, 1, 3, 3;
2, 1, 2, 1, 1, 2,-0.2884, 2, 1, 1, 3, 3;
2, 3, 2, 3, 2, 3,-0.1006, 1, 2, 2, 3, 3;
1, 3, 1, 2, 2, 3, 0.5710, 1, 3, 2, 3, 3;
1, 1, 2, 1, 1, 2, 0.2776, 2, 3, 2, 3, 3;
2, 3, 2, 2, 1, 3,-0.7561, 4, 4, 2, 3, 3;
1, 2, 2, 2, 1, 2, 1.5549, 1, 4, 2, 3, 3];

vpr = [int64(1):7];
gamma = zeros(8, 1);
gamma(1) = -1; % Estimate initial values for the variance components

% Call the initialisation routine once to get lrcomm and licomm
lrcomm = int64(0);
licomm = int64(2);
[nff, nlsv, nrf, rcomm, icomm, ifail] = ...
g02jc( ...
dat, levels, y, fixed, rndm, lrcomm, licomm);
licomm = icomm(1);
lrcomm = icomm(2);

% Pre-process the data
[nff, nlsv, nrf, rcomm, icomm, ifail] = ...
g02jc( ...
dat, levels, y, fixed, rndm, lrcomm, licomm);

% Use default options
iopt = zeros(0, 0, 'int64');
ropt = zeros(0, 0);

lb = nff + nrf*nlsv;

% Perform the analysis
[gamma, effn, rnkx, ncov, lnlike, id, b, se, czz, cxx, cxz, ifail] = ...
g02je( ...
vpr, nvpr, gamma, lb, rcomm, icomm, iopt, ropt);

% Print results
fprintf('Number of observations (n)                    = %d\n', n);
fprintf('Number of random factors (nrf)                = %d\n', nrf);
fprintf('Number of fixed factors (nff)                 = %d\n', nff);
fprintf('Number of subject levels (nlsv)               = %d\n', nlsv);
fprintf('Rank of X (rnkx)                              = %d\n', rnkx);
fprintf('Effective n (effn)                            = %d\n', effn);
fprintf('Number of non-zero variance components (ncov) = %d\n', ncov);

fprintf('\nParameter Estimates\n');

if nrf > 0
fprintf('\nRandom Effects\n');
end
pb = -999;
pfmt = ' ';
for k = 1:nrf*nlsv
tb = id(1,k);
if tb ~= -999
vid = id(2,k);
nv = rndm(1,tb);
ns = rndm(3+nv,tb);
tfmt = sprintf('%d ', id(3+1:3+ns,k));
if ( (pb ~= tb) || (strcmp(tfmt, pfmt) == 0) )
if (k ~= 1)
fprintf('\n');
end
fprintf('  Subject: ');
for l=1:ns
fprintf(' Variable %2d (Level %2d)',rndm(3+nv+l,tb), id(3+l,k));
end
fprintf('\n');
end
if (vid==0)
% Intercept
fprintf('    Intercept%18s%10.4f %10.4f\n', ' ', b(k), se(k));
else
% variable vid specified in rndm
aid = rndm(2+vid,tb);
if (id(3,k)==0)
fprintf('    Variable %2d%16s%10.4f %10.4f\n', aid, ' ', b(k), se(k));
else
fprintf('    Variable %2d (Level %2d)     %10.4f %10.4f\n', ...
aid, id(3,k), b(k), se(k));
end
end
pfmt = tfmt;
end
pb = tb;
end

if nff>0
fprintf('\nFixed Effects\n');
end
for k = (nrf*nlsv+1):(nrf*nlsv+nff)
if vid~=-999
vid = id(2,k);
if vid==0
% Intercept
fprintf('    Intercept%18s%10.4f %10.4f\n', ' ', b(k), se(k));
else
% vid'th variable specified in fixed
aid = fixed(2+vid);
if (id(3,k)==0)
fprintf('    Variable %2d%16s%10.4f %10.4f\n', aid, ' ', b(k), se(k));
else
fprintf('    Variable %2d (Level %2d)     %10.4f %10.4f\n', ...
aid, id(3,k), b(k), se(k));
end
end
end
end

fprintf('\nVariance Components\n');
fprintf('  Estimate     Parameter       Subject\n');
for k = 1:nvpr
fprintf('%10.5f     ', gamma(k));
p = 0;
for tb = 1:nrndm
nv = rndm(1,tb);
ns = rndm(3+nv,tb);
if (rndm(2,tb)==1)
p = p + 1;
if (vpr(p)==k)
fprintf('Intercept       Variables ');
fprintf('%2d ', rndm(3+nv+1:3+nv+ns, tb));
fprintf('\n');
end
end
for i = 1:nv
p = p + 1;
if (vpr(p)==k)
fprintf('Variable %2d     Variables %2d ', rndm(2+i,tb));
fprintf('%2d ', rndm(3+nv+1:3+nv+ns, tb));
end
end
end
fprintf('\n');
end

fprintf('\nsigma^2          = %15.5f\n', gamma(nvpr+1));
fprintf('-2log likelihood = %15.5f\n', lnlike);

```
```g02je example results

Number of observations (n)                    = 90
Number of random factors (nrf)                = 55
Number of fixed factors (nff)                 = 4
Number of subject levels (nlsv)               = 3
Rank of X (rnkx)                              = 4
Effective n (effn)                            = 90
Number of non-zero variance components (ncov) = 7

Parameter Estimates

Random Effects
Subject:  Variable 10 (Level  1) Variable 11 (Level  1) Variable 12 (Level  1)
Variable  3 (Level  1)         2.1566     3.7320
Variable  3 (Level  2)         1.7769     3.8543
Variable  4 (Level  1)         0.5583     3.0508

Subject:  Variable 10 (Level  1) Variable 11 (Level  1) Variable 12 (Level  1)
Variable  4 (Level  3)         0.6776     3.0358

Subject:  Variable 10 (Level  2) Variable 11 (Level  1) Variable 12 (Level  1)
Variable  3 (Level  1)         1.4448     3.3293
Variable  3 (Level  2)        -2.8634     3.3533
Variable  4 (Level  1)         3.6811     2.2253
Variable  4 (Level  2)        -1.9988     2.2929
Variable  4 (Level  3)        -2.1281     1.9896

Subject:  Variable 10 (Level  1) Variable 11 (Level  2) Variable 12 (Level  1)
Variable  3 (Level  1)        -3.1562     3.8624
Variable  3 (Level  2)         2.8856     4.6985
Variable  4 (Level  1)        -4.6811     2.2236
Variable  4 (Level  2)         5.5794     2.1390
Variable  4 (Level  3)        -0.9832     2.2841

Subject:  Variable 10 (Level  2) Variable 11 (Level  2) Variable 12 (Level  1)
Variable  3 (Level  1)         4.3449     3.6258
Variable  3 (Level  2)        -4.4285     3.4096
Variable  4 (Level  1)        -1.0798     3.1008
Variable  4 (Level  2)         1.0536     2.9612

Subject:  Variable 10 (Level  1) Variable 11 (Level  3) Variable 12 (Level  1)
Variable  3 (Level  1)         0.4216     4.0146
Variable  3 (Level  2)         0.2268     3.4265
Variable  4 (Level  1)        -1.0626     2.3505

Subject:  Variable 10 (Level  1) Variable 11 (Level  3) Variable 12 (Level  1)
Variable  4 (Level  3)         1.2664     2.5276

Subject:  Variable 10 (Level  2) Variable 11 (Level  3) Variable 12 (Level  1)
Variable  3 (Level  1)         1.2785     3.4331
Variable  3 (Level  2)        -1.6652     3.8605

Subject:  Variable 10 (Level  2) Variable 11 (Level  3) Variable 12 (Level  1)
Variable  4 (Level  2)         0.7332     2.6958
Variable  4 (Level  3)        -0.8547     2.7819

Subject:  Variable 11 (Level  1) Variable 12 (Level  1)
Variable  5 (Level  1)        -0.5540     2.8120
Variable  5 (Level  2)         1.9179     2.7500
Variable  6 (Level  1)         0.6925     3.6813
Variable  6 (Level  2)        -2.2632     3.1202
Variable  6 (Level  3)         4.3216     3.1131

Subject:  Variable 11 (Level  2) Variable 12 (Level  1)
Variable  5 (Level  1)         1.5151     2.9154
Variable  5 (Level  2)        -1.7072     2.8715
Variable  6 (Level  1)         0.2154     3.9398
Variable  6 (Level  2)        -3.7591     4.2153
Variable  6 (Level  3)         3.1563     4.7621

Subject:  Variable 11 (Level  3) Variable 12 (Level  1)
Variable  5 (Level  1)         1.7892     3.1214
Variable  5 (Level  2)        -1.6473     3.1579
Variable  6 (Level  1)        -1.2268     3.8853
Variable  6 (Level  2)         4.6247     3.6412
Variable  6 (Level  3)        -3.1117     3.1648

Subject:  Variable 12 (Level  1)
Variable  7                    0.6016     0.4634
Variable  8 (Level  1)         1.5887     1.2518
Variable  8 (Level  2)        -0.7951     1.4856
Variable  8 (Level  3)         0.3798     1.6037
Variable  8 (Level  4)        -0.8295     1.6629
Variable  9 (Level  1)         0.5197     1.5510
Variable  9 (Level  2)         0.0156     1.8248
Variable  9 (Level  3)        -0.1723     1.8271
Variable  9 (Level  4)         0.4305     1.9494
Variable  9 (Level  5)        -0.1412     2.0379

Subject:  Variable 10 (Level  1) Variable 11 (Level  1) Variable 12 (Level  2)
Variable  3 (Level  1)         6.3424     3.3173
Variable  3 (Level  2)         5.7538     3.3626

Subject:  Variable 10 (Level  1) Variable 11 (Level  1) Variable 12 (Level  2)
Variable  4 (Level  2)         2.5053     2.6520
Variable  4 (Level  3)         1.2953     2.6978

Subject:  Variable 10 (Level  2) Variable 11 (Level  1) Variable 12 (Level  2)
Variable  3 (Level  1)         1.6342     3.7874
Variable  3 (Level  2)        -2.8693     3.8549
Variable  4 (Level  1)        -0.9274     2.7266

Subject:  Variable 10 (Level  2) Variable 11 (Level  1) Variable 12 (Level  2)
Variable  4 (Level  3)         0.5394     2.7100

Subject:  Variable 10 (Level  1) Variable 11 (Level  2) Variable 12 (Level  2)
Variable  3 (Level  1)       -10.2379     3.2977
Variable  3 (Level  2)         3.2457     4.0593
Variable  4 (Level  1)        -2.8362     2.2599
Variable  4 (Level  2)         0.2805     2.9513
Variable  4 (Level  3)         0.3587     2.8663

Subject:  Variable 10 (Level  2) Variable 11 (Level  2) Variable 12 (Level  2)
Variable  3 (Level  1)        -1.3161     3.1545
Variable  3 (Level  2)         8.2719     3.9322
Variable  4 (Level  1)        -0.4813     2.3705
Variable  4 (Level  2)         2.6668     2.4832

Subject:  Variable 10 (Level  1) Variable 11 (Level  3) Variable 12 (Level  2)
Variable  3 (Level  1)         4.9485     3.9465
Variable  3 (Level  2)         0.0987     3.5531
Variable  4 (Level  1)         3.0791     2.1790
Variable  4 (Level  2)        -1.9469     2.3796
Variable  4 (Level  3)         0.4536     2.1984

Subject:  Variable 10 (Level  2) Variable 11 (Level  3) Variable 12 (Level  2)
Variable  3 (Level  1)        -4.5419     3.2940
Variable  3 (Level  2)        -3.9095     4.0163
Variable  4 (Level  1)        -0.4456     2.6194
Variable  4 (Level  2)        -1.5462     2.6514
Variable  4 (Level  3)        -0.6636     2.8738

Subject:  Variable 11 (Level  1) Variable 12 (Level  2)
Variable  5 (Level  1)         4.9921     3.0570
Variable  5 (Level  2)         0.8986     3.0576
Variable  6 (Level  1)         7.0091     3.7851
Variable  6 (Level  2)        -1.3173     3.1348
Variable  6 (Level  3)         6.1881     3.4928

Subject:  Variable 11 (Level  2) Variable 12 (Level  2)
Variable  5 (Level  1)        -0.3947     3.0751
Variable  5 (Level  2)         0.3750     3.0579
Variable  6 (Level  1)         6.9902     3.2654
Variable  6 (Level  2)        -1.0683     3.5699
Variable  6 (Level  3)        -5.9617     3.6688

Subject:  Variable 11 (Level  3) Variable 12 (Level  2)
Variable  5 (Level  1)        -1.0471     3.0732
Variable  5 (Level  2)        -0.7991     2.9597
Variable  6 (Level  1)         2.7549     3.8142
Variable  6 (Level  2)        -6.3441     3.2624
Variable  6 (Level  3)        -0.1341     3.5956

Subject:  Variable 12 (Level  2)
Variable  7                    0.1533     0.5196
Variable  8 (Level  1)         1.6630     1.8224
Variable  8 (Level  2)        -0.6835     1.6502
Variable  8 (Level  3)        -0.0959     1.5604
Variable  8 (Level  4)         0.1696     1.4537
Variable  9 (Level  1)         1.0203     2.2901
Variable  9 (Level  2)         6.4354     1.7420
Variable  9 (Level  3)        -1.5942     1.7761
Variable  9 (Level  4)         0.0955     1.9436
Variable  9 (Level  5)        -3.9588     1.7124

Subject:  Variable 10 (Level  1) Variable 11 (Level  1) Variable 12 (Level  3)
Variable  3 (Level  1)        10.9751     3.2085
Variable  3 (Level  2)        -1.0674     3.7219
Variable  4 (Level  1)        -2.8350     2.2037
Variable  4 (Level  2)         3.7075     2.7912
Variable  4 (Level  3)         2.2405     2.2796

Subject:  Variable 10 (Level  2) Variable 11 (Level  1) Variable 12 (Level  3)
Variable  3 (Level  1)        -6.2719     3.3190
Variable  3 (Level  2)        -9.2923     3.7884
Variable  4 (Level  1)        -2.8586     2.3728

Subject:  Variable 10 (Level  2) Variable 11 (Level  1) Variable 12 (Level  3)
Variable  4 (Level  3)        -2.0316     2.2895

Subject:  Variable 10 (Level  1) Variable 11 (Level  2) Variable 12 (Level  3)
Variable  3 (Level  1)        -3.3222     3.4246
Variable  3 (Level  2)        -0.3111     3.2221
Variable  4 (Level  1)         1.6131     2.3970
Variable  4 (Level  2)        -3.0099     2.9300
Variable  4 (Level  3)         0.2552     2.7229

Subject:  Variable 10 (Level  2) Variable 11 (Level  2) Variable 12 (Level  3)
Variable  3 (Level  1)         6.6372     3.9751
Variable  3 (Level  2)        -5.4249     3.4039
Variable  4 (Level  1)        -3.2357     2.8565
Variable  4 (Level  2)         1.5313     2.8232
Variable  4 (Level  3)         2.0854     3.0661

Subject:  Variable 10 (Level  1) Variable 11 (Level  3) Variable 12 (Level  3)
Variable  3 (Level  1)         8.5902     4.0894
Variable  3 (Level  2)        -1.6058     3.2906
Variable  4 (Level  1)         3.2575     2.5450

Subject:  Variable 10 (Level  1) Variable 11 (Level  3) Variable 12 (Level  3)
Variable  4 (Level  3)        -1.0630     2.8692

Subject:  Variable 10 (Level  2) Variable 11 (Level  3) Variable 12 (Level  3)
Variable  3 (Level  1)        -4.5747     3.9475
Variable  3 (Level  2)        -4.1752     3.0911
Variable  4 (Level  1)         1.0578     2.5496
Variable  4 (Level  2)        -4.4284     2.2029
Variable  4 (Level  3)         0.6214     2.5884

Subject:  Variable 11 (Level  1) Variable 12 (Level  3)
Variable  5 (Level  1)         5.4387     3.0091
Variable  5 (Level  2)        -8.5065     3.1099
Variable  6 (Level  1)        -0.9179     3.7257
Variable  6 (Level  2)        -2.4920     3.1176
Variable  6 (Level  3)        -2.7772     3.4083

Subject:  Variable 11 (Level  2) Variable 12 (Level  3)
Variable  5 (Level  1)         4.4193     3.1282
Variable  5 (Level  2)        -5.7324     3.1435
Variable  6 (Level  1)        -5.9992     3.1431
Variable  6 (Level  2)         5.5657     3.2599
Variable  6 (Level  3)        -2.2147     3.1758

Subject:  Variable 11 (Level  3) Variable 12 (Level  3)
Variable  5 (Level  1)         0.3594     2.9017
Variable  5 (Level  2)        -1.3169     3.0004
Variable  6 (Level  1)        14.5815     3.8519
Variable  6 (Level  2)        -5.2262     3.2578
Variable  6 (Level  3)       -11.2864     3.1821

Subject:  Variable 12 (Level  3)
Variable  7                   -0.2970     0.5930
Variable  8 (Level  1)         2.6255     1.5201
Variable  8 (Level  2)         0.5048     1.7865
Variable  8 (Level  3)        -0.1518     1.8905
Variable  8 (Level  4)        -4.3754     1.4651
Variable  9 (Level  1)        -4.4219     2.0532
Variable  9 (Level  2)         3.7058     1.9085
Variable  9 (Level  3)        -1.7524     1.7894
Variable  9 (Level  4)         0.4339     1.8210
Variable  9 (Level  5)        -0.6161     2.3700

Fixed Effects
Intercept                      1.5913     2.4106
Variable  1 (Level  2)        -1.5994     0.8183
Variable  2 (Level  2)        -2.3793     1.0996
Variable  2 (Level  3)         0.5328     1.1677

Variance Components
Estimate     Parameter       Subject
36.38867     Variable  3     Variables 10 11 12
11.43322     Variable  4     Variables 10 11 12
19.73586     Variable  5     Variables 11 12
39.80174     Variable  6     Variables 11 12
0.41583     Variable  7     Variables 12
5.16442     Variable  8     Variables 12
9.79904     Variable  9     Variables 12

sigma^2          =         0.00042
-2log likelihood =       617.11969
```