naginterfaces.library.correg.lmm_​init

naginterfaces.library.correg.lmm_init(hlmm, hddesc, hfixed, y, dat, hrndm=None, wt=None)

lmm_init preprocesses a dataset prior to fitting a linear mixed effects regression model via lmm_fit().

Note: this function uses optional algorithmic parameters, see also: blgm.optset, blgm.optget, lmm_fit().

For full information please refer to the NAG Library document for g02jf

https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/g02/g02jff.html

Parameters

hlmmHandle, modified in place

On entry: must be set to a null Handle or, alternatively, an existing G22 handle may be supplied in which case lmm_init will destroy the supplied G22 handle as if blgm.handle_free had been called.

On exit: holds a G22 handle to the internal data structure containing a description of the model. You must not change the G22 handle other than through the functions in submodule correg or submodule blgm.

hddescHandle

A G22 handle to the internal data structure containing a description of the data matrix, $D$, as returned in $\text{hddesc}$ by blgm.lm_describe_data.

hfixedHandle

A G22 handle to the internal data structure containing a description of the fixed part of the model $M_f$ as returned in $\text{hform}$ by blgm.lm_formula.

If $hfixed$ is a null Handle then the model is assumed to not have a fixed part.

yfloat, array-like, shape $\left(n\right)$

$y$, the vector of observations on the dependent variable.

datfloat, array-like, shape $(:,:)$

The data matrix, $D$. By default, $D_{ij}$, the $\text{i}$th value for the $\text{j}$th variable, for $\text{j}=1,2,\dots ,m_d$, for $\text{i}=1,2,\dots ,n$, should be supplied in $dat[i-1,j-1]$.

If the option ‘Storage Order’, described in blgm.lm_describe_data, is set to ‘VAROBS’, $D_{ij}$ should be supplied in $dat[j-1,i-1]$.

If either $y_i$, $w_i$ or $D_{ij}$, for a variable $j$ used in the model, is NaN (Not A Number) then that value is treated as missing and the whole observation is excluded from the analysis.

hrndmNone or Handle, list, shape $\left(\text{nrndm}\right)$, optional

A series of G22 handles to internal data structures containing a description of the random part of the model $M_r$ as returned in $\text{hform}$ by blgm.lm_formula.

wtNone or float, array-like, shape $\left(n\right)$, optional

Optionally, the diagonal elements of the weight matrix $W_c$.

If $wt[i-1]=0.0$, the $i$th observation is not included in the model and the effective number of observations is the number of observations with nonzero weights.

If weights are not provided then $wt$ must be set to None, and the effective number of observations is $n$.

Returns

fnlsvint

The number of levels for the overall subject variable in $M_f$. If there is no overall subject variable, $fnlsv=1$.

nffint

The number of fixed effects estimated in each of the $fnlsv$ subject blocks. The number of columns, $p$, in the design matrix $X$ is given by $p=nff\times fnlsv$.

rnlsvint

The number of levels for the overall subject variable in $M_r$. If there is no overall subject variable, $rnlsv=1$.

nrfint

The number of random effects estimated in each of the $rnlsv$ subject blocks. The number of columns, $q$, in the design matrix $Z$ is given by $q=nrf\times rnlsv$.

nvprint

$g$, the number of variance components being estimated (excluding the overall variance, ${\sigma }_R^2$). This is defined by the number of terms in the random part of the model, $M_r$ (see Algorithmic Details for details).

commdict, communication object

Communication structure.

Other Parameters

‘Gamma Lower Bound’float

Default $=\sqrt{\text{machine precision}}/100$

A lower bound for the elements of ${\gamma }^{\ast }$, where ${\gamma }^{\ast }=\gamma /{\sigma }_R^2$.

‘Gamma Upper Bound’float

Default $={10}^{20}$

An upper bound for the elements of ${\gamma }^{\ast }$, where ${\gamma }^{\ast }=\gamma /{\sigma }_R^2$.

‘Initial Distance’float

Default $=100000.0$

The initial distance from the solution.

When $\text{‘Solver'}=\text{'E04LB'}$, lmm_fit() passes ‘Initial Distance’ to the solver as $\text{stepmx}$.

When $\text{‘Solver'}=\text{'E04UC'}$, this option is ignored.

‘Initial Value Strategy’int

Default $=\text{special}$

Controls how lmm_fit() will choose the initial values for the variance components, $\gamma$, if not supplied.

$\text{‘Initial Value Strategy'}=0$

The MIVQUE0 estimates of the variance components based on the likelihood specified by ‘Likelihood’ are used.

$\text{‘Initial Value Strategy'}=1$

The MIVQUE0 estimates based on the maximum likelihood are used, irrespective of the value of ‘Likelihood’.

See Rao (1972) for a description of the minimum variance quadratic unbiased estimators (MIVQUE0).

By default, for small problems, $\text{‘Initial Value Strategy'}=0$ and for large problems $\text{‘Initial Value Strategy'}=1$.

‘Likelihood’str

Default $=\text{'REML'}$

‘Likelihood’ defines whether lmm_fit() will use the restricted maximum likelihood (REML) or the maximum likelihood (ML) when fitting the model.

‘Linear Minimization Accuracy’float

Default $=0.9$

The accuracy of the linear minimizations.

When $\text{‘Solver'}=\text{'E04LB'}$, lmm_fit() passes ‘Linear Minimization Accuracy’ to the solver as $\text{eta}$.

When $\text{‘Solver'}=\text{'E04UC'}$, this option is ignored.

‘Line Search Tolerance’float

Default $=0.9$

The line search tolerance.

When $\text{‘Solver'}=\text{'E04LB'}$, this option is ignored.

When $\text{‘Solver'}=\text{'E04UC'}$, lmm_fit() passes ‘Line Search Tolerance’ to the solver as ‘Line Search Tolerance’.

‘List’valueless

Option ‘List’ enables printing of each option specification as it is supplied. ‘NoList’ suppresses this printing.

‘NoList’valueless

Default

Option ‘List’ enables printing of each option specification as it is supplied. ‘NoList’ suppresses this printing.

‘Major Iteration Limit’int

Default $=\text{special}$

The number of major iterations.

When $\text{‘Solver'}=\text{'E04LB'}$, lmm_fit() passes ‘Major Iteration Limit’ to the solver as $\text{maxcal}$. In this case, the default value used is $1000$.

When $\text{‘Solver'}=\text{'E04UC'}$, lmm_fit() passes ‘Major Iteration Limit’ to the solver as ‘Major Iteration Limit’. In this case, the default value used is $max(50,3\times g)$, where $g$ is the number of variance components being estimated (excluding the overall variance, ${\sigma }_R^2$).

‘Major Print Level’int

Default $=\text{special}$

The frequency that monitoring information is output to ‘Unit Number’.

When $\text{‘Solver'}=\text{'E04LB'}$, lmm_fit() passes ‘Major Print Level’ to the solver as $\text{iprint}$. In this case, the default value used is $-1$ and hence no monitoring information will be output.

When $\text{‘Solver'}=\text{'E04UC'}$, lmm_fit() passes ‘Major Print Level’ to the solver as ‘Major Print Level’. In this case, the default value used is $0$ and hence no monitoring information will be output.

‘Maximum Number of Threads’int

Default $=\text{special}$

Controls the maximum number of threads used by lmm_fit() in a multithreaded library. By default, the maximum number of available threads are used.

In a library that is not multithreaded, this option has no effect.

‘Minor Iteration Limit’int

Default $=max(50,3\times g)$

The number of minor iterations.

When $\text{‘Solver'}=\text{'E04LB'}$, this option is ignored.

When $\text{‘Solver'}=\text{'E04UC'}$, lmm_fit() passes ‘Minor Iteration Limit’ to the solver as ‘Minor Iteration Limit’. In this case, the default value used is $max(50,3\times g)$, where $g$ is the number of variance components being estimated (excluding the overall variance, ${\sigma }_R^2$).

‘Minor Print Level’int

Default $=0$

The frequency that additional monitoring information is output to ‘Unit Number’.

When $\text{‘Solver'}=\text{'E04LB'}$, this option is ignored.

When $\text{‘Solver'}=\text{'E04UC'}$, lmm_fit() passes ‘Minor Print Level’ to the solver as ‘Minor Print Level’. The default value of $0$ means that no additional monitoring information will be output.

‘Optimality Tolerance’float

Default $={\text{machine precision}}^{0.72}$

The optimality tolerance.

When $\text{‘Solver'}=\text{'E04LB'}$, this option is ignored.

When $\text{‘Solver'}=\text{'E04UC'}$, lmm_fit() passes ‘Optimality Tolerance’ to the solver as ‘Optimality Tolerance’.

‘Parallelisation Strategy’int

Default $=\text{special}$

If $\text{‘Maximum Number of Threads'}>0$ then ‘Parallelisation Strategy’ controls how lmm_fit() is parallelised in a multithreaded library.

$\text{‘Parallelisation Strategy'}=1$

lmm_fit() will attempt to parallelise operations involving $Z$, even if $rnlsv=1$.

$\text{‘Parallelisation Strategy'}=2$

lmm_fit() will only attempt to parallelise operations involving $Z$, if $rnlsv>1$.

By default, $\text{‘Parallelisation Strategy'}=1$, however, for some models / datasets, this may be slower than using $\text{‘Parallelisation Strategy'}=2$ when $rnlsv=1$.

In a library that is not multithreaded, this option has no effect.

‘Solution Accuracy’float

Default $=0.0$

The accuracy to which the solution is required.

When $\text{‘Solver'}=\text{'E04LB'}$, lmm_fit() passes ‘Solution Accuracy’ to the solver as $\text{xtol}$.

When $\text{‘Solver'}=\text{'E04UC'}$, this option is ignored.

‘Solver’str

Default $=\text{special}$

Controls which solver lmm_fit() will use when fitting the model. By default, $\text{‘Solver'}=\text{'E04LB'}$ is used for small problems and $\text{‘Solver'}=\text{'E04UC'}$, otherwise.

If $\text{‘Solver'}=\text{'E04LB'}$, then the solver used is the one implemented in opt.bounds_mod_deriv2_comp and if $\text{‘Solver'}=\text{'E04UC'}$, then the solver used is the one implemented in opt.nlp1_solve.

‘Sweep Tolerance’float

Default $=\text{special}$

The sweep tolerance used by lmm_fit() when performing the sweep operation Wolfinger et al. (1994). The default value used is $\text{‘Sweep Tolerance'}=max(ϵ,ϵ\times \left({\text{max}}_i\left({\left(Z^T\right)}_{ii}\right)\right))$, where $ϵ=\sqrt{\text{machine precision}}$.

‘Unit Number’int

Default $=\text{advisory message unit number}$

The monitoring unit number to which lmm_fit() will send any monitoring information.

Raises

NagValueError

(errno $11$)

On entry, $hlmm$ is not a null Handle or a recognised G22 handle.

(errno $21$)

$hddesc$ has not been initialized or is corrupt.

(errno $22$)

$hddesc$ is not a G22 handle as generated by blgm.lm_describe_data.

(errno $31$)

$hfixed$ has not been initialized or is corrupt.

(errno $32$)

$hfixed$ is not a G22 handle as generated by blgm.lm_formula.

(errno $33$)

A variable name used when creating $hfixed$ is not present in $hddesc$.

Variable name: $⟨value⟩$.

(errno $41$)

On entry, $\text{nrndm}=⟨value⟩$.

Constraint: $\text{nrndm}\ge 0$.

(errno $51$)

$i=⟨value⟩$.

$hrndm[i-1]$ has not been initialized or is corrupt.

(errno $52$)

$i=⟨value⟩$.

$hrndm[i-1]$ is not a G22 handle as generated by blgm.lm_formula.

(errno $53$)

No model has been specified.

(errno $54$)

A variable name used when creating $hrndm$ is not present in $hddesc$.

Variable name: $⟨value⟩$.

(errno $61$)

On entry, $\text{lwt}=⟨value⟩$ and $n=⟨value⟩$.

Constraint: $\text{lwt}=0$ or $n$.

(errno $71$)

On entry, $n=⟨value⟩$.

Constraint: $n\ge 1$.

(errno $72$)

On entry, $n=⟨value⟩$ and $n_d=⟨value⟩$.

Constraint: $n\le n_d$.

(errno $73$)

On entry, no observations due to zero weights or missing values.

(errno $91$)

On entry, $i=⟨value⟩$ and $wt[i-1]=⟨value⟩$.

Constraint: $wt[i-1]\ge 0.0$.

(errno $101$)

On entry, column $j$ of the data matrix, $D$, is not consistent with information supplied in $hddesc$, $j=⟨value⟩$.

(errno $111$)

On entry, $n=⟨value⟩$ and $\text{lddat}=⟨value⟩$.

Constraint: $\text{lddat}\ge n$.

(errno $112$)

On entry, $m_d=⟨value⟩$ and $\text{lddat}=⟨value⟩$.

Constraint: $\text{lddat}\ge m_d$.

(errno $121$)

On entry, $m_d=⟨value⟩$ and $\text{sddat}=⟨value⟩$.

Constraint: $\text{sddat}\ge m_d$.

(errno $122$)

On entry, $n=⟨value⟩$ and $\text{sddat}=⟨value⟩$.

Constraint: $\text{sddat}\ge n$.

Warns

NagAlgorithmicWarning

(errno $34$)

The fixed part of the model contains categorical variables, but no intercept or main effects terms have been requested.

(errno $102$)

Column $j$ of the data matrix, $D$, required rounding more than expected when being treated as a categorical variable, $j=⟨value⟩$.

Notes

lmm_init must be called prior to fitting a linear mixed effects regression model via lmm_fit().

The model is of the form:

$$y=X\beta +Z\nu +ϵ$$

where	$y$ is a vector of $n$ observations on the dependent variable,
	$X$ is an $n\times p$ design matrix of fixed independent variables,
	$\beta$ is a vector of $p$ unknown fixed effects,
	$Z$ is an $n\times q$ design matrix of random independent variables,
	$\nu$ is a vector of length $q$ of unknown random effects,
	$ϵ$ is a vector of length $n$ of unknown random errors.

Both $\nu$ and $ϵ$ are assumed to have a Gaussian distribution with expectation zero and variance/covariance matrix defined by

$$\begin{array}{c}Var\left[\begin{array}{c}\nu \\ ϵ\end{array}\right]=\left[\begin{array}{cc}G& 0\\ 0& R\end{array}\right]\end{array}$$

where $R={\sigma }_R^{2}I$, $I$ is the $n\times n$ identity matrix and $G$ is a diagonal matrix. It is assumed that the random variables, $Z$, can be subdivided into $g\le q$ groups with each group being identically distributed with expectation zero and variance ${\sigma }_i^2$. The diagonal elements of matrix $G$, therefore, take one of the values $\{{\sigma }_i^2:i=1,2,\dots ,g\}$, depending on which group the associated random variable belongs to.

The model, therefore, contains three sets of unknowns: the fixed effects $\beta$, the random effects $\nu$ and a vector of $g+1$ variance components $\gamma$, where $\gamma =\{{\sigma }_1^2,{\sigma }_2^2,\dots ,{\sigma }_{g-1}^2,{\sigma }_g^2,{\sigma }_R^2\}$.

Case weights can be incorporated into the model by replacing $X$ and $Z$ with $W_c^{1/2}X$ and $W_c^{1/2}Z$ respectively where $W_c$ is a diagonal weight matrix.

The design matrices, $X$ and $Z$, are constructed from an $n\times m_d$ data matrix, $D$, a description of the fixed independent variables, $M_f$, and a description of the random independent variables, $M_r$. See Algorithmic Details for further details.

References

Rao, C R, 1972, Estimation of variance and covariance components in a linear model, J. Am. Stat. Assoc. (67), 112–115

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

See also

naginterfaces.library.examples.correg.lmm_init_combine_ex.main()