hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
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ν + ε
y=Xβ+Zν+ε
where yy is a vector of nn observations on the dependent variable,
XX is a known nn by pp design matrix for the fixed independent variables,
ββ is a vector of length pp of unknown fixed effects,
ZZ is a known nn by qq design matrix for the random independent variables,
νν is a vector of length qq of unknown random effects,
and εε is a vector of length nn of unknown random errors.
Both ν ν  and ε ε  are assumed to have a Gaussian distribution with expectation zero and variance/covariance matrix defined by
Var
[ ν ε ]
=
[ G 0 0 R ]
Var[ ν ε ] = [ G 0 0 R ]
where R = σR2 I R= σ R 2 I , I I  is the n × n n×n  identity matrix and G G  is a diagonal matrix. It is assumed that the random variables, Z Z , can be subdivided into g q g q  groups with each group being identically distributed with expectation zero and variance σi2 σi2 . The diagonal elements of matrix G G  therefore take one of the values {σi2 : i = 1,2,,g} { σi2 : i=1,2,,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 g+1  variance components γ γ , where γ = {σ12,σ22,,σg12,σg2,σR2} γ = {σ12,σ22,, σ g-1 2 ,σg2,σR2} . Rather than working directly with γ γ , nag_correg_mixeff_hier_ml (g02je) uses an iterative process to estimate γ* = { σ12 / σR2 , σ22 / σR2 ,, σg12 / σR2 , σg2 / σR2 ,1} γ* = { σ12 / σR2 , σ22 / σR2 ,, σg-12 / σR2 , σg2 / σR2 ,1} . Due to the iterative nature of the estimation a set of initial values, γ0 γ0 , for γ* γ*  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 lR = log(|V|) + n log(rTV1r) + log(2π / n)
-2 l R = log( |V| ) + n log( rT V-1 r ) + log( 2 π / n )
where
V = ZG ZT + R,   r = yXb   and   b = (XTV1X)1 XT V1 y .
V = ZG ZT + R,   r=y-Xb   and   b = ( XT V-1 X )-1 XT V-1 y .
Once the final estimates for γ* γ *  have been obtained, the value of σR2 σR2  is given by
σR2 = (rTV1r) / (np) .
σR2 = ( rT V-1 r ) / ( n - p ) .
Case weights, Wc Wc , can be incorporated into the model by replacing XTX XTX  and ZTZ ZTZ  with XTWcX XTWcX  and ZTWcZ ZTWcZ  respectively, for a diagonal weight matrix Wc Wc .
The log-likelihood, lRlR, 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 parameter 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:     vpr(lvpr) – int64int32nag_int array
lvpr, the dimension of the array, must satisfy the constraint lvpr = iRNDM(1,i) + RNDM(2,i)lvpr=iRNDM(1,i)+RNDM(2,i).
A vector of flags indicating the mapping between the random variables specified in rndm and the variance components, σi2σi2. See Section [Further Comments] for more details.
Constraint: 1vpr(i)nvpr1vprinvpr, for i = 1,2,,lvpri=1,2,,lvpr.
2:     nvpr – int64int32nag_int scalar
gg, the number of variance components being estimated (excluding the overall variance, σR2σR2).
Constraint: 1nvprlvpr 1nvprlvpr .
3:     gamma(nvpr + 1nvpr+1) – double array
Holds the initial values of the variance components, γ0 γ0 , with gamma(i)gammai the initial value for σi2 / σR2σi2/σR2, for i = 1,2,,nvpri=1,2,,nvpr.
If gamma(1) = 1.0gamma1=-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: gamma(1) = 1.0 ​ or ​ gamma(i)0.0gamma1=-1.0 ​ or ​ gammai0.0, for i = 1,2,,gi=1,2,,g.
4:     lb – int64int32nag_int scalar
The dimension of the arrays b and se and the second dimension of the array id as declared in the (sub)program from which nag_correg_mixeff_hier_ml (g02je) is called.
Constraint: lbnff + nrf × nlsvlbnff+nrf×nlsv.
5:     rcomm( : :) – double array
Note: the dimension of the array rcomm must be at least lrcommlrcomm (see nag_correg_mixeff_hier_init (g02jc)).
Communication array initialized by a call to nag_correg_mixeff_hier_init (g02jc).
6:     icomm( : :) – int64int32nag_int array
Note: the dimension of the array icomm must be at least licommlicomm (see nag_correg_mixeff_hier_init (g02jc)).
Communication array initialized by a call to nag_correg_mixeff_hier_init (g02jc).
7:     iopt(liopt) – 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 (e04uc) can be chosen by setting iopt(5) = 1iopt5=1. If liopt < 4liopt<4 or iopt(5)1iopt51 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-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 11.
nag_opt_bounds_mod_deriv2_comp (e04lb) is being used


ii


Description
Equivalent
nag_opt_bounds_mod_deriv2_comp (e04lb)
parameter


Default Value
11 Number of iterations maxcal 10001000
22 Unit number for monitoring information n/a As returned by nag_file_set_unit_advisory (x04ab)
33 Print optional parameters (1 = 1= print) n/a 1-1 (no printing performed)
44 Frequency that monitoring information is printed iprint 1-1
55 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).
(e04uc) is being used


ii


Description
Equivalent
(e04uc)
parameter


Default Value
11 Number of iterations Major Iteration Limit max(50,3 × nvpr)max(50,3×nvpr)
22 Unit number for monitoring information n/a As returned by nag_file_set_unit_advisory (x04ab)
33 Print optional parameters (1 = 1= print, otherwise no print) List/Nolist 1-1 (no printing performed)
44 Frequency that monitoring information is printed Major Print Level 00
55 Optimizer used n/a n/a
66 Number of minor iterations Minor Iteration Limit max(50,3 × nvpr)max(50,3×nvpr)
77 Frequency that additional monitoring information is printed Minor Print Level 00
8:     ropt(lropt) – 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-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 11.
nag_opt_bounds_mod_deriv2_comp (e04lb) is being used


ii


Description
Equivalent
nag_opt_bounds_mod_deriv2_comp (e04lb)
parameter


Default Value
11 Sweep tolerance n/a max(sqrt(eps),sqrt(eps) × maxi (zzii))max(eps,eps×maxi(zzii))
22 Accuracy of linear minimizations eta 0.90.9
33 Accuracy to which solution is required xtol 0.00.0
44 Initial distance from solution stepmx 100000.0100000.0
(e04uc) is being used


ii


Description
Equivalent
(e04uc)
parameter


Default Value
11 Sweep tolerance n/a max(sqrt(eps),sqrt(eps) × maxi (zzii))max(eps,eps×maxi(zzii))
22 Lower bound for γ*γ* n/a eps / 100eps/100
33 Upper bound for γ*γ* n/a 10201020
44 Line search tolerance Line Search Tolerance 0.90.9
55 Optimality tolerance Optimality Tolerance eps0.72eps0.72
where epseps is the machine precision returned by nag_machine_precision (x02aj) and zziizzii denotes the ii diagonal element of ZTZZTZ.

Optional Input Parameters

1:     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: lvpr = iRNDM(1,i) + RNDM(2,i)lvpr=iRNDM(1,i)+RNDM(2,i).
2:     liopt – int64int32nag_int scalar
Default: The dimension of the array iopt.
Length of the options array iopt. If liopt0liopt0 then iopt is not referenced and default values are used for all optional parameters.
3:     lropt – int64int32nag_int scalar
Default: The dimension of the array ropt.
Length of the options array ropt. If lropt0lropt0 then ropt is not referenced and default values are used for all optional parameters.

Input Parameters Omitted from the MATLAB Interface

ldid ldczz ldcxx ldcxz

Output Parameters

1:     gamma(nvpr + 1nvpr+1) – double array
gamma(i)gammai, for i = 1,2,,nvpri=1,2,,nvpr, holds the final estimate of σi2σi2 and gamma(nvpr + 1)gammanvpr+1 holds the final estimate for σR2σR2.
2:     effn – int64int32nag_int scalar
Effective number of observations. If there are no weights (i.e., weight = 'U'weight='U'), or all weights are nonzero, then effn = neffn=n.
3:     rnkx – int64int32nag_int scalar
The rank of the design matrix, XX, for the fixed effects.
4:     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 ncov = nvprncov=nvpr.
5:     lnlike – double scalar
2 lR (γ̂) - 2 lR (γ^) where lR lR is the log of the maximum likelihood calculated at γ̂ γ^ , the estimated variance components returned in gamma.
6:     id(ldid,lb) – int64int32nag_int array
ldid = icomm(5)ldid=icomm 5 .
An array describing the parameter estimates returned in b. The first nlsv × nrfnlsv×nrf columns of id describe the parameter estimates for the random effects and the last nff columns the parameter estimates for the fixed effects.
A print function for decoding the parameter estimates given in b using information from id is supplied with the example program for this function.
For fixed effects:
  • for l = nrf × nlsv + 1 ,, nrf × nlsv + nff l=nrf×nlsv+1 ,, nrf×nlsv+nff
    • if b(l)bl contains the parameter estimate for the intercept then
      id(1,l) = id(2,l) = id(3,l) = 0 ;
      id1l = id2l = id3l = 0 ;
    • if b(l)bl contains the parameter estimate for the iith level of the jjth fixed variable, that is the vector of values held in the kkth column of dat when fixed(j + 2) = kfixedj+2=k then
      id(1,l) = 0, 
      id(2,l) = j, 
      id(3,l) = i;
      id1l=0,  id2l=j,  id3l=i;
    • if the jjth variable is continuous or binary, that is levels(fixed(j + 2)) = 1levelsfixedj+2=1, then id(3,l) = 0id3l=0;
    • any remaining rows of the llth column of id are set to 00.
For random effects:
  • let
    • NRbNRb denote the number of random variables in the bbth random statement, that is NRb = RNDM(1,b)NRb=RNDM(1,b);
    • RjbRjb denote the jjth random variable from the bbth random statement, that is the vector of values held in the kkth column of dat when RNDM(2 + j,b) = kRNDM(2+j,b)=k;
    • NSbNSb denote the number of subject variables in the bbth random statement, that is NSb = RNDM(3 + NRb,b)NSb=RNDM(3+NRb,b);
    • SjbSjb denote the jjth subject variable from the bbth random statement, that is the vector of values held in the kkth column of dat when RNDM(3 + NRb + j,b) = kRNDM(3+NRb+j,b)=k;
    • L(Sjb)L(Sjb) denote the number of levels for SjbSjb, that is L(Sjb) = levels(RNDM(3 + NRb + j,b))L(Sjb)=levelsRNDM(3+NRb+j,b);
  • then
    • for l = 1,2, nrf × nlsvl=1,2, nrf×nlsv, if b(l)bl contains the parameter estimate for the iith level of RjbRjb when Skb = skSkb=sk, for k = 1,2,,NSbk=1,2,,NSb and 1skL(Sjb)1skL(Sjb), i.e., sksk is a valid value for the kkth subject variable, then
      id(1,l) = b, 
      id(2,l) = j, 
      id(3,l) = i, 
      id(3 + k,l) = sk, ​k = 1,2,,NSb;
      id1l=b,  id2l=j,  id3l=i,  id3+kl=sk, ​k=1,2,,NSb;
    • if the parameter being estimated is for the intercept then id(2,l) = id(3,l) = 0id2l=id3l=0;
    • if the jjth variable is continuous, or binary, that is L(Sjb) = 1L(Sjb)=1, then id(3,l) = 0id3l=0;
    • the remaining rows of the llth column of id are set to 00.
In some situations, certain combinations of variables are never observed. In such circumstances all elements of the llth row of id are set to 999-999.
7:     b(lb) – double array
The parameter estimates, with the first nrf × nlsvnrf×nlsv elements of bb containing the parameter estimates for the random effects, νν, and the remaining nff elements containing the parameter estimates for the fixed effects, ββ. The order of these estimates are described by the id parameter.
8:     se(lb) – double array
The standard errors of the parameter estimates given in b.
9:     czz(ldczz, : :) – double array
The first dimension of the array czz will be nrfnrf
The second dimension of the array will be nrf × nlsvnrf×nlsv
ldczznrfldczznrf.
If nlsv = 1nlsv=1, then czz holds the lower triangular portion of the matrix (1 / σ2) (ZT1Z + 1) ( 1/ σ 2 ) ( ZT R^-1 Z + G^-1 ) , where R^ and G^ are the estimates of RR and GG respectively. If nlsv > 1nlsv>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:   cxx(ldcxx, : :) – double array
The first dimension of the array cxx will be nffnff
The second dimension of the array will be nffnff
ldcxxnffldcxxnff.
cxx holds the lower triangular portion of the matrix (1 / σ2) XT 1 X ( 1/ σ2 ) XT V^-1 X , where V^ is the estimated value of VV.
11:   cxz(ldcxz, : :) – double array
The first dimension of the array cxz will be nffnff
The second dimension of the array will be nlsv × nrfnlsv×nrf
ldcxznffldcxznff.
If nlsv = 1nlsv=1, then cxz holds the matrix (1 / σ2) (XT1Z) (1/σ2) ( XT V^-1 Z ) G^ , where V^ and G^ are the estimates of VV and GG respectively. If nlsv > 1nlsv>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:   ifail – int64int32nag_int scalar
ifail = 0ifail=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.

  ifail = 1ifail=1
lvpr is too small.
  ifail = 2ifail=2
Constraint: 1vpr(i)nvpr1vprinvpr.
  ifail = 3ifail=3
Constraint: .
  ifail = 4ifail=4
Constraint: gamma(i)0gammai0.
  ifail = 9ifail=9
lb is too small. lb is too small.
  ifail = 11ifail=11
ldid is too small.
  ifail = 15ifail=15
ldczz is too small.
  ifail = 17ifail=17
ldcxx is too small.
  ifail = 19ifail=19
ldcxz is too small.
  ifail = 21ifail=21
On entry, icomm has not been initialized correctly. On entry, icomm has not been initialized correctly. On entry, icomm has not been initialized correctly. On entry, icomm has not been initialized correctly. On entry, icomm has not been initialized correctly. On entry, icomm has not been initialized correctly.
  ifail = 32ifail=32
On entry, at least one value of i, for i = 1,2,,nvpri=1,2,,nvpr, does not appear in vpr.
W ifail = 101ifail=101
Optimal solution found, but requested accuracy not achieved.
  ifail = 102ifail=102
Too many major iterations.
W ifail = 103ifail=103
Current point cannot be improved upon.
W ifail = 104ifail=104
At least one negative estimate for gamma was obtained. All negative estimates have been set to zero.
  ifail = 999ifail=-999
Dynamic memory allocation failed.

Accuracy

Not applicable.

Further Comments

The parameter vpr gives the mapping between the random variables and the variance components. In most cases vpr(i) = ivpri=i, for i = 1,2,,iRNDM(1,i) + RNDM(2,i)i=1,2,,iRNDM(1,i)+RNDM(2,i). 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  
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)
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
RNDM =
(00)
fixed = ;
(30123)
T
.
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)
.
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  
DAT= 100103.6 010104.5 001101.1 100018.3 010017.2 001016.1
where each column only has one level
levels =
(111111)
.
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
RNDM =
(00)
fixed = ;
(60123456)
T
.
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)
.
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

function nag_correg_mixeff_hier_ml_example
% Problem size
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] = ...
  nag_correg_mixeff_hier_init(dat, levels, y, fixed, rndm, lrcomm, licomm);
licomm = icomm(1);
lrcomm = icomm(2);

% Pre-process the data
[nff, nlsv, nrf, rcomm, icomm, ifail] = ...
  nag_correg_mixeff_hier_init(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] = ...
    nag_correg_mixeff_hier_ml(vpr, nvpr, gamma, lb, rcomm, icomm, iopt, ropt);

% Print results
fprintf('\nNumber 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');

tdid = double(nff + nrf*nlsv);

if nrf > 0
  fprintf('\nRandom Effects\n');
end
pb = -999;
pfmt = ' ';
for k = 1:double(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+double(ns),k));
    if ( (pb ~= tb) || (strcmp(tfmt, pfmt) == 0) )
      if (k ~= 1)
        fprintf('\n');
      end
      fprintf('  Subject: ');
      for l=1:double(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                  %10.4f %10.4f\n', b(k), se(k));
    else
      % vid'th variable specified in rndm
      aid = rndm(2+vid,tb);
      if (id(3,k)==0)
        fprintf('    Variable %2d                %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 = double(nrf*nlsv + 1):tdid
  if vid~=-999
    vid = id(2,k);
    if vid==0
      % Intercept
      fprintf('    Intercept                  %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                %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:double(nvpr)
  fprintf('%10.5f     ', gamma(k));
  p = 0;
  for tb = 1:double(nrndm)
    nv = double(rndm(1,tb));
    ns = double(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);
 

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

function g02je_example
% Problem size
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('\nNumber 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');

tdid = double(nff + nrf*nlsv);

if nrf > 0
  fprintf('\nRandom Effects\n');
end
pb = -999;
pfmt = ' ';
for k = 1:double(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+double(ns),k));
    if ( (pb ~= tb) || (strcmp(tfmt, pfmt) == 0) )
      if (k ~= 1)
        fprintf('\n');
      end
      fprintf('  Subject: ');
      for l=1:double(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                  %10.4f %10.4f\n', b(k), se(k));
    else
      % vid'th variable specified in rndm
      aid = rndm(2+vid,tb);
      if (id(3,k)==0)
        fprintf('    Variable %2d                %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 = double(nrf*nlsv + 1):tdid
  if vid~=-999
    vid = id(2,k);
    if vid==0
      % Intercept
      fprintf('    Intercept                  %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                %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:double(nvpr)
  fprintf('%10.5f     ', gamma(k));
  p = 0;
  for tb = 1:double(nrndm)
    nv = double(rndm(1,tb));
    ns = double(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);
 

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


PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013