g02jb fits a linear mixed effects regression model using maximum likelihood (ML).

# Syntax

C#
```public static void g02jb(
int n,
int ncol,
double[,] dat,
int[] levels,
int yvid,
int cwid,
int nfv,
int[] fvid,
int fint,
int nrv,
int[] rvid,
int nvpr,
int[] vpr,
int rint,
int svid,
double[] gamma,
out int nff,
out int nrf,
out int df,
out double ml,
double[] b,
double[] se,
int maxit,
double tol,
out int warn,
out int ifail
)```
Visual Basic
```Public Shared Sub g02jb ( _
n As Integer, _
ncol As Integer, _
dat As Double(,), _
levels As Integer(), _
yvid As Integer, _
cwid As Integer, _
nfv As Integer, _
fvid As Integer(), _
fint As Integer, _
nrv As Integer, _
rvid As Integer(), _
nvpr As Integer, _
vpr As Integer(), _
rint As Integer, _
svid As Integer, _
gamma As Double(), _
<OutAttribute> ByRef nff As Integer, _
<OutAttribute> ByRef nrf As Integer, _
<OutAttribute> ByRef df As Integer, _
<OutAttribute> ByRef ml As Double, _
b As Double(), _
se As Double(), _
maxit As Integer, _
tol As Double, _
<OutAttribute> ByRef warn As Integer, _
<OutAttribute> ByRef ifail As Integer _
)```
Visual C++
```public:
static void g02jb(
int n,
int ncol,
array<double,2>^ dat,
array<int>^ levels,
int yvid,
int cwid,
int nfv,
array<int>^ fvid,
int fint,
int nrv,
array<int>^ rvid,
int nvpr,
array<int>^ vpr,
int rint,
int svid,
array<double>^ gamma,
[OutAttribute] int% nff,
[OutAttribute] int% nrf,
[OutAttribute] int% df,
[OutAttribute] double% ml,
array<double>^ b,
array<double>^ se,
int maxit,
double tol,
[OutAttribute] int% warn,
[OutAttribute] int% ifail
)```
F#
```static member g02jb :
n : int *
ncol : int *
dat : float[,] *
levels : int[] *
yvid : int *
cwid : int *
nfv : int *
fvid : int[] *
fint : int *
nrv : int *
rvid : int[] *
nvpr : int *
vpr : int[] *
rint : int *
svid : int *
gamma : float[] *
nff : int byref *
nrf : int byref *
df : int byref *
ml : float byref *
b : float[] *
se : float[] *
maxit : int *
tol : float *
warn : int byref *
ifail : int byref -> unit
```

#### Parameters

n
Type: System..::..Int32
On entry: $n$, the number of observations.
Constraint: ${\mathbf{n}}\ge 1$.
ncol
Type: System..::..Int32
On entry: the number of columns in the data matrix, dat.
Constraint: ${\mathbf{ncol}}\ge 1$.
dat
Type: array<System..::..Double,2>[,](,)[,][,]
An array of size [dim1, dim2]
Note: dim1 must satisfy the constraint: $\mathrm{dim1}\ge {\mathbf{n}}$
Note: the second dimension of the array dat must be at least ${\mathbf{ncol}}$ if $\mathbf{_sorder}=1$, and at least ${\mathbf{n}}$ otherwise.
On entry: array containing all of the data. For the $i$th observation:
• ${\mathbf{dat}}\left[i-1,{\mathbf{yvid}}-1\right]$ holds the dependent variable, $y$;
• if ${\mathbf{cwid}}\ne 0$, ${\mathbf{dat}}\left[i-1,{\mathbf{cwid}}-1\right]$ holds the case weights;
• if ${\mathbf{svid}}\ne 0$, ${\mathbf{dat}}\left[i-1,{\mathbf{svid}}-1\right]$ holds the subject variable.
The remaining columns hold the values of the independent variables.
Constraints:
• if ${\mathbf{cwid}}\ne 0$, ${\mathbf{dat}}\left[i-1,{\mathbf{cwid}}-1\right]\ge 0.0$;
• if ${\mathbf{levels}}\left[j-1\right]\ne 1$, $1\le {\mathbf{dat}}\left[i-1,j-1\right]\le {\mathbf{levels}}\left[j-1\right]$.
levels
Type: array<System..::..Int32>[]()[][]
An array of size [ncol]
On entry: ${\mathbf{levels}}\left[i-1\right]$ contains the number of levels associated with the $i$th variable of the data matrix dat. If this variable is continuous or binary (i.e., only takes the values zero or one) then ${\mathbf{levels}}\left[i-1\right]$ should be $1$; if the variable is discrete then ${\mathbf{levels}}\left[i-1\right]$ is the number of levels associated with it and ${\mathbf{dat}}\left[\mathit{j}-1,i-1\right]$ is assumed to take the values $1$ to ${\mathbf{levels}}\left[i-1\right]$, for $\mathit{j}=1,2,\dots ,{\mathbf{n}}$.
Constraint: ${\mathbf{levels}}\left[\mathit{i}-1\right]\ge 1$, for $\mathit{i}=1,2,\dots ,{\mathbf{ncol}}$.
yvid
Type: System..::..Int32
On entry: the column of dat holding the dependent, $y$, variable.
Constraint: $1\le {\mathbf{yvid}}\le {\mathbf{ncol}}$.
cwid
Type: System..::..Int32
On entry: the column of dat holding the case weights.
If ${\mathbf{cwid}}=0$, no weights are used.
Constraint: $0\le {\mathbf{cwid}}\le {\mathbf{ncol}}$.
nfv
Type: System..::..Int32
On entry: the number of independent variables in the model which are to be treated as being fixed.
Constraint: $0\le {\mathbf{nfv}}<{\mathbf{ncol}}$.
fvid
Type: array<System..::..Int32>[]()[][]
An array of size [nfv]
On entry: the columns of the data matrix dat holding the fixed independent variables with ${\mathbf{fvid}}\left[i-1\right]$ holding the column number corresponding to the $i$th fixed variable.
Constraint: $1\le {\mathbf{fvid}}\left[\mathit{i}-1\right]\le {\mathbf{ncol}}$, for $\mathit{i}=1,2,\dots ,{\mathbf{nfv}}$.
fint
Type: System..::..Int32
On entry: flag indicating whether a fixed intercept is included (${\mathbf{fint}}=1$).
Constraint: ${\mathbf{fint}}=0$ or $1$.
nrv
Type: System..::..Int32
On entry: the number of independent variables in the model which are to be treated as being random.
Constraints:
• $0\le {\mathbf{nrv}}<{\mathbf{ncol}}$;
• ${\mathbf{nrv}}+{\mathbf{rint}}>0$.
rvid
Type: array<System..::..Int32>[]()[][]
An array of size [nrv]
On entry: the columns of the data matrix ${\mathbf{dat}}$ holding the random independent variables with ${\mathbf{rvid}}\left[i-1\right]$ holding the column number corresponding to the $i$th random variable.
Constraint: $1\le {\mathbf{rvid}}\left[\mathit{i}-1\right]\le {\mathbf{ncol}}$, for $\mathit{i}=1,2,\dots ,{\mathbf{nrv}}$.
nvpr
Type: System..::..Int32
On entry: if ${\mathbf{rint}}=1$ and ${\mathbf{svid}}\ne 0$, nvpr is the number of variance components being $\text{estimated}-2$, ($g-1$), else ${\mathbf{nvpr}}=g$.
If ${\mathbf{nrv}}=0$, ${\mathbf{nvpr}}$ is not referenced.
Constraint: if ${\mathbf{nrv}}\ne 0$, $1\le {\mathbf{nvpr}}\le {\mathbf{nrv}}$.
vpr
Type: array<System..::..Int32>[]()[][]
An array of size [nrv]
On entry: ${\mathbf{vpr}}\left[i-1\right]$ holds a flag indicating the variance of the $i$th random variable. The variance of the $i$th random variable is ${\sigma }_{j}^{2}$, where $j={\mathbf{vpr}}\left[i-1\right]+1$ if ${\mathbf{rint}}=1$ and ${\mathbf{svid}}\ne 0$ and $j={\mathbf{vpr}}\left[i-1\right]$ otherwise. Random variables with the same value of $j$ are assumed to be taken from the same distribution.
Constraint: $1\le {\mathbf{vpr}}\left[\mathit{i}-1\right]\le {\mathbf{nvpr}}$, for $\mathit{i}=1,2,\dots ,{\mathbf{nrv}}$.
rint
Type: System..::..Int32
On entry: flag indicating whether a random intercept is included (${\mathbf{rint}}=1$).
If ${\mathbf{svid}}=0$, rint is not referenced.
Constraint: ${\mathbf{rint}}=0$ or $1$.
svid
Type: System..::..Int32
On entry: the column of dat holding the subject variable.
If ${\mathbf{svid}}=0$, no subject variable is used.
Specifying a subject variable is equivalent to specifying the interaction between that variable and all of the random-effects. Letting the notation ${Z}_{1}×{Z}_{S}$ denote the interaction between variables ${Z}_{1}$ and ${Z}_{S}$, fitting a model with ${\mathbf{rint}}=0$, random-effects ${Z}_{1}+{Z}_{2}$ and subject variable ${Z}_{S}$ is equivalent to fitting a model with random-effects ${Z}_{1}×{Z}_{S}+{Z}_{2}×{Z}_{S}$ and no subject variable. If ${\mathbf{rint}}=1$ the model is equivalent to fitting ${Z}_{S}+{Z}_{1}×{Z}_{S}+{Z}_{2}×{Z}_{S}$ and no subject variable.
Constraint: $0\le {\mathbf{svid}}\le {\mathbf{ncol}}$.
gamma
Type: array<System..::..Double>[]()[][]
An array of size [${\mathbf{nvpr}}+2$]
On entry: holds the initial values of the variance components, ${\gamma }_{0}$, with ${\mathbf{gamma}}\left[\mathit{i}-1\right]$ the initial value for ${\sigma }_{\mathit{i}}^{2}/{\sigma }_{R}^{2}$, for $\mathit{i}=1,2,\dots ,g$. If ${\mathbf{rint}}=1$ and ${\mathbf{svid}}\ne 0$, $g={\mathbf{nvpr}}+1$, else $g={\mathbf{nvpr}}$.
If ${\mathbf{gamma}}\left[0\right]=-1.0$, the remaining elements of gamma are ignored and the initial values for the variance components are estimated from the data using MIVQUE0.
On exit: ${\mathbf{gamma}}\left[\mathit{i}-1\right]$, for $\mathit{i}=1,2,\dots ,g$, holds the final estimate of ${\sigma }_{\mathit{i}}^{2}$ and ${\mathbf{gamma}}\left[g\right]$ holds the final estimate for ${\sigma }_{R}^{2}$.
Constraint: ${\mathbf{gamma}}\left[0\right]=-1.0\text{​ or ​}{\mathbf{gamma}}\left[\mathit{i}-1\right]\ge 0.0$, for $\mathit{i}=1,2,\dots ,g$.
nff
Type: System..::..Int32%
On exit: the number of fixed effects estimated (i.e., the number of columns, $p$, in the design matrix $X$).
nrf
Type: System..::..Int32%
On exit: the number of random effects estimated (i.e., the number of columns, $q$, in the design matrix $Z$).
df
Type: System..::..Int32%
On exit: the degrees of freedom.
ml
Type: System..::..Double%
On exit: $-2{l}_{R}\left(\stackrel{^}{\gamma }\right)$ where ${l}_{R}$ is the log of the maximum likelihood calculated at $\stackrel{^}{\gamma }$, the estimated variance components returned in gamma.
b
Type: array<System..::..Double>[]()[][]
An array of size [dim1]
Note: dim1 must satisfy the constraint: $\mathbf{_lb}\ge {\mathbf{fint}}+\sum _{i=1}^{{\mathbf{nfv}}}\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{levels}}\left[{\mathbf{fvid}}\left[i-1\right]-1\right]-1,1\right)+{L}_{S}×\left({\mathbf{rint}}+\sum _{i=1}^{{\mathbf{nrv}}}{\mathbf{levels}}\left[{\mathbf{rvid}}\left[i-1\right]-1\right]\right)$ where ${L}_{S}={\mathbf{levels}}\left[{\mathbf{svid}}-1\right]$ if ${\mathbf{svid}}\ne 0$ and $1$ otherwise
On exit: the parameter estimates, $\left(\beta ,\nu \right)$, with the first nff elements of b containing the fixed effect parameter estimates, $\beta$ and the next nrf elements of b containing the random effect parameter estimates, $\nu$.
Fixed effects
If ${\mathbf{fint}}=1$, ${\mathbf{b}}\left[0\right]$ contains the estimate of the fixed intercept. Let ${L}_{i}$ denote the number of levels associated with the $i$th fixed variable, that is ${L}_{i}={\mathbf{levels}}\left[{\mathbf{fvid}}\left[i-1\right]-1\right]$. Define
• if ${\mathbf{fint}}=1$, ${F}_{1}=2$ else if ${\mathbf{fint}}=0$, ${F}_{1}=1$;
• ${F}_{i+1}={F}_{i}+\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({L}_{i}-1,1\right)$, $i\ge 1$.
Then for $i=1,2,\dots ,{\mathbf{nfv}}$:
• if ${L}_{i}>1$, ${\mathbf{b}}\left[{F}_{i}+\mathit{j}-3\right]$ contains the parameter estimate for the $\mathit{j}$th level of the $i$th fixed variable, for $\mathit{j}=2,3,\dots ,{L}_{i}$;
• if ${L}_{i}\le 1$, ${\mathbf{b}}\left[{F}_{i}-1\right]$ contains the parameter estimate for the $i$th fixed variable.
Random effects
Redefining ${L}_{i}$ to denote the number of levels associated with the $i$th random variable, that is ${L}_{i}={\mathbf{levels}}\left[{\mathbf{rvid}}\left[i-1\right]-1\right]$. Define
• if ${\mathbf{rint}}=1$, ${R}_{1}=2$ else if ${\mathbf{rint}}=0$, ${R}_{1}=1$;
${R}_{i+1}={R}_{i}+{L}_{i}$, $i\ge 1$.
Then for $i=1,2,\dots ,{\mathbf{nrv}}$:
• if ${\mathbf{svid}}=0$,
• if ${L}_{i}>1$, ${\mathbf{b}}\left[{\mathbf{nff}}+{R}_{i}+\mathit{j}-2\right]$ contains the parameter estimate for the $\mathit{j}$th level of the $i$th random variable, for $\mathit{j}=1,2,\dots ,{L}_{i}$;
• if ${L}_{i}\le 1$, ${\mathbf{b}}\left[{\mathbf{nff}}+{R}_{i}-1\right]$ contains the parameter estimate for the $i$th random variable;
• if ${\mathbf{svid}}\ne 0$,
• let ${L}_{S}$ denote the number of levels associated with the subject variable, that is ${L}_{S}={\mathbf{levels}}\left[{\mathbf{svid}}-1\right]$;
• if ${L}_{i}>1$, ${\mathbf{b}}\left[{\mathbf{nff}}+\left(\mathit{s}-1\right){L}_{S}+{R}_{i}+\mathit{j}-2\right]$ contains the parameter estimate for the interaction between the $\mathit{s}$th level of the subject variable and the $\mathit{j}$th level of the $i$th random variable, for $\mathit{s}=1,2,\dots ,{L}_{S}$ and $\mathit{j}=1,2,\dots ,{L}_{i}$;
• if ${L}_{i}\le 1$, ${\mathbf{b}}\left[{\mathbf{nff}}+\left(\mathit{s}-1\right){L}_{S}+{R}_{i}-1\right]$ contains the parameter estimate for the interaction between the $\mathit{s}$th level of the subject variable and the $i$th random variable, for $\mathit{s}=1,2,\dots ,{L}_{S}$;
• if ${\mathbf{rint}}=1$, ${\mathbf{b}}\left[{\mathbf{nff}}\right]$ contains the estimate of the random intercept.
se
Type: array<System..::..Double>[]()[][]
An array of size [dim1]
Note: dim1 must satisfy the constraint: $\mathbf{_lb}\ge {\mathbf{fint}}+\sum _{i=1}^{{\mathbf{nfv}}}\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{levels}}\left[{\mathbf{fvid}}\left[i-1\right]-1\right]-1,1\right)+{L}_{S}×\left({\mathbf{rint}}+\sum _{i=1}^{{\mathbf{nrv}}}{\mathbf{levels}}\left[{\mathbf{rvid}}\left[i-1\right]-1\right]\right)$ where ${L}_{S}={\mathbf{levels}}\left[{\mathbf{svid}}-1\right]$ if ${\mathbf{svid}}\ne 0$ and $1$ otherwise
On exit: the standard errors of the parameter estimates given in b.
maxit
Type: System..::..Int32
On entry: the maximum number of iterations.
If ${\mathbf{maxit}}<0$, the default value of $100$ is used.
If ${\mathbf{maxit}}=0$, the parameter estimates $\left(\beta ,\nu \right)$ and corresponding standard errors are calculated based on the value of ${\gamma }_{0}$ supplied in gamma.
tol
Type: System..::..Double
On entry: the tolerance used to assess convergence.
If ${\mathbf{tol}}\le 0.0$, the default value of ${\epsilon }^{0.7}$ is used, where $\epsilon$ is the machine precision.
warn
Type: System..::..Int32%
On exit: is set to $1$ if a variance component was estimated to be a negative value during the fitting process. Otherwise warn is set to $0$.
If ${\mathbf{warn}}=1$, the negative estimate is set to zero and the estimation process allowed to continue.
ifail
Type: System..::..Int32%
On exit: ${\mathbf{ifail}}={0}$ unless the method detects an error or a warning has been flagged (see [Error Indicators and Warnings]).

# Description

g02jb fits a model of the form:
 $y=Xβ+Zν+ε$
where
• $y$ is a vector of $n$ observations on the dependent variable,
• $X$ is a known $n$ by $p$ design matrix for the fixed independent variables,
• $\beta$ is a vector of length $p$ of unknown fixed effects,
• $Z$ is a known $n$ by $q$ design matrix for the random independent variables,
• $\nu$ is a vector of length $q$ of unknown random effects;
and
• $\epsilon$ is a vector of length $n$ of unknown random errors.
Both $\nu$ and $\epsilon$ are assumed to have a Gaussian distribution with expectation zero and
 $Varνε=G00R$
where $R={\sigma }_{R}^{2}I$, $I$ is the $n×n$ identity matrix and $G$ is a diagonal matrix. It is assumed that the random variables, $Z$, can be subdivided into $g\le q$ groups with each group being identically distributed with expectations zero and variance ${\sigma }_{i}^{2}$. The diagonal elements of matrix $G$ therefore take one of the values $\left\{{\sigma }_{i}^{2}:i=1,2,\dots ,g\right\}$, depending on which group the associated random variable belongs to.
The model therefore contains three sets of unknowns, the fixed effects, $\beta$, the random effects $\nu$ and a vector of $g+1$ variance components, $\gamma$, where $\gamma =\left\{{\sigma }_{1}^{2},{\sigma }_{2}^{2},\dots ,{\sigma }_{g-1}^{2},{\sigma }_{g}^{2},{\sigma }_{R}^{2}\right\}$. Rather than working directly with $\gamma$, g02jb uses an iterative process to estimate ${\gamma }^{*}=\left\{{\sigma }_{1}^{2}/{\sigma }_{R}^{2},{\sigma }_{2}^{2}/{\sigma }_{R}^{2},\dots ,{\sigma }_{g-1}^{2}/{\sigma }_{R}^{2},{\sigma }_{g}^{2}/{\sigma }_{R}^{2},1\right\}$. Due to the iterative nature of the estimation a set of initial values, ${\gamma }_{0}$, for ${\gamma }^{*}$ is required. g02jb 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).
g02jb fits the model using a quasi-Newton algorithm to maximize the log-likelihood function:
 $-2lR=logV+nlogr′V-1r+log2π/n$
where
 $V=ZGZ′+R, r=y-Xb and b=X′V-1X-1X′V-1y.$
Once the final estimates for ${\gamma }^{*}$ have been obtained, the value of ${\sigma }_{R}^{2}$ is given by:
 $σR2=r′V-1r/n-p.$
Case weights, ${W}_{c}$, can be incorporated into the model by replacing ${X}^{\prime }X$ and ${Z}^{\prime }Z$ with ${X}^{\prime }{W}_{c}X$ and ${Z}^{\prime }{W}_{c}Z$ respectively, for a diagonal weight matrix ${W}_{c}$.
The log-likelihood, ${l}_{R}$, is calculated using the sweep algorithm detailed in Wolfinger et al. (1994).

# References

Goodnight J H (1979) A tutorial on the SWEEP operator The American Statistician 33(3) 149–158
Harville D A (1977) Maximum likelihood approaches to variance component estimation and to related problems JASA 72 320–340
Rao C R (1972) Estimation of variance and covariance components in a linear model J. Am. Stat. Assoc. 67 112–115
Stroup W W (1989) Predictable functions and prediction space in the mixed model procedure Applications of Mixed Models in Agriculture and Related Disciplines Southern Cooperative Series Bulletin No. 343 39–48
Wolfinger R, Tobias R and Sall J (1994) Computing Gaussian likelihoods and their derivatives for general linear mixed models SIAM Sci. Statist. Comput. 15 1294–1310

# Error Indicators and Warnings

Errors or warnings detected by the method:
Some error messages may refer to parameters that are dropped from this interface (LDDAT, LB) In these cases, an error in another parameter has usually caused an incorrect value to be inferred.
${\mathbf{ifail}}=1$
 On entry, ${\mathbf{n}}<2$, or ${\mathbf{ncol}}<1$, or ${\mathbf{yvid}}<1$ or ${\mathbf{yvid}}>{\mathbf{ncol}}$, or ${\mathbf{cwid}}<0$ or ${\mathbf{cwid}}>{\mathbf{ncol}}$, or ${\mathbf{nfv}}<0$ or ${\mathbf{nfv}}\ge {\mathbf{ncol}}$, or ${\mathbf{fint}}\ne 0$ and ${\mathbf{fint}}\ne 1$, or ${\mathbf{nrv}}<0$ or ${\mathbf{nrv}}>{\mathbf{ncol}}$ or ${\mathbf{nrv}}+{\mathbf{rint}}<1$, or ${\mathbf{nvpr}}<0$ or ${\mathbf{nvpr}}>{\mathbf{nrv}}$, or ${\mathbf{rint}}\ne 0$ and ${\mathbf{rint}}\ne 1$, or ${\mathbf{svid}}<0$ or ${\mathbf{svid}}>{\mathbf{ncol}}$,
${\mathbf{ifail}}=2$
 On entry, ${\mathbf{levels}}\left[i-1\right]<1$, for at least one $i$, or ${\mathbf{fvid}}\left[i-1\right]<1$, or ${\mathbf{fvid}}\left[i-1\right]>{\mathbf{ncol}}$, for at least one $i$, or ${\mathbf{rvid}}\left[i-1\right]<1$, or ${\mathbf{rvid}}\left[i-1\right]>{\mathbf{ncol}}$, for at least one $i$, or ${\mathbf{vpr}}\left[i-1\right]<1$ or ${\mathbf{vpr}}\left[i-1\right]>{\mathbf{nvpr}}$, for at least one $i$, or at least one discrete variable in array dat has a value greater than that specified in levels, or ${\mathbf{gamma}}\left[i-1\right]<0$, for at least one $i$, and ${\mathbf{gamma}}\left[0\right]\ne -1$.
${\mathbf{ifail}}=3$
Degrees of freedom $<1$. The number of parameters exceed the effective number of observations.
${\mathbf{ifail}}=4$
The method failed to converge to the specified tolerance in maxit iterations. See [Further Comments] for advice.
${\mathbf{ifail}}=-9000$
An error occured, see message report.
${\mathbf{ifail}}=-6000$
Invalid Parameters $〈\mathit{\text{value}}〉$
${\mathbf{ifail}}=-4000$
Invalid dimension for array $〈\mathit{\text{value}}〉$
${\mathbf{ifail}}=-8000$
Negative dimension for array $〈\mathit{\text{value}}〉$
${\mathbf{ifail}}=-6000$
Invalid Parameters $〈\mathit{\text{value}}〉$

# Accuracy

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

# Parallelism and Performance

None.

Wherever possible any block structure present in the design matrix $Z$ should be modelled through a subject variable, specified via svid, rather than being explicitly entered into dat.
g02jb uses an iterative process to fit the specified model and for some problems this process may fail to converge (see ${\mathbf{ifail}}={4}$). If the method fails to converge then the maximum number of iterations (see maxit) or tolerance (see tol) may require increasing; try a different starting estimate in gamma. Alternatively, the model can be fit using restricted maximum likelihood (see g02ja) or using the noniterative MIVQUE0.
To fit the model just using MIVQUE0, the first element of gamma should be set to $-1$ and maxit should be set to zero.
Although the quasi-Newton algorithm used in g02jb tends to require more iterations before converging compared to the Newton–Raphson algorithm recommended by Wolfinger et al. (1994), it does not require the second derivatives of the likelihood function to be calculated and consequentially takes significantly less time per iteration.

# Example

The following dataset is taken from Stroup (1989) and arises from a balanced split-plot design with the whole plots arranged in a randomized complete block-design.
In this example the full design matrix for the random independent variable, $Z$, is given by:
 $Z=110000000000000010100000000000001001000000000000000011000000000000001010000000000000100100000000000000001100000000000000101000000000000010010000000000000000110000000000000010100000000000001001110000000000000010100000000000001001000000000000000011000000000000001010000000000000100100000000000000001100000000000000101000000000000010010000000000000000110000000000000010100000000000001001$
 $=A0000A0000A0000AA0000A0000A0000A,$ (1)
where
 $A=110010101001.$
The block structure evident in (1) is modelled by specifying a four-level subject variable, taking the values $\left\{1,1,1,2,2,2,3,3,3,4,4,4,1,1,1,2,2,2,3,3,3,4,4,4\right\}$. The first column of $1\mathrm{s}$ is added to $A$ by setting ${\mathbf{rint}}=1$. The remaining columns of $A$ are specified by a three level factor, taking the values, $\left\{1,2,3,1,2,3,1,\dots \right\}$.

Example program (C#): g02jbe.cs

Example program data: g02jbe.d

Example program results: g02jbe.r