naginterfaces.library.correg.mixeff_​hier_​reml

naginterfaces.library.correg.mixeff_hier_reml(vpr, gamma, comm, iopt=None, ropt=None, io_manager=None)

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

Deprecated since version 27.0.0.0: mixeff_hier_reml is deprecated. Please use lmm_fit() instead. See also the Replacement Calls document.

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

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

Parameters

vprint, array-like, shape $\left(\text{lvpr}\right)$

A vector of flags indicating the mapping between the random variables specified in $\text{rndm}$ and the variance components, ${\sigma }_i^2$. See Further Comments for more details.

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 mixeff_hier_init().

ioptNone or int, array-like, shape $\left(\text{liopt}\right)$, optional

Options passed to the optimization function.

By default mixeff_hier_reml fits the specified model using a modified Newton optimization algorithm as implemented in opt.bounds_mod_deriv2_comp. In some cases, where the calculation of the derivatives is computationally expensive it may be more efficient to use a sequential QP algorithm.

The sequential QP algorithm as implemented in opt.nlp1_solve can be chosen by setting $iopt\left[4\right]=1$.

If $\text{liopt}<5$ or $iopt\left[4\right]\ne 1$, then the modified Newton algorithm will be used.

Different options are available depending on the optimization function used.

In all cases, using a value of $-1$ will cause the default value to be used.

In addition only the first $\text{liopt}$ values of $iopt$ are used, so for example, if only the first element of $iopt$ needs changing and default values for all other options are sufficient $\text{liopt}$ can be set to $1$.

The following table lists the association between elements of $iopt$ and arguments in the optimizer when the modified Newton algorithm is being used.

$i$	Description	Equivalent argument	Default Value
$0$	Number of iterations	$\text{maxcal}$	$1000$
$1$	Unit number for monitoring information	n/a	The advisory unit number
$2$	Print options ($1=$ print)	n/a	$-1$ (no printing performed)
$3$	Frequency that monitoring information is printed	$\text{iprint}$	$-1$
$4$	Optimizer used	n/a	n/a

If requested, monitoring information is displayed in a similar format to that given by the modified Newton optimizer.

The following table lists the association between elements of $iopt$ and options in the optimizer when the sequential QP algorithm is being used.

$i$	Description	Equivalent option	Default Value
$0$	Number of iterations	‘Major Iteration Limit’	$max(50,3\times \text{nvpr})$
$1$	Unit number for monitoring information	n/a	The advisory unit number
$2$	Print options ($1=$ print, otherwise no print)	‘List’/’Nolist’	$-1$ (no printing performed)
$3$	Frequency that monitoring information is printed	‘Major Print Level’	$0$
$4$	Optimizer used	n/a	n/a
$5$	Number of minor iterations	‘Minor Iteration Limit’	$max(50,3\times \text{nvpr})$
$6$	Frequency that additional monitoring information is printed	‘Minor Print Level’	$0$

If $\text{liopt}\le 0$, default values are used for all options and $iopt$ is not referenced.

roptNone or float, array-like, shape $\left(\text{lropt}\right)$, optional

Options passed to the optimization function.

Different options are available depending on the optimization function used.

In all cases, using a value of $-1.0$ will cause the default value to be used.

In addition only the first $\text{lropt}$ values of $ropt$ are used, so for example, if only the first element of $ropt$ needs changing and default values for all other options are sufficient $\text{lropt}$ can be set to $1$.

The following table lists the association between elements of $ropt$ and arguments in the optimizer when the modified Newton algorithm is being used.

$i$	Description	Equivalent argument	Default Value
$0$	Sweep tolerance	n/a	$max(\sqrt{eps},\sqrt{eps}\times {max}_i\left({\text{zz}}_{ii}\right))$
$1$	Lower bound for ${\gamma }^{\ast }$	n/a	$eps/100$
$2$	Upper bound for ${\gamma }^{\ast }$	n/a	${10}^{20}$
$3$	Accuracy of linear minimizations	$\text{eta}$	$0.9$
$4$	Accuracy to which solution is required	$\text{xtol}$	$0.0$
$5$	Initial distance from solution	$\text{stepmx}$	$100000.0$

The following table lists the association between elements of $ropt$ and options in the optimizer when the sequential QP algorithm is being used.

$i$	Description	Equivalent option	Default Value
$0$	Sweep tolerance	n/a	$max(\sqrt{eps},\sqrt{eps}\times {max}_i\left({\text{zz}}_{ii}\right))$
$1$	Lower bound for ${\gamma }^{\ast }$	n/a	$eps/100$
$2$	Upper bound for ${\gamma }^{\ast }$	n/a	${10}^{20}$
$3$	Line search tolerance	‘Line Search Tolerance’	$0.9$
$4$	Optimality tolerance	‘Optimality Tolerance’	${eps}^{0.72}$

where $eps$ is the machine precision returned by machine.precision and ${\text{zz}}_{ii}$ denotes the $i$ diagonal element of $Z^{T}Z$.

If $\text{lropt}\le 0$, then default values are used for all options and $ropt$ may be set to None.

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

effnint

Effective number of observations. If no weights were supplied to mixeff_hier_init() or all supplied weights were 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, then $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$.

estidint, ndarray, shape $(:,\text{lb})$

An array describing the parameter estimates returned in $b$. The first $\text{nlsv}\times \text{nrf}$ columns of $estid$ describe the parameter estimates for the random effects and the last $\text{nff}$ columns the parameter estimates for the fixed effects.

For fixed effects:

for $l=\text{nrf}\times \text{nlsv}+1,\dots ,\text{nrf}\times \text{nlsv}+\text{nff}$

if $b[l-1]$ contains the parameter estimate for the intercept, then

$$estid[0,l-1]=estid[1,l-1]=estid[2,l-1]=0\text{;}$$

if $b[l-1]$ contains the parameter estimate for the $i$th level of the $j$th fixed variable, that is the vector of values held in the $k$th column of $\text{dat}$ when $\text{fixed}[j+1]=k$ then

$$\begin{array}{c}\begin{array}{c}estid[0,l-1]=0\text{,}\\ estid[1,l-1]=j\text{,}\\ estid[2,l-1]=i\text{;}\end{array}\end{array}$$

if the $j$th variable is continuous or binary, that is $\text{levels}\left[\text{fixed}\right[j+1]-1]=1$, $estid[2,l-1]=0$;

any remaining rows of the $l$th column of $estid$ are set to $0$.

For random effects:

let

$N_{R_b}$ denote the number of random variables in the $b$th random statement, that is $N_{R_b}=\text{rndm}(1,b)$;

$R_{jb}$ denote the $j$th random variable from the $b$th random statement, that is the vector of values held in the $k$th column of $\text{dat}$ when $\text{rndm}(2+j,b)=k$;

$N_{S_b}$ denote the number of subject variables in the $b$th random statement, that is $N_{S_b}=\text{rndm}(3+N_{R_b},b)$;

$S_{jb}$ denote the $j$th subject variable from the $b$th random statement, that is the vector of values held in the $k$th column of $\text{dat}$ when $\text{rndm}(3+N_{R_b}+j,b)=k$;

$L\left(S_{jb}\right)$ denote the number of levels for $S_{jb}$, that is $L\left(S_{jb}\right)=\text{levels}[\text{rndm}(3+N_{R_b}+j,b)-1]$;

then

for $l=1,2,\dots \text{nrf}\times \text{nlsv}$, if $b[l-1]$ contains the parameter estimate for the $i$th level of $R_{jb}$ when $S_{\text{k}b}=s_{\text{k}}$, for $\text{k}=1,2,\dots ,N_{S_b}$ and $1\le s_{\text{k}}\le L\left(S_{jb}\right)$, i.e., $s_k$ is a valid value for the $k$th subject variable, then

$$\begin{array}{c}\begin{array}{c}estid[0,l-1]=b\text{,}\\ estid[1,l-1]=j\text{,}\\ estid[2,l-1]=i\text{,}\\ estid[3+k-1,l-1]=s_k\text{,}k=1,2,\dots ,N_{S_b}\text{;}\end{array}\end{array}$$

if the parameter being estimated is for the intercept, then $estid[1,l-1]=estid[2,l-1]=0$;

if the $j$th variable is continuous, or binary, that is $L\left(S_{jb}\right)=1$, then $estid[2,l-1]=0$;

the remaining rows of the $l$th column of $estid$ are set to $0$.

In some situations, certain combinations of variables are never observed.

In such circumstances all elements of the $l$th row of $estid$ are set to $-999$.

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

The parameter estimates, with the first $\text{nrf}\times \text{nlsv}$ elements of $b$ containing the parameter estimates for the random effects, $\nu$, and the remaining $\text{nff}$ elements containing the parameter estimates for the fixed effects, $\beta$. The order of these estimates are described by the $estid$ argument.

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

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

czzfloat, ndarray, shape $(\text{nrf},:)$

If $\text{nlsv}=1$, then $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{nlsv}>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 subject variable, the next $\text{nrf}$ columns the part corresponding to the second level of the overall subject variable etc.

cxxfloat, ndarray, shape $(\text{nff},:)$

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

cxzfloat, ndarray, shape $(\text{nff},:)$

If $\text{nlsv}=1$, then $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 $\text{nlsv}>1$, then $cxz$ 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 subject variable, the next $\text{nrf}$ columns the part corresponding to the second level of the overall subject variable etc.

Raises

NagValueError

(errno $1$)

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

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

(errno $2$)

On entry, $vpr[⟨value⟩]=⟨value⟩$ and $\text{nvpr}=⟨value⟩$.

Constraint: $1\le vpr\left[i\right]\le \text{nvpr}$.

(errno $3$)

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

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

(errno $4$)

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

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

(errno $9$)

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

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

(errno $11$)

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

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

(errno $21$)

On entry, $comm$[‘iopts’] has not been initialized correctly.

(errno $32$)

On entry, at least one value of i, for $\text{i}=1,2,\dots ,\text{nvpr}$, does not appear in $vpr$.

Warns

NagAlgorithmicWarning

(errno $101$)

Optimal solution found, but requested accuracy not achieved.

(errno $102$)

Too many major iterations.

(errno $103$)

Current point cannot be improved upon.

(errno $104$)

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

Notes

mixeff_hier_reml 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$, mixeff_hier_reml 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. mixeff_hier_reml 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).

mixeff_hier_reml 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))$$

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{.}$$

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