naginterfaces.library.correg.lmm_​fit

naginterfaces.library.correg.lmm_fit(hlmm, gamma, comm, wantc=False, io_manager=None)

lmm_fit fits a multi-level linear mixed effects regression model using restricted maximum likelihood (REML) or maximum likelihood (ML). Prior to calling lmm_fit the initialization function lmm_init() must be called.

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

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

Parameters

hlmmHandle

A G22 handle to the internal data structure containing a description of the required model as returned in $\text{hlmm}$ by lmm_init().

gammafloat, array-like, shape $(\text{nvpr}+1)$

Holds the initial values of the variance components, ${\gamma }_0$, with $gamma[\text{i}-1]$ the initial value for ${\sigma }_{\text{i}}^2/{\sigma }_R^2$, for $\text{i}=1,2,\dots ,\text{nvpr}$.

If $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.

commdict, communication object, modified in place

Communication structure.

This argument must have been initialized by a prior call to lmm_init().

wantcbool, optional

Flag indicating whether the arrays $czz$, $cxx$ and $cxz$ are required.

io_managerFileObjManager, optional

Manager for I/O in this routine.

Returns

gammafloat, ndarray, shape $(\text{nvpr}+1)$

$gamma[\text{i}-1]$, for $\text{i}=1,2,\dots ,\text{nvpr}$, holds the final estimate of ${\sigma }_{\text{i}}^2$ and $gamma\left[\text{nvpr}\right]$ holds the final estimate for ${\sigma }_R^2$.

Labels for the variance components can be obtained using blgm.lm_submodel.

effnint

Effective number of observations. If there are no weights, or all weights are nonzero, $effn=\text{n}$.

rnkxint

The rank of the design matrix, $X$, for the fixed effects.

ncovint

Number of variance components not estimated to be zero. If none of the variance components are estimated to be zero, $ncov=\text{nvpr}$.

lnlikefloat

$-2l_R\left(\hat{\gamma }\right)$ where $l_R$ is the log of the restricted maximum likelihood calculated at $\hat{\gamma }$, the estimated variance components returned in $gamma$.

bfloat, ndarray, shape $(\text{fnlsv}\times \text{nff}+\text{rnlsv}\times \text{nrf})$

The parameter estimates, with the first $\text{nrf}\times \text{rnlsv}$ elements of $b$ containing the parameter estimates for the random effects, $\nu$, and the remaining $\text{nff}\times \text{fnlsv}$ elements containing the parameter estimates for the fixed effects, $\beta$.

Labels for the parameter estimates can be obtained using blgm.lm_submodel.

sefloat, ndarray, shape $(\text{fnlsv}\times \text{nff}+\text{rnlsv}\times \text{nrf})$

The standard errors of the parameter estimates given in $b$.

czzNone or float, ndarray, shape $(:,:)$

If $\text{rnlsv}=1$, $czz$ holds the lower triangular portion of the matrix $\left(1/{\sigma }^2\right)(Z^T{\hat{R}}^{-1}Z+{\hat{G}}^{-1})$, where $\hat{R}$ and $\hat{G}$ are the estimates of $R$ and $G$ respectively.

If $\text{rnlsv}>1$, then $czz$ holds this matrix in compressed form, with the first $\text{nrf}$ columns holding the part of the matrix corresponding to the first level of the overall random subject variable, the next $\text{nrf}$ columns holding the part of the matrix corresponding to the second level of the overall random subject variable etc.

If $wantc=False$, $czz$ is returned as None.

cxxNone or float, ndarray, shape $(:,:)$

If $\text{fnlsv}=1$, $cxx$ holds the lower triangular portion of the matrix $\left(1/{\sigma }^2\right)\left(X^T{\hat{V}}^{-1}X\right)$, where $\hat{V}$ is the estimated value of $V$.

If $\text{fnlsv}>1$, then $cxx$ holds this matrix in compressed form, with the first $\text{nff}$ columns holding the part of the matrix corresponding to the first level of the overall fixed subject variable, the next $\text{nff}$ columns holding the part of the matrix corresponding to the second level of the overall fixed subject variable, etc.

If $wantc=False$, $cxx$ is returned as None.

cxzNone or float, ndarray, shape $(:,:)$

$cxz$ holds the matrix $\left(1/{\sigma }^2\right)\left(X^T{\hat{V}}^{-1}Z\right)\hat{G}$, where $\hat{V}$ and $\hat{G}$ are the estimates of $V$ and $G$ respectively.

If $wantc=False$, $cxz$ is returned as None.

Raises

NagValueError

(errno $11$)

$hlmm$ has not been initialized or is corrupt.

(errno $12$)

$hlmm$ is not a G22 handle as generated by lmm_init().

(errno $21$)

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

Constraint: $\text{nvpr}\ge ⟨value⟩$.

(errno $31$)

On entry, $gamma[⟨value⟩]=⟨value⟩$.

Constraint: $gamma\left[0\right]=-1.0$ or $gamma[i-1]\ge 0.0$.

(errno $81$)

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

Constraint: $\text{lb}\ge ⟨value⟩$.

(errno $171$)

On entry, the communication arrays, $comm$[‘iopts’] and $comm$[‘opts’], have not been initialized correctly.

Warns

NagAlgorithmicWarning

(errno $1001$)

Optimal solution found, but requested accuracy not achieved.

(errno $1002$)

Too many major iterations.

(errno $1003$)

Current point cannot be improved upon.

(errno $1004$)

At least one negative estimate for $gamma$ was obtained. All negative estimates have been set to zero.

Notes

lmm_fit fits a model of the form:

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

where	$y$ is a vector of $n$ observations on the dependent variable,
	$X$ is a known $n\times p$ design matrix for the fixed independent variables,
	$\beta$ is a vector of length $p$ of unknown fixed effects,
	$Z$ is a known $n\times q$ design matrix for the random independent variables,
	$\nu$ is a vector of length $q$ of unknown random effects,
and	$ϵ$ 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\}$. Rather than working directly with $\gamma$, lmm_fit uses an iterative process to estimate ${\gamma }^{\ast }=\{{\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\}$. Due to the iterative nature of the estimation a set of initial values, ${\gamma }_0$, for ${\gamma }^{\ast }$ is required. lmm_fit 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).

lmm_fit fits the model by maximizing the restricted log-likelihood function:

$$-2l_R=\log \left(\left|V\right|\right)+(n-p)\log \left(r^T{V}^{-1}r\right)+\log \left(\left|X^T{V}^{-1}X\right|\right)+(n-p)(1+\log \left(2\pi /(n-p)\right))$$

or the log-likelihood function:

$$-2l_R=\log \left(\left|V\right|\right)+n\log \left(r^T{V}^{-1}r\right)+\log \left(2\pi /n\right)$$

where

$$V=ZG{Z}^T+R\text{,}r=y-Xb\text{and}b={\left(X^T{V}^{-1}X\right)}^{-1}{X}^T{V}^{-1}y\text{.}$$

By default the restricted log-likelihood function is used, the log-likelihood function can be chosen through the option ‘Solver’ as detailed in the documentation for lmm_init().

Once the final estimates for ${\gamma }^{\ast }$ have been obtained, the value of ${\sigma }_R^2$ is given by

$${\sigma }_R^2=\left(r^T{V}^{-1}r\right)/(n-p)\text{.}$$

Case weights, $W_c$, can be incorporated into the model by replacing $X^{T}X$ and $Z^{T}Z$ with $X^T{W}_{c}X$ and $Z^T{W}_{c}Z$ respectively, for a diagonal weight matrix $W_c$.

The log-likelihood, $l_R$, is calculated using the sweep algorithm detailed in Wolfinger et al. (1994).

References

Goodnight, J H, 1979, A tutorial on the SWEEP operator, The American Statistician (33(3)), 149–158

Harville, D A, 1977, Maximum likelihood approaches to variance component estimation and to related problems, JASA (72), 320–340

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

Stroup, W W, 1989, Predictable functions and prediction space in the mixed model procedure, Applications of Mixed Models in Agriculture and Related Disciplines (Southern Cooperative Series Bulletin No. 343), 39–48

Wolfinger, R, Tobias, R and Sall, J, 1994, Computing Gaussian likelihoods and their derivatives for general linear mixed models, SIAM Sci. Statist. Comput. (15), 1294–1310

See also

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