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_init (g02jc)

## Purpose

nag_correg_mixeff_hier_init (g02jc) preprocesses a dataset prior to fitting a linear mixed effects regression model of the following form via either nag_correg_mixeff_hier_reml (g02jd) or nag_correg_mixeff_hier_ml (g02je).

## Syntax

[nff, nlsv, nrf, rcomm, icomm, ifail] = g02jc(dat, levels, y, fixed, rndm, lrcomm, licomm, 'n', n, 'ncol', ncol, 'wt', wt, 'lfixed', lfixed, 'nrndm', nrndm)
[nff, nlsv, nrf, rcomm, icomm, ifail] = nag_correg_mixeff_hier_init(dat, levels, y, fixed, rndm, lrcomm, licomm, 'n', n, 'ncol', ncol, 'wt', wt, 'lfixed', lfixed, 'nrndm', nrndm)
Note: the interface to this routine has changed since earlier releases of the toolbox:
Mark 23: weight omitted from interface
.

## Description

nag_correg_mixeff_hier_init (g02jc) must be called prior to fitting a linear mixed effects regression model with either nag_correg_mixeff_hier_reml (g02jd) or nag_correg_mixeff_hier_ml (g02je).
The model fitting functions nag_correg_mixeff_hier_reml (g02jd) and nag_correg_mixeff_hier_ml (g02je) fit a model of the following form:
 y = Xβ + Zν + ε $y=Xβ+Zν+ε$
 where y$y$ is a vector of n$n$ observations on the dependent variable, X$X$ is an n$n$ by p$p$ design matrix of fixed independent variables, β$\beta$ is a vector of p$p$ unknown fixed effects, Z$Z$ is an n$n$ by q$q$ design matrix of random independent variables, ν$\nu$ is a vector of length q$q$ of unknown random effects, ε$\epsilon$ is a vector of length n$n$ of unknown random errors,
and ν$\nu$ and ε$\epsilon$ are Normally distributed with expectation zero and variance/covariance matrix defined by
Var
 [ ν ε ]
=
 [ G 0 0 R ]
$Var[ ν ε ] = [ G 0 0 R ]$
where R = σR2 I $R={\sigma }_{R}^{2}I$, I$I$ is the n × n$n×n$ identity matrix and G$G$ is a diagonal matrix.
Case weights can be incorporated into the model by replacing X$X$ and Z$Z$ with Wc1 / 2X${W}_{c}^{1/2}X$ and Wc1 / 2Z${W}_{c}^{1/2}Z$ respectively where Wc${W}_{c}$ is a diagonal weight matrix.

None.

## Parameters

### Compulsory Input Parameters

1:     dat(lddat,ncol) – double array
lddat, the first dimension of the array, must satisfy the constraint lddatn$\mathit{lddat}\ge {\mathbf{n}}$.
A matrix of data, with dat(i,j)${\mathbf{dat}}\left(i,j\right)$ holding the i$i$th observation on the j$j$th variable. The two design matrices X$X$ and Z$Z$ are constructed from dat and the information given in fixed (for X$X$) and rndm (for Z$Z$).
Constraint: if levels(j)1,1dat(i,j)levels(j)${\mathbf{levels}}\left(j\right)\ne 1,1\le {\mathbf{dat}}\left(i,j\right)\le {\mathbf{levels}}\left(j\right)$.
2:     levels(ncol) – int64int32nag_int array
ncol, the dimension of the array, must satisfy the constraint ncol0${\mathbf{ncol}}\ge 0$.
levels(i)${\mathbf{levels}}\left(i\right)$ contains the number of levels associated with the i$i$th variable held in dat.
If the i$i$th variable is continuous or binary (i.e., only takes the values zero or one) then levels(i)${\mathbf{levels}}\left(i\right)$ must be set to 1$1$. Otherwise the i$i$th variable is assumed to take an integer value between 1$1$ and levels(i)${\mathbf{levels}}\left(i\right)$, (i.e., the i$i$th variable is discrete with levels(i)${\mathbf{levels}}\left(i\right)$ levels).
Constraint: levels(i)1${\mathbf{levels}}\left(\mathit{i}\right)\ge 1$, for i = 1,2,,ncol$\mathit{i}=1,2,\dots ,{\mathbf{ncol}}$.
3:     y(n) – double array
n, the dimension of the array, must satisfy the constraint n1${\mathbf{n}}\ge 1$.
y$y$, the vector of observations on the dependent variable.
4:     fixed(lfixed) – int64int32nag_int array
lfixed, the dimension of the array, must satisfy the constraint lfixed2 + fixed(1)${\mathbf{lfixed}}\ge 2+{\mathbf{fixed}}\left(1\right)$.
Defines the structure of the fixed effects design matrix, X$X$.
fixed(1)${\mathbf{fixed}}\left(1\right)$
The number of variables, NF${N}_{F}$, to include as fixed effects (not including the intercept if present).
fixed(2)${\mathbf{fixed}}\left(2\right)$
The fixed intercept flag which must contain 1$1$ if a fixed intercept is to be included and 0$0$ otherwise.
fixed(2 + i)${\mathbf{fixed}}\left(2+i\right)$
The column of dat holding the i$\mathit{i}$th fixed variable, for i = 1,2,,fixed(1)$\mathit{i}=1,2,\dots ,{\mathbf{fixed}}\left(1\right)$.
See Section [Construction of the design matrix, ] for more details on the construction of X$X$.
Constraints:
• fixed(1)0${\mathbf{fixed}}\left(1\right)\ge 0$;
• fixed(2) = 0​ or ​1${\mathbf{fixed}}\left(2\right)=0\text{​ or ​}1$;
• 1fixed(2 + i)ncol$1\le {\mathbf{fixed}}\left(2+\mathit{i}\right)\le {\mathbf{ncol}}$, for i = 1,2,,fixed(1)$\mathit{i}=1,2,\dots ,{\mathbf{fixed}}\left(1\right)$.
5:     rndm(ldrndm,nrndm) – int64int32nag_int array
ldrndm, the first dimension of the array, must satisfy the constraint ldrndm max b  (3 + NRb + NSb) $\mathit{ldrndm}\ge \underset{b}{max}\phantom{\rule{0.25em}{0ex}}\left(3+{N}_{{R}_{b}}+{N}_{{S}_{b}}\right)$.
rndm(i,j)${\mathbf{rndm}}\left(i,j\right)$ defines the structure of the random effects design matrix, Z$Z$. The b$b$th column of rndm defines a block of columns in the design matrix Z$Z$.
rndm(1,b)${\mathbf{rndm}}\left(1,b\right)$
The number of variables, NRb${N}_{{R}_{b}}$, to include as random effects in the b$b$th block (not including the random intercept if present).
rndm(2,b)${\mathbf{rndm}}\left(2,b\right)$
The random intercept flag which must contain 1$1$ if block b$b$ includes a random intercept and 0$0$ otherwise.
rndm(2 + i,b)${\mathbf{rndm}}\left(2+i,b\right)$
The column of dat holding the i$\mathit{i}$th random variable in the b$b$th block, for i = 1,2,,rndm(1,b)$\mathit{i}=1,2,\dots ,{\mathbf{rndm}}\left(1,b\right)$.
rndm(3 + NRb,b)${\mathbf{rndm}}\left(3+{N}_{{R}_{b}},b\right)$
The number of subject variables, NSb${N}_{{S}_{b}}$, for the b$b$th block. The subject variables define the nesting structure for this block.
rndm(3 + NRb + i,b)${\mathbf{rndm}}\left(3+{N}_{{R}_{b}}+i,b\right)$
The column of dat holding the i$\mathit{i}$th subject variable in the b$b$th block, for i = 1,2,,rndm(3 + NRb,b)$\mathit{i}=1,2,\dots ,{\mathbf{rndm}}\left(3+{N}_{{R}_{b}},b\right)$.
See Section [Construction of design matrix, ] for more details on the construction of Z$Z$.
Constraints:
• rndm(1,b)0${\mathbf{rndm}}\left(1,b\right)\ge 0$;
• rndm(2,b) = 0​ or ​1${\mathbf{rndm}}\left(2,b\right)=0\text{​ or ​}1$;
• at least one random variable or random intercept must be specified in each block, i.e., rndm(1,b) + rndm(2,b) > 0 ${\mathbf{rndm}}\left(1,b\right)+{\mathbf{rndm}}\left(2,b\right)>0$;
• the column identifiers associated with the random variables must be in the range 1$1$ to ncol, i.e., 1 rndm(2 + i,b) ncol $1\le {\mathbf{rndm}}\left(2+\mathit{i},b\right)\le {\mathbf{ncol}}$, for i = 1,2,,rndm(1,b)$\mathit{i}=1,2,\dots ,{\mathbf{rndm}}\left(1,b\right)$;
• rndm(3 + NRb,b) 0 ${\mathbf{rndm}}\left(3+{N}_{{R}_{b}},b\right)\ge 0$;
• the column identifiers associated with the subject variables must be in the range 1$1$ to ncol, i.e., 1 rndm(3 + NRb + i ,b) ncol $1\le {\mathbf{rndm}}\left(3+{N}_{{R}_{b}}+\mathit{i},b\right)\le {\mathbf{ncol}}$, for i = 1,2,,rndm(3 + NRb,b)$\mathit{i}=1,2,\dots ,{\mathbf{rndm}}\left(3+{N}_{{R}_{b}},b\right)$.
6:     lrcomm – int64int32nag_int scalar
The dimension of the array rcomm as declared in the (sub)program from which nag_correg_mixeff_hier_init (g02jc) is called.
Constraint: lrcommnrf × nlsv + nff + nff × nlsv + nrf × nlsv + nff + 2${\mathbf{lrcomm}}\ge {\mathbf{nrf}}×{\mathbf{nlsv}}+{\mathbf{nff}}+{\mathbf{nff}}×{\mathbf{nlsv}}+{\mathbf{nrf}}×{\mathbf{nlsv}}+{\mathbf{nff}}+2$.
7:     licomm – int64int32nag_int scalar
The dimension of the array icomm as declared in the (sub)program from which nag_correg_mixeff_hier_init (g02jc) is called.
Constraint: licomm = 2${\mathbf{licomm}}=2$ or
licomm34 + NF × (MFL + 1) + nrndm × MNR × MRL + (LRNDM + 2) × nrndm +   ncol + LDID × LB, ${\mathbf{licomm}}\ge 34+{N}_{F}×\left(\text{MFL}+1\right)+{\mathbf{nrndm}}×\text{MNR}×\text{MRL}+\left(\text{LRNDM}+2\right)×{\mathbf{nrndm}}+\phantom{\rule{0ex}{0ex}}{\mathbf{ncol}}+\text{LDID}×\text{LB,}$.
where
• MNR = maxb  (NRb) $\text{MNR}=\underset{b}{\mathrm{max}}\phantom{\rule{0.25em}{0ex}}\left({N}_{{R}_{b}}\right)$,
• MFL = maxi  (levels(fixed(2 + i))) $\text{MFL}=\underset{i}{\mathrm{max}}\phantom{\rule{0.25em}{0ex}}\left({\mathbf{levels}}\left({\mathbf{fixed}}\left(2+i\right)\right)\right)$,
• MRL = maxb,i  (levels(rndm(2 + i,b))) $\text{MRL}=\underset{b,i}{\mathrm{max}}\phantom{\rule{0.25em}{0ex}}\left({\mathbf{levels}}\left({\mathbf{rndm}}\left(2+i,b\right)\right)\right)$,
• LDID = maxb  NSb $\text{LDID}=\underset{b}{\mathrm{max}}\phantom{\rule{0.25em}{0ex}}{N}_{{S}_{b}}$,
• LB = nff + nrf × nlsv$\text{LB}={\mathbf{nff}}+{\mathbf{nrf}}×{\mathbf{nlsv}}$, and
• LRNDM = max b  (3 + NRb + NSb) $\text{LRNDM}=\underset{b}{max}\phantom{\rule{0.25em}{0ex}}\left(3+{N}_{{R}_{b}}+{N}_{{S}_{b}}\right)$

### Optional Input Parameters

1:     n – int64int32nag_int scalar
Default: The dimension of the array y and the first dimension of the array dat. (An error is raised if these dimensions are not equal.)
n$n$, the number of observations.
The effective number of observations, that is the number of observations with nonzero weight (see wt for more detail), must be greater than the number of fixed effects in the model (as returned in nff).
Constraint: n1${\mathbf{n}}\ge 1$.
2:     ncol – int64int32nag_int scalar
Default: The dimension of the arrays dat, levels and the second dimension of the array dat. (An error is raised if these dimensions are not equal.)
The number of columns in the data matrix, dat.
Constraint: ncol0${\mathbf{ncol}}\ge 0$.
3:     wt( : $:$) – double array
Note: the dimension of the array wt must be at least n${\mathbf{n}}$ if weight = 'W'$\mathit{weight}=\text{'W'}$.
If weight = 'W'$\mathit{weight}=\text{'W'}$, wt must contain the diagonal elements of the weight matrix Wc${W}_{c}$.
If wt(i) = 0.0${\mathbf{wt}}\left(i\right)=0.0$, the i$i$th observation is not included in the model and the effective number of observations is the number of observations with nonzero weights.
If weight = 'U'$\mathit{weight}=\text{'U'}$, wt is not referenced and the effective number of observations is n$n$.
Constraint: if weight = 'W'$\mathit{weight}=\text{'W'}$, wt(i)0.0${\mathbf{wt}}\left(\mathit{i}\right)\ge 0.0$, for i = 1,2,,n$\mathit{i}=1,2,\dots ,n$.
4:     lfixed – int64int32nag_int scalar
Default: The dimension of the array fixed.
Length of the vector fixed.
Constraint: lfixed2 + fixed(1)${\mathbf{lfixed}}\ge 2+{\mathbf{fixed}}\left(1\right)$.
5:     nrndm – int64int32nag_int scalar
Default: The second dimension of the array rndm.
The second dimension of the array rndm as declared in the (sub)program from which nag_correg_mixeff_hier_init (g02jc) is called.
Constraint: nrndm > 0${\mathbf{nrndm}}>0$.

### Input Parameters Omitted from the MATLAB Interface

weight lddat ldrndm

### Output Parameters

1:     nff – int64int32nag_int scalar
p$p$, the number of fixed effects estimated, i.e., the number of columns in the design matrix X$X$.
2:     nlsv – int64int32nag_int scalar
The number of levels for the overall subject variable (see Section [Construction of design matrix, ] for a description of what this means). If there is no overall subject variable, nlsv = 1${\mathbf{nlsv}}=1$.
3:     nrf – int64int32nag_int scalar
The number of random effects estimated in each of the overall subject blocks. The number of columns in the design matrix Z$Z$ is given by q = nrf × nlsv$q={\mathbf{nrf}}×{\mathbf{nlsv}}$.
4:     rcomm(lrcomm) – double array
Communication array as required by the analysis functions nag_correg_mixeff_hier_reml (g02jd) and nag_correg_mixeff_hier_ml (g02je).
5:     icomm(licomm) – int64int32nag_int array
If licomm = 2${\mathbf{licomm}}=2$, icomm(1)${\mathbf{icomm}}\left(1\right)$ holds the minimum required value for licomm and icomm(2)${\mathbf{icomm}}\left(2\right)$ holds the minimum required value for lrcomm, otherwise icomm is a communication array as required by the analysis functions nag_correg_mixeff_hier_reml (g02jd) and nag_correg_mixeff_hier_ml (g02je).
6:     ifail – int64int32nag_int scalar
${\mathrm{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:
ifail = 1${\mathbf{ifail}}=1$
On entry, weight had an illegal value.
ifail = 2${\mathbf{ifail}}=2$
Constraint: n1${\mathbf{n}}\ge 1$.
ifail = 3${\mathbf{ifail}}=3$
Constraint: ncol0${\mathbf{ncol}}\ge 0$.
ifail = 4${\mathbf{ifail}}=4$
On entry, variable j$j$ of observation i$i$ is less than 1$1$ or greater than levels(j)${\mathbf{levels}}\left(j\right)$.
ifail = 5${\mathbf{ifail}}=5$
Constraint: lddatn$\mathit{lddat}\ge {\mathbf{n}}$.
ifail = 6${\mathbf{ifail}}=6$
Constraint: levels(i)1${\mathbf{levels}}\left(i\right)\ge 1$.
ifail = 8${\mathbf{ifail}}=8$
Constraint: wt(i)0.0${\mathbf{wt}}\left(i\right)\ge 0.0$.
ifail = 9${\mathbf{ifail}}=9$
On entry, number of fixed parameters, _$_$ is less than zero.
ifail = 10${\mathbf{ifail}}=10$
lfixed is too small.
ifail = 11${\mathbf{ifail}}=11$
Constraint: nrndm > 0${\mathbf{nrndm}}>0$.
ifail = 12${\mathbf{ifail}}=12$
On entry, number of random parameters for random statement i$i$ is less than 0$0$.
ifail = 13${\mathbf{ifail}}=13$
ldrndm is too small.
ifail = 18${\mathbf{ifail}}=18$
lrcomm is too small.
ifail = 20${\mathbf{ifail}}=20$
licomm is too small.
ifail = 102${\mathbf{ifail}}=102$
n is too small.
ifail = 108${\mathbf{ifail}}=108$
On entry, no observations due to zero weights.
ifail = 109${\mathbf{ifail}}=109$
On entry, invalid value for fixed intercept flag.
ifail = 112${\mathbf{ifail}}=112$
On entry, invalid value for random intercept flag for random statement i$i$.
ifail = 209${\mathbf{ifail}}=209$
On entry, index of fixed variable j$j$ is less than 1$1$ or greater than ncol${\mathbf{ncol}}$.
ifail = 212${\mathbf{ifail}}=212$
On entry, must be at least one parameter, or an intercept in each random statement i$i$:
ifail = 312${\mathbf{ifail}}=312$
On entry, index of random variable j$j$ in random statement i$i$ is less than 1$1$ or greater than ncol${\mathbf{ncol}}$.
ifail = 412${\mathbf{ifail}}=412$
On entry, number of subject parameters for random statement i$i$ is less than 0$0$.
ifail = 512${\mathbf{ifail}}=512$
On entry, nesting variable j$j$ in random statement i$i$ has one level.
ifail = 999${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.

## Accuracy

Not applicable.

### Construction of the fixed effects design matrix, X

Let
• NF${N}_{F}$ denote the number of fixed variables, that is fixed(1) = NF${\mathbf{fixed}}\left(1\right)={N}_{F}$;
• Fj${F}_{j}$ denote the j$j$th fixed variable, that is the vector of values held in the k$k$th column of dat when fixed(2 + j) = k${\mathbf{fixed}}\left(2+j\right)=k$;
• Fij${F}_{ij}$ denote the i$i$th element of Fj${F}_{j}$;
• L(Fj)$L\left({F}_{j}\right)$ denote the number of levels for Fj${F}_{j}$, that is L(Fj) = levels(fixed(2 + j))$L\left({F}_{j}\right)={\mathbf{levels}}\left({\mathbf{fixed}}\left(2+j\right)\right)$;
• Dv(Fj)${D}_{v}\left({F}_{j}\right)$ denoted an indicator function that returns a vector of values whose i$i$th element is 1$1$ if Fij = v${F}_{ij}=v$ and 0$0$ otherwise.
The design matrix for the fixed effects, X$X$, is constructed as follows:
• set k$k$ to one and the flag done_first$\text{done_first}$ to false;
• if a fixed intercept is included, that is fixed(2) = 1${\mathbf{fixed}}\left(2\right)=1$,
• set the first column of X$X$ to a vector of 1$1$s;
• set k = k + 1$k=k+1$;
• set done_first$\text{done_first}$ to true;
• loop over each fixed variable, so for each j = 1,2,,NF$j=1,2,\dots ,{N}_{F}$,
• if L(Fj) = 1$L\left({F}_{j}\right)=1$,
• set the k$k$th column of X$X$ to be Fj${F}_{j}$;
• set k = k + 1$k=k+1$;
• else
• if done_first$\text{done_first}$ is false then
• set the L(Fj)$L\left({F}_{j}\right)$ columns, k$k$ to k + L(Fj)1$k+L\left({F}_{j}\right)-1$, of X$X$ to Dv(Fj)${D}_{\mathit{v}}\left({F}_{j}\right)$, for v = 1,2,,L(Fj)$\mathit{v}=1,2,\dots ,L\left({F}_{j}\right)$;
• set k = k + L(Fj)$k=k+L\left({F}_{j}\right)$;
• set done_first$\text{done_first}$ to true;
• else
• set the L(Fj)1$L\left({F}_{j}\right)-1$ columns, k$k$ to k + L(Fj)2$k+L\left({F}_{j}\right)-2$, of X$X$ to Dv(Fj)${D}_{\mathit{v}}\left({F}_{j}\right)$, for v = 2,3,,L(Fj)$\mathit{v}=2,3,\dots ,L\left({F}_{j}\right)$;
• set k = k + L(Fj)1$k=k+L\left({F}_{j}\right)-1$.
The number of columns in the design matrix, X$X$, is therefore given by
 p = 1 + ∑ j = 1 NF (levels(fixed( 2 + j )) − 1) . $p= 1+ ∑ j=1 N F ( levels fixed 2+j -1 ) .$
This quantity is returned in nff.
In summary, nag_correg_mixeff_hier_init (g02jc) converts all non-binary categorical variables (i.e., where L(Fj) > 1$L\left({F}_{j}\right)>1$) to dummy variables. If a fixed intercept is included in the model then the first level of all such variables is dropped. If a fixed intercept is not included in the model then the first level of all such variables, other than the first, is dropped. The variables are added into the model in the order they are specified in fixed.

### Construction of random effects design matrix, Z

Let
• NRb${N}_{{R}_{b}}$ denote the number of random variables in the b$b$th random statement, that is NRb = rndm(1,b)${N}_{{R}_{b}}={\mathbf{rndm}}\left(1,b\right)$;
• Rjb${R}_{jb}$ denote the j$j$th random variable from the b$b$th random statement, that is the vector of values held in the k$k$th column of dat when rndm(2 + j,b) = k${\mathbf{rndm}}\left(2+j,b\right)=k$;
• Rijb${R}_{ijb}$ denote the i$i$th element of Rjb${R}_{jb}$;
• L(Rjb)$L\left({R}_{jb}\right)$ denote the number of levels for Rjb${R}_{jb}$, that is L(Rjb) = levels(rndm(2 + j,b))$L\left({R}_{jb}\right)={\mathbf{levels}}\left({\mathbf{rndm}}\left(2+j,b\right)\right)$;
• Dv(Rjb)${D}_{v}\left({R}_{jb}\right)$ denoted an indicator function that returns a vector of values whose i$i$th element is 1$1$ if Rijb = v${R}_{ijb}=v$ and 0$0$ otherwise;
• NSb${N}_{{S}_{b}}$ denote the number of subject variables in the b$b$th random statement, that is NSb = rndm(3 + NRb,b)${N}_{{S}_{b}}={\mathbf{rndm}}\left(3+{N}_{{R}_{b}},b\right)$;
• Sjb${S}_{jb}$ denote the j$j$th subject variable from the b$b$th random statement, that is the vector of values held in the k$k$th column of dat when rndm(3 + NRb + j,b) = k${\mathbf{rndm}}\left(3+{N}_{{R}_{b}}+j,b\right)=k$;
• Sijb${S}_{ijb}$ denote the i$i$th element of Sjb${S}_{jb}$;
• L(Sjb)$L\left({S}_{jb}\right)$ denote the number of levels for Sjb${S}_{jb}$, that is L(Sjb) = levels(rndm(3 + NRb + j,b))$L\left({S}_{jb}\right)={\mathbf{levels}}\left({\mathbf{rndm}}\left(3+{N}_{{R}_{b}}+j,b\right)\right)$;
• Ib (s1,s2,,sNSb) ${I}_{b}\left({s}_{1},{s}_{2},\dots ,{s}_{{N}_{{S}_{b}}}\right)$ denoted an indicator function that returns a vector of values whose i$i$th element is 1$1$ if Sijb = sj ${S}_{ijb}={s}_{j}$ for all j = 1,2,,NSb$j=1,2,\dots ,{N}_{{S}_{b}}$ and 0$0$ otherwise.
The design matrix for the random effects, Z$Z$, is constructed as follows:
• set k$k$ to one;
• loop over each random statement, so for each b = 1,2,,nrndm$b=1,2,\dots ,{\mathbf{nrndm}}$,
• loop over each level of the last subject variable, so for each sNSb = 1,2,,L(RNSb b)${s}_{{N}_{{S}_{b}}}=1,2,\dots ,L\left({R}_{{N}_{{S}_{b}}b}\right)$,
• $⋮$
• loop over each level of the second subject variable, so for each s2 = 1,2,,L(R2b)${s}_{2}=1,2,\dots ,L\left({R}_{2b}\right)$,
• loop over each level of the first subject variable, so for each s1 = 1,2,,L(R1b)${s}_{1}=1,2,\dots ,L\left({R}_{1b}\right)$,
• if a random intercept is included, that is rndm(2,b) = 1${\mathbf{rndm}}\left(2,b\right)=1$,
• set the k$k$th column of Z$Z$ to Ib(s1,s2,,sNSb)${I}_{b}\left({s}_{1},{s}_{2},\dots ,{s}_{{N}_{{S}_{b}}}\right)$;
• set k = k + 1$k=k+1$;
• loop over each random variable in the b$b$th random statement, so for each j = 1,2,,NRb$j=1,2,\dots ,{N}_{{R}_{b}}$,
• if L(Rjb) = 1$L\left({R}_{jb}\right)=1$,
• set the k$k$th column of Z$Z$ to Rjb × Ib (s1,s2,,sNSb) ${R}_{jb}×{I}_{b}\left({s}_{1},{s}_{2},\dots ,{s}_{{N}_{{S}_{b}}}\right)$ where × $×$ indicates an element-wise multiplication between the two vectors, Rjb${R}_{jb}$ and Ib()${I}_{b}\left(\dots \right)$;
• set k = k + 1$k=k+1$;
• else
• set the L(Rbj)$L\left({R}_{bj}\right)$ columns, k$k$ to k + L(Rbj)$k+L\left({R}_{bj}\right)$, of Z$Z$ to Dv(Rjb) × Ib(s1,s2,,sNSb)${D}_{\mathit{v}}\left({R}_{jb}\right)×{I}_{b}\left({s}_{1},{s}_{2},\dots ,{s}_{{N}_{{S}_{b}}}\right)$, for v = 1,2,,L(Rjb)$\mathit{v}=1,2,\dots ,L\left({R}_{jb}\right)$. As before, × $×$ indicates an element-wise multiplication between the two vectors, Dv()${D}_{v}\left(\dots \right)$ and Ib()${I}_{b}\left(\dots \right)$;
• set k = k + L(Rjb)$k=k+L\left({R}_{jb}\right)$.
In summary, each column of rndm defines a block of consecutive columns in Z$Z$. nag_correg_mixeff_hier_init (g02jc) converts all non-binary categorical variables (i.e., where L(Rjb)$L\left({R}_{jb}\right)$ or L(Sjb) > 1$L\left({S}_{jb}\right)>1$) to dummy variables. All random variables defined within a column of rndm are nested within all subject variables defined in the same column of rndm. In addition each of the subject variables are nested within each other, starting with the first (i.e., each of the Rjb,j = 1,2,,NRb${R}_{jb},j=1,2,\dots ,{N}_{{R}_{b}}$ are nested within S1b${S}_{1b}$ which in turn is nested within S2b${S}_{2b}$, which in turn is nested within S3b${S}_{3b}$, etc.).
If the last subject variable in each column of rndm are the same (i.e., SNS11 = SNS22 = = SNSbb ${S}_{{N}_{{S}_{1}}1}={S}_{{N}_{{S}_{2}}2}=\dots ={S}_{{N}_{{S}_{b}}b}$) then all random effects in the model are nested within this variable. In such instances the last subject variable ( SNS11 ${S}_{{N}_{{S}_{1}}1}$) is called the overall subject variable. The fact that all of the random effects in the model are nested within the overall subject variable means that ZTZ${Z}^{\mathrm{T}}Z$ is block diagonal in structure. This fact can be utilised to improve the efficiency of the underlying computation and reduce the amount of internal storage required. The number of levels in the overall subject variable is returned in nlsv = L(SNS11)${\mathbf{nlsv}}=L\left({S}_{{N}_{{S}_{1}}1}\right)$.
If the last k$k$ subject variables in each column of rndm are the same, for k > 1$k>1$ then the overall subject variable is defined as the interaction of these k$k$ variables and
 NS1 nlsv = ∏ L(Sj1). j = NS1 − k + 1
$nlsv= ∏ j=NS1-k+1 NS1 L(Sj1) .$
If there is no overall subject variable then nlsv = 1${\mathbf{nlsv}}=1$.
The number of columns in the design matrix Z$Z$ is given by q = nrf × nlsv$q={\mathbf{nrf}}×{\mathbf{nlsv}}$.

### The rndm parameter

To illustrate some additional points about the rndm parameter, we assume that we have a dataset with three discrete variables, V1${V}_{1}$, V2${V}_{2}$ and V3${V}_{3}$, with 2,4$2,4$ and 3$3$ levels respectively, and that V1${V}_{1}$ is in the first column of dat, V2${V}_{2}$ in the second and V3${V}_{3}$ the third. Also assume that we wish to fit a model containing V1${V}_{1}$ along with V2${V}_{2}$ nested within V3${V}_{3}$, as random effects. In order to do this the rndm matrix requires two columns:
RNDM =
 1 1 0 0 1 2 0 1 0 3
$RNDM= 1 1 0 0 1 2 0 1 0 3$
The first column, (1,0,1,0,0)$\left(1,0,1,0,0\right)$, indicates one random variable (rndm(1,1) = 1${\mathbf{rndm}}\left(1,1\right)=1$), no intercept (rndm(2,1) = 0${\mathbf{rndm}}\left(2,1\right)=0$), the random variable is in the first column of dat (rndm(3,1) = 1${\mathbf{rndm}}\left(3,1\right)=1$), there are no subject variables; as no nesting is required for V1${V}_{1}$ (rndm(4,1) = 0${\mathbf{rndm}}\left(4,1\right)=0$). The last element in this column is ignored.
The second column, (1,0,2,1,3)$\left(1,0,2,1,3\right)$, indicates one random variable (rndm(1,2) = 1${\mathbf{rndm}}\left(1,2\right)=1$), no intercept (rndm(2,2) = 0${\mathbf{rndm}}\left(2,2\right)=0$), the random variable is in the second column of dat (rndm(3,2) = 2)$\left({\mathbf{rndm}}\left(3,2\right)=2\right)$, there is one subject variable (rndm(4,2) = 1${\mathbf{rndm}}\left(4,2\right)=1$), and the subject variable is in the third column of dat (rndm(5,2) = 3)$\left({\mathbf{rndm}}\left(5,2\right)=3\right)$.
The corresponding Z$Z$ matrix would have 14$14$ columns, with 2$2$ coming from V1${V}_{1}$ and 12$12$ (4 × 3$4×3$) from V2${V}_{2}$ nested within V3${V}_{3}$. The, symmetric, ZTZ${Z}^{\mathrm{T}}Z$ matrix has the form
 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 0 0 0 0 0 0 0 0 - - - - - - 0 0 0 0 0 0 0 0 - - - - - - 0 0 0 0 0 0 0 0 - - - - - - 0 0 0 0 0 0 0 0 - - 0 0 0 0 - - - - 0 0 0 0 - - 0 0 0 0 - - - - 0 0 0 0 - - 0 0 0 0 - - - - 0 0 0 0 - - 0 0 0 0 - - - - 0 0 0 0 - - 0 0 0 0 0 0 0 0 - - - - - - 0 0 0 0 0 0 0 0 - - - - - - 0 0 0 0 0 0 0 0 - - - - - - 0 0 0 0 0 0 0 0 - - - -
$- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 0 0 0 0 0 0 0 0 - - - - - - 0 0 0 0 0 0 0 0 - - - - - - 0 0 0 0 0 0 0 0 - - - - - - 0 0 0 0 0 0 0 0 - - 0 0 0 0 - - - - 0 0 0 0 - - 0 0 0 0 - - - - 0 0 0 0 - - 0 0 0 0 - - - - 0 0 0 0 - - 0 0 0 0 - - - - 0 0 0 0 - - 0 0 0 0 0 0 0 0 - - - - - - 0 0 0 0 0 0 0 0 - - - - - - 0 0 0 0 0 0 0 0 - - - - - - 0 0 0 0 0 0 0 0 - - - -$
where 0$0$ indicates a structural zero, i.e., it always takes the value 0$0$, irrespective of the data, and -$-$ a value that is not a structural zero. The first two rows and columns of ZTZ${Z}^{\mathrm{T}}Z$ correspond to V1${V}_{1}$. The block diagonal matrix in the 12 rows and columns in the bottom right correspond to V2${V}_{2}$ nested within V3${V}_{3}$. With the 4 × 4$4×4$ blocks corresponding to the levels of V2${V}_{2}$. There are three blocks as the subject variable (V3${V}_{3}$) has three levels.
The model fitting functions, nag_correg_mixeff_hier_reml (g02jd) and nag_correg_mixeff_hier_ml (g02je), use the sweep algorithm to calculate the log likelihood function for a given set of variance components. This algorithm consists of moving down the diagonal elements (called pivots) of a matrix which is similar in structure to ZTZ${Z}^{\mathrm{T}}Z$, and updating each element in that matrix. When using the k$k$ diagonal element of a matrix A$A$, an element a i j ,ik,jk ${a}_{ij},i\ne k,j\ne k$, is adjusted by an amount equal to a i k a i j / a k k ${a}_{ik}{a}_{ij}/{a}_{kk}$. This process can be referred to as sweeping on the k$k$th pivot. As there are no structural zeros in the first row or column of the above ZTZ${Z}^{\mathrm{T}}Z$, sweeping on the first pivot of ZTZ${Z}^{\mathrm{T}}Z$ would alter each element of the matrix and therefore destroy the structural zeros, i.e., we could no longer guarantee they would be zero.
Reordering the rndm matrix to
RNDM =
 1 1 0 0 2 1 1 0 3 0
$RNDM= 1 1 0 0 2 1 1 0 3 0$
i.e., the swapping the two columns, results in a ZTZ${Z}^{\mathrm{T}}Z$ matrix of the form
 - - - - 0 0 0 0 0 0 0 0 - - - - - - 0 0 0 0 0 0 0 0 - - - - - - 0 0 0 0 0 0 0 0 - - - - - - 0 0 0 0 0 0 0 0 - - 0 0 0 0 - - - - 0 0 0 0 - - 0 0 0 0 - - - - 0 0 0 0 - - 0 0 0 0 - - - - 0 0 0 0 - - 0 0 0 0 - - - - 0 0 0 0 - - 0 0 0 0 0 0 0 0 - - - - - - 0 0 0 0 0 0 0 0 - - - - - - 0 0 0 0 0 0 0 0 - - - - - - 0 0 0 0 0 0 0 0 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
$- - - - 0 0 0 0 0 0 0 0 - - - - - - 0 0 0 0 0 0 0 0 - - - - - - 0 0 0 0 0 0 0 0 - - - - - - 0 0 0 0 0 0 0 0 - - 0 0 0 0 - - - - 0 0 0 0 - - 0 0 0 0 - - - - 0 0 0 0 - - 0 0 0 0 - - - - 0 0 0 0 - - 0 0 0 0 - - - - 0 0 0 0 - - 0 0 0 0 0 0 0 0 - - - - - - 0 0 0 0 0 0 0 0 - - - - - - 0 0 0 0 0 0 0 0 - - - - - - 0 0 0 0 0 0 0 0 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -$
This matrix is identical to the previous one, except the first two rows and columns have become the last two rows and columns. Sweeping a matrix, A = {aij}$A=\left\{{a}_{ij}\right\}$, of this form on the first pivot will only affect those elements aij${a}_{ij}$, where ai10​ and ​a1j0${a}_{i1}\ne 0\text{​ and ​}{a}_{1j}\ne 0$, which is only the 13$13$th and 14$14$th row and columns, and the top left hand block of 4$4$ rows and columns. The block diagonal nature of the first 12$12$ rows and columns therefore greatly reduces the amount of work the algorithm needs to perform.
nag_correg_mixeff_hier_init (g02jc) constructs the ZTZ${Z}^{\mathrm{T}}Z$ as specified by the rndm matrix, and does not attempt to reorder it to improve performance. Therefore for best performance some thought is required on what ordering to use. In general it is more efficient to structure rndm in such a way that the first row relates to the deepest level of nesting, the second to the next level, etc..

## Example

```function nag_correg_mixeff_hier_init_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_reml(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.1561     3.7946
Variable  3 (Level  2)         1.8951     3.9284
Variable  4 (Level  1)         0.6496     3.1617

Subject:  Variable 10 (Level  1) Variable 11 (Level  1) Variable 12 (Level  1)
Variable  4 (Level  3)         0.7390     3.1424

Subject:  Variable 10 (Level  2) Variable 11 (Level  1) Variable 12 (Level  1)
Variable  3 (Level  1)         1.4216     3.3773
Variable  3 (Level  2)        -2.8921     3.3953
Variable  4 (Level  1)         3.6789     2.3162
Variable  4 (Level  2)        -1.9742     2.3887
Variable  4 (Level  3)        -2.2088     2.0697

Subject:  Variable 10 (Level  1) Variable 11 (Level  2) Variable 12 (Level  1)
Variable  3 (Level  1)        -2.9659     3.9127
Variable  3 (Level  2)         2.7951     4.7183
Variable  4 (Level  1)        -4.7330     2.3094
Variable  4 (Level  2)         5.5161     2.2330
Variable  4 (Level  3)        -0.8417     2.3826

Subject:  Variable 10 (Level  2) Variable 11 (Level  2) Variable 12 (Level  1)
Variable  3 (Level  1)         4.2202     3.6675
Variable  3 (Level  2)        -4.3883     3.4424
Variable  4 (Level  1)        -1.1391     3.2187
Variable  4 (Level  2)         1.0814     3.0654

Subject:  Variable 10 (Level  1) Variable 11 (Level  3) Variable 12 (Level  1)
Variable  3 (Level  1)         0.3391     4.0647
Variable  3 (Level  2)         0.1502     3.4787
Variable  4 (Level  1)        -1.0026     2.4363

Subject:  Variable 10 (Level  1) Variable 11 (Level  3) Variable 12 (Level  1)
Variable  4 (Level  3)         1.1703     2.6365

Subject:  Variable 10 (Level  2) Variable 11 (Level  3) Variable 12 (Level  1)
Variable  3 (Level  1)         1.2658     3.4819
Variable  3 (Level  2)        -1.5356     3.9097

Subject:  Variable 10 (Level  2) Variable 11 (Level  3) Variable 12 (Level  1)
Variable  4 (Level  2)         0.7992     2.7902
Variable  4 (Level  3)        -0.8916     2.8763

Subject:  Variable 11 (Level  1) Variable 12 (Level  1)
Variable  5 (Level  1)        -0.4885     2.8206
Variable  5 (Level  2)         1.8829     2.7530
Variable  6 (Level  1)         0.9249     3.7747
Variable  6 (Level  2)        -2.3568     3.1624
Variable  6 (Level  3)         4.3117     3.1474

Subject:  Variable 11 (Level  2) Variable 12 (Level  1)
Variable  5 (Level  1)         1.3898     2.9362
Variable  5 (Level  2)        -1.5729     2.8909
Variable  6 (Level  1)         0.2111     3.9967
Variable  6 (Level  2)        -3.7083     4.2866
Variable  6 (Level  3)         3.1190     4.7983

Subject:  Variable 11 (Level  3) Variable 12 (Level  1)
Variable  5 (Level  1)         1.7352     3.1370
Variable  5 (Level  2)        -1.6165     3.1713
Variable  6 (Level  1)        -1.1102     3.9374
Variable  6 (Level  2)         4.4877     3.6980
Variable  6 (Level  3)        -3.1325     3.1966

Subject:  Variable 12 (Level  1)
Variable  7                    0.6827     0.5060
Variable  8 (Level  1)         1.5964     1.3206
Variable  8 (Level  2)        -0.7533     1.5663
Variable  8 (Level  3)         0.4035     1.6840
Variable  8 (Level  4)        -0.8523     1.7518
Variable  9 (Level  1)         0.5699     1.6236
Variable  9 (Level  2)         0.0012     1.9111
Variable  9 (Level  3)        -0.2850     1.9245
Variable  9 (Level  4)         0.4468     2.0329
Variable  9 (Level  5)         0.0030     2.1390

Subject:  Variable 10 (Level  1) Variable 11 (Level  1) Variable 12 (Level  2)
Variable  3 (Level  1)         6.2551     3.3595
Variable  3 (Level  2)         5.6085     3.4127

Subject:  Variable 10 (Level  1) Variable 11 (Level  1) Variable 12 (Level  2)
Variable  4 (Level  2)         2.6922     2.7542
Variable  4 (Level  3)         1.3742     2.8068

Subject:  Variable 10 (Level  2) Variable 11 (Level  1) Variable 12 (Level  2)
Variable  3 (Level  1)         1.5647     3.8353
Variable  3 (Level  2)        -2.7565     3.9041
Variable  4 (Level  1)        -0.8621     2.8257

Subject:  Variable 10 (Level  2) Variable 11 (Level  1) Variable 12 (Level  2)
Variable  4 (Level  3)         0.4536     2.8070

Subject:  Variable 10 (Level  1) Variable 11 (Level  2) Variable 12 (Level  2)
Variable  3 (Level  1)       -10.1544     3.3433
Variable  3 (Level  2)         3.2446     4.1221
Variable  4 (Level  1)        -2.9419     2.3508
Variable  4 (Level  2)         0.2510     3.0675
Variable  4 (Level  3)         0.3224     2.9710

Subject:  Variable 10 (Level  2) Variable 11 (Level  2) Variable 12 (Level  2)
Variable  3 (Level  1)        -1.3577     3.1925
Variable  3 (Level  2)         8.1277     3.9975
Variable  4 (Level  1)        -0.4290     2.4578
Variable  4 (Level  2)         2.7495     2.5821

Subject:  Variable 10 (Level  1) Variable 11 (Level  3) Variable 12 (Level  2)
Variable  3 (Level  1)         4.8432     4.0069
Variable  3 (Level  2)         0.0370     3.6006
Variable  4 (Level  1)         3.0713     2.2706
Variable  4 (Level  2)        -1.8899     2.4756
Variable  4 (Level  3)         0.4914     2.2914

Subject:  Variable 10 (Level  2) Variable 11 (Level  3) Variable 12 (Level  2)
Variable  3 (Level  1)        -4.4766     3.3355
Variable  3 (Level  2)        -3.7936     4.0759
Variable  4 (Level  1)        -0.5459     2.7097
Variable  4 (Level  2)        -1.5619     2.7412
Variable  4 (Level  3)        -0.7269     2.9735

Subject:  Variable 11 (Level  1) Variable 12 (Level  2)
Variable  5 (Level  1)         4.8653     3.0706
Variable  5 (Level  2)         0.9011     3.0696
Variable  6 (Level  1)         6.9277     3.8411
Variable  6 (Level  2)        -1.3108     3.1667
Variable  6 (Level  3)         6.2916     3.5327

Subject:  Variable 11 (Level  2) Variable 12 (Level  2)
Variable  5 (Level  1)        -0.4047     3.0956
Variable  5 (Level  2)         0.3291     3.0784
Variable  6 (Level  1)         6.9096     3.3073
Variable  6 (Level  2)        -1.0680     3.6213
Variable  6 (Level  3)        -5.9977     3.7299

Subject:  Variable 11 (Level  3) Variable 12 (Level  2)
Variable  5 (Level  1)        -1.0925     3.0994
Variable  5 (Level  2)        -0.7392     2.9900
Variable  6 (Level  1)         2.7758     3.8748
Variable  6 (Level  2)        -6.3526     3.3014
Variable  6 (Level  3)        -0.2060     3.6481

Subject:  Variable 12 (Level  2)
Variable  7                    0.1711     0.5785
Variable  8 (Level  1)         1.7186     1.9143
Variable  8 (Level  2)        -0.6768     1.7352
Variable  8 (Level  3)        -0.0439     1.6395
Variable  8 (Level  4)         0.1463     1.5358
Variable  9 (Level  1)         0.9761     2.3930
Variable  9 (Level  2)         6.5436     1.8193
Variable  9 (Level  3)        -1.5504     1.8527
Variable  9 (Level  4)         0.1047     2.0244
Variable  9 (Level  5)        -3.9386     1.7937

Subject:  Variable 10 (Level  1) Variable 11 (Level  1) Variable 12 (Level  3)
Variable  3 (Level  1)        10.6802     3.2596
Variable  3 (Level  2)        -1.0290     3.7842
Variable  4 (Level  1)        -2.8612     2.2917
Variable  4 (Level  2)         3.9265     2.8934
Variable  4 (Level  3)         2.2427     2.3737

Subject:  Variable 10 (Level  2) Variable 11 (Level  1) Variable 12 (Level  3)
Variable  3 (Level  1)        -6.2076     3.3642
Variable  3 (Level  2)        -8.7670     3.8463
Variable  4 (Level  1)        -2.9251     2.4657

Subject:  Variable 10 (Level  2) Variable 11 (Level  1) Variable 12 (Level  3)
Variable  4 (Level  3)        -2.2077     2.3743

Subject:  Variable 10 (Level  1) Variable 11 (Level  2) Variable 12 (Level  3)
Variable  3 (Level  1)        -3.3334     3.4665
Variable  3 (Level  2)        -0.3111     3.2650
Variable  4 (Level  1)         1.5131     2.4890
Variable  4 (Level  2)        -3.0345     3.0562
Variable  4 (Level  3)         0.2722     2.8300

Subject:  Variable 10 (Level  2) Variable 11 (Level  2) Variable 12 (Level  3)
Variable  3 (Level  1)         6.5905     4.0386
Variable  3 (Level  2)        -5.3168     3.4549
Variable  4 (Level  1)        -3.5280     2.9663
Variable  4 (Level  2)         1.7056     2.9293
Variable  4 (Level  3)         2.2590     3.1780

Subject:  Variable 10 (Level  1) Variable 11 (Level  3) Variable 12 (Level  3)
Variable  3 (Level  1)         8.1889     4.1429
Variable  3 (Level  2)        -1.5388     3.3333
Variable  4 (Level  1)         3.4338     2.6376

Subject:  Variable 10 (Level  1) Variable 11 (Level  3) Variable 12 (Level  3)
Variable  4 (Level  3)        -1.1544     2.9885

Subject:  Variable 10 (Level  2) Variable 11 (Level  3) Variable 12 (Level  3)
Variable  3 (Level  1)        -4.4243     4.0049
Variable  3 (Level  2)        -4.1349     3.1248
Variable  4 (Level  1)         1.0460     2.6550
Variable  4 (Level  2)        -4.4844     2.2843
Variable  4 (Level  3)         0.5046     2.6926

Subject:  Variable 11 (Level  1) Variable 12 (Level  3)
Variable  5 (Level  1)         5.3030     3.0278
Variable  5 (Level  2)        -8.1794     3.1335
Variable  6 (Level  1)        -0.8188     3.7810
Variable  6 (Level  2)        -2.5078     3.1514
Variable  6 (Level  3)        -2.6138     3.4600

Subject:  Variable 11 (Level  2) Variable 12 (Level  3)
Variable  5 (Level  1)         4.3331     3.1489
Variable  5 (Level  2)        -5.6142     3.1649
Variable  6 (Level  1)        -5.8804     3.1770
Variable  6 (Level  2)         5.4265     3.3006
Variable  6 (Level  3)        -2.1917     3.2156

Subject:  Variable 11 (Level  3) Variable 12 (Level  3)
Variable  5 (Level  1)         0.4305     2.9144
Variable  5 (Level  2)        -1.4620     3.0119
Variable  6 (Level  1)        14.3595     3.9254
Variable  6 (Level  2)        -5.2399     3.3099
Variable  6 (Level  3)       -11.2498     3.2212

Subject:  Variable 12 (Level  3)
Variable  7                   -0.3839     0.6755
Variable  8 (Level  1)         2.7549     1.6017
Variable  8 (Level  2)         0.4377     1.8826
Variable  8 (Level  3)        -0.2261     1.9909
Variable  8 (Level  4)        -4.5051     1.5398
Variable  9 (Level  1)        -4.7091     2.1458
Variable  9 (Level  2)         3.7940     1.9872
Variable  9 (Level  3)        -1.7994     1.8614
Variable  9 (Level  4)         0.4480     1.9016
Variable  9 (Level  5)        -0.6047     2.4729

Fixed Effects
Intercept                      1.6433     2.4596
Variable  1 (Level  2)        -1.6224     0.8549
Variable  2 (Level  2)        -2.4817     1.1414
Variable  2 (Level  3)         0.4624     1.2133

Variance Components
Estimate     Parameter       Subject
36.32491     Variable  3     Variables 10 11 12
12.45090     Variable  4     Variables 10 11 12
19.62767     Variable  5     Variables 11 12
40.53480     Variable  6     Variables 11 12
0.56320     Variable  7     Variables 12
5.81968     Variable  8     Variables 12
10.86069     Variable  9     Variables 12

sigma^2          =         0.00239
-2log likelihood =       608.19449

```
```function g02jc_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] = ...
g02jd(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.1561     3.7946
Variable  3 (Level  2)         1.8951     3.9284
Variable  4 (Level  1)         0.6496     3.1617

Subject:  Variable 10 (Level  1) Variable 11 (Level  1) Variable 12 (Level  1)
Variable  4 (Level  3)         0.7390     3.1424

Subject:  Variable 10 (Level  2) Variable 11 (Level  1) Variable 12 (Level  1)
Variable  3 (Level  1)         1.4216     3.3773
Variable  3 (Level  2)        -2.8921     3.3953
Variable  4 (Level  1)         3.6789     2.3162
Variable  4 (Level  2)        -1.9742     2.3887
Variable  4 (Level  3)        -2.2088     2.0697

Subject:  Variable 10 (Level  1) Variable 11 (Level  2) Variable 12 (Level  1)
Variable  3 (Level  1)        -2.9659     3.9127
Variable  3 (Level  2)         2.7951     4.7183
Variable  4 (Level  1)        -4.7330     2.3094
Variable  4 (Level  2)         5.5161     2.2330
Variable  4 (Level  3)        -0.8417     2.3826

Subject:  Variable 10 (Level  2) Variable 11 (Level  2) Variable 12 (Level  1)
Variable  3 (Level  1)         4.2202     3.6675
Variable  3 (Level  2)        -4.3883     3.4424
Variable  4 (Level  1)        -1.1391     3.2187
Variable  4 (Level  2)         1.0814     3.0654

Subject:  Variable 10 (Level  1) Variable 11 (Level  3) Variable 12 (Level  1)
Variable  3 (Level  1)         0.3391     4.0647
Variable  3 (Level  2)         0.1502     3.4787
Variable  4 (Level  1)        -1.0026     2.4363

Subject:  Variable 10 (Level  1) Variable 11 (Level  3) Variable 12 (Level  1)
Variable  4 (Level  3)         1.1703     2.6365

Subject:  Variable 10 (Level  2) Variable 11 (Level  3) Variable 12 (Level  1)
Variable  3 (Level  1)         1.2658     3.4819
Variable  3 (Level  2)        -1.5356     3.9097

Subject:  Variable 10 (Level  2) Variable 11 (Level  3) Variable 12 (Level  1)
Variable  4 (Level  2)         0.7992     2.7902
Variable  4 (Level  3)        -0.8916     2.8763

Subject:  Variable 11 (Level  1) Variable 12 (Level  1)
Variable  5 (Level  1)        -0.4885     2.8206
Variable  5 (Level  2)         1.8829     2.7530
Variable  6 (Level  1)         0.9249     3.7747
Variable  6 (Level  2)        -2.3568     3.1624
Variable  6 (Level  3)         4.3117     3.1474

Subject:  Variable 11 (Level  2) Variable 12 (Level  1)
Variable  5 (Level  1)         1.3898     2.9362
Variable  5 (Level  2)        -1.5729     2.8909
Variable  6 (Level  1)         0.2111     3.9967
Variable  6 (Level  2)        -3.7083     4.2866
Variable  6 (Level  3)         3.1190     4.7983

Subject:  Variable 11 (Level  3) Variable 12 (Level  1)
Variable  5 (Level  1)         1.7352     3.1370
Variable  5 (Level  2)        -1.6165     3.1713
Variable  6 (Level  1)        -1.1102     3.9374
Variable  6 (Level  2)         4.4877     3.6980
Variable  6 (Level  3)        -3.1325     3.1966

Subject:  Variable 12 (Level  1)
Variable  7                    0.6827     0.5060
Variable  8 (Level  1)         1.5964     1.3206
Variable  8 (Level  2)        -0.7533     1.5663
Variable  8 (Level  3)         0.4035     1.6840
Variable  8 (Level  4)        -0.8523     1.7518
Variable  9 (Level  1)         0.5699     1.6236
Variable  9 (Level  2)         0.0012     1.9111
Variable  9 (Level  3)        -0.2850     1.9245
Variable  9 (Level  4)         0.4468     2.0329
Variable  9 (Level  5)         0.0030     2.1390

Subject:  Variable 10 (Level  1) Variable 11 (Level  1) Variable 12 (Level  2)
Variable  3 (Level  1)         6.2551     3.3595
Variable  3 (Level  2)         5.6085     3.4127

Subject:  Variable 10 (Level  1) Variable 11 (Level  1) Variable 12 (Level  2)
Variable  4 (Level  2)         2.6922     2.7542
Variable  4 (Level  3)         1.3742     2.8068

Subject:  Variable 10 (Level  2) Variable 11 (Level  1) Variable 12 (Level  2)
Variable  3 (Level  1)         1.5647     3.8353
Variable  3 (Level  2)        -2.7565     3.9041
Variable  4 (Level  1)        -0.8621     2.8257

Subject:  Variable 10 (Level  2) Variable 11 (Level  1) Variable 12 (Level  2)
Variable  4 (Level  3)         0.4536     2.8070

Subject:  Variable 10 (Level  1) Variable 11 (Level  2) Variable 12 (Level  2)
Variable  3 (Level  1)       -10.1544     3.3433
Variable  3 (Level  2)         3.2446     4.1221
Variable  4 (Level  1)        -2.9419     2.3508
Variable  4 (Level  2)         0.2510     3.0675
Variable  4 (Level  3)         0.3224     2.9710

Subject:  Variable 10 (Level  2) Variable 11 (Level  2) Variable 12 (Level  2)
Variable  3 (Level  1)        -1.3577     3.1925
Variable  3 (Level  2)         8.1277     3.9975
Variable  4 (Level  1)        -0.4290     2.4578
Variable  4 (Level  2)         2.7495     2.5821

Subject:  Variable 10 (Level  1) Variable 11 (Level  3) Variable 12 (Level  2)
Variable  3 (Level  1)         4.8432     4.0069
Variable  3 (Level  2)         0.0370     3.6006
Variable  4 (Level  1)         3.0713     2.2706
Variable  4 (Level  2)        -1.8899     2.4756
Variable  4 (Level  3)         0.4914     2.2914

Subject:  Variable 10 (Level  2) Variable 11 (Level  3) Variable 12 (Level  2)
Variable  3 (Level  1)        -4.4766     3.3355
Variable  3 (Level  2)        -3.7936     4.0759
Variable  4 (Level  1)        -0.5459     2.7097
Variable  4 (Level  2)        -1.5619     2.7412
Variable  4 (Level  3)        -0.7269     2.9735

Subject:  Variable 11 (Level  1) Variable 12 (Level  2)
Variable  5 (Level  1)         4.8653     3.0706
Variable  5 (Level  2)         0.9011     3.0696
Variable  6 (Level  1)         6.9277     3.8411
Variable  6 (Level  2)        -1.3108     3.1667
Variable  6 (Level  3)         6.2916     3.5327

Subject:  Variable 11 (Level  2) Variable 12 (Level  2)
Variable  5 (Level  1)        -0.4047     3.0956
Variable  5 (Level  2)         0.3291     3.0784
Variable  6 (Level  1)         6.9096     3.3073
Variable  6 (Level  2)        -1.0680     3.6213
Variable  6 (Level  3)        -5.9977     3.7299

Subject:  Variable 11 (Level  3) Variable 12 (Level  2)
Variable  5 (Level  1)        -1.0925     3.0994
Variable  5 (Level  2)        -0.7392     2.9900
Variable  6 (Level  1)         2.7758     3.8748
Variable  6 (Level  2)        -6.3526     3.3014
Variable  6 (Level  3)        -0.2060     3.6481

Subject:  Variable 12 (Level  2)
Variable  7                    0.1711     0.5785
Variable  8 (Level  1)         1.7186     1.9143
Variable  8 (Level  2)        -0.6768     1.7352
Variable  8 (Level  3)        -0.0439     1.6395
Variable  8 (Level  4)         0.1463     1.5358
Variable  9 (Level  1)         0.9761     2.3930
Variable  9 (Level  2)         6.5436     1.8193
Variable  9 (Level  3)        -1.5504     1.8527
Variable  9 (Level  4)         0.1047     2.0244
Variable  9 (Level  5)        -3.9386     1.7937

Subject:  Variable 10 (Level  1) Variable 11 (Level  1) Variable 12 (Level  3)
Variable  3 (Level  1)        10.6802     3.2596
Variable  3 (Level  2)        -1.0290     3.7842
Variable  4 (Level  1)        -2.8612     2.2917
Variable  4 (Level  2)         3.9265     2.8934
Variable  4 (Level  3)         2.2427     2.3737

Subject:  Variable 10 (Level  2) Variable 11 (Level  1) Variable 12 (Level  3)
Variable  3 (Level  1)        -6.2076     3.3642
Variable  3 (Level  2)        -8.7670     3.8463
Variable  4 (Level  1)        -2.9251     2.4657

Subject:  Variable 10 (Level  2) Variable 11 (Level  1) Variable 12 (Level  3)
Variable  4 (Level  3)        -2.2077     2.3743

Subject:  Variable 10 (Level  1) Variable 11 (Level  2) Variable 12 (Level  3)
Variable  3 (Level  1)        -3.3334     3.4665
Variable  3 (Level  2)        -0.3111     3.2650
Variable  4 (Level  1)         1.5131     2.4890
Variable  4 (Level  2)        -3.0345     3.0562
Variable  4 (Level  3)         0.2722     2.8300

Subject:  Variable 10 (Level  2) Variable 11 (Level  2) Variable 12 (Level  3)
Variable  3 (Level  1)         6.5905     4.0386
Variable  3 (Level  2)        -5.3168     3.4549
Variable  4 (Level  1)        -3.5280     2.9663
Variable  4 (Level  2)         1.7056     2.9293
Variable  4 (Level  3)         2.2590     3.1780

Subject:  Variable 10 (Level  1) Variable 11 (Level  3) Variable 12 (Level  3)
Variable  3 (Level  1)         8.1889     4.1429
Variable  3 (Level  2)        -1.5388     3.3333
Variable  4 (Level  1)         3.4338     2.6376

Subject:  Variable 10 (Level  1) Variable 11 (Level  3) Variable 12 (Level  3)
Variable  4 (Level  3)        -1.1544     2.9885

Subject:  Variable 10 (Level  2) Variable 11 (Level  3) Variable 12 (Level  3)
Variable  3 (Level  1)        -4.4243     4.0049
Variable  3 (Level  2)        -4.1349     3.1248
Variable  4 (Level  1)         1.0460     2.6550
Variable  4 (Level  2)        -4.4844     2.2843
Variable  4 (Level  3)         0.5046     2.6926

Subject:  Variable 11 (Level  1) Variable 12 (Level  3)
Variable  5 (Level  1)         5.3030     3.0278
Variable  5 (Level  2)        -8.1794     3.1335
Variable  6 (Level  1)        -0.8188     3.7810
Variable  6 (Level  2)        -2.5078     3.1514
Variable  6 (Level  3)        -2.6138     3.4600

Subject:  Variable 11 (Level  2) Variable 12 (Level  3)
Variable  5 (Level  1)         4.3331     3.1489
Variable  5 (Level  2)        -5.6142     3.1649
Variable  6 (Level  1)        -5.8804     3.1770
Variable  6 (Level  2)         5.4265     3.3006
Variable  6 (Level  3)        -2.1917     3.2156

Subject:  Variable 11 (Level  3) Variable 12 (Level  3)
Variable  5 (Level  1)         0.4305     2.9144
Variable  5 (Level  2)        -1.4620     3.0119
Variable  6 (Level  1)        14.3595     3.9254
Variable  6 (Level  2)        -5.2399     3.3099
Variable  6 (Level  3)       -11.2498     3.2212

Subject:  Variable 12 (Level  3)
Variable  7                   -0.3839     0.6755
Variable  8 (Level  1)         2.7549     1.6017
Variable  8 (Level  2)         0.4377     1.8826
Variable  8 (Level  3)        -0.2261     1.9909
Variable  8 (Level  4)        -4.5051     1.5398
Variable  9 (Level  1)        -4.7091     2.1458
Variable  9 (Level  2)         3.7940     1.9872
Variable  9 (Level  3)        -1.7994     1.8614
Variable  9 (Level  4)         0.4480     1.9016
Variable  9 (Level  5)        -0.6047     2.4729

Fixed Effects
Intercept                      1.6433     2.4596
Variable  1 (Level  2)        -1.6224     0.8549
Variable  2 (Level  2)        -2.4817     1.1414
Variable  2 (Level  3)         0.4624     1.2133

Variance Components
Estimate     Parameter       Subject
36.32491     Variable  3     Variables 10 11 12
12.45090     Variable  4     Variables 10 11 12
19.62767     Variable  5     Variables 11 12
40.53480     Variable  6     Variables 11 12
0.56320     Variable  7     Variables 12
5.81968     Variable  8     Variables 12
10.86069     Variable  9     Variables 12

sigma^2          =         0.00239
-2log likelihood =       608.19449

```