naginterfaces.library.correg.mixeff_​reml

naginterfaces.library.correg.mixeff_reml(dat, levels, yvid, cwid, fvid, fint, rvid, nvpr, vpr, rint, svid, gamma, lb, maxit=- 1, tol=0.0, io_manager=None)

mixeff_reml fits a linear mixed effects regression model using restricted maximum likelihood (REML).

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

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

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

Parameters

datfloat, array-like, shape $(n,\text{ncol})$

Array containing all of the data. For the $i$th observation:

$dat[i-1,yvid-1]$ holds the dependent variable, $y$;

if $cwid\ne 0$, $dat[i-1,cwid-1]$ holds the case weights;

if $svid\ne 0$, $dat[i-1,svid-1]$ holds the subject variable.

The remaining columns hold the values of the independent variables.

levelsint, array-like, shape $\left(\text{ncol}\right)$

$levels[i-1]$ 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 $levels[i-1]$ should be $1$; if the variable is discrete then $levels[i-1]$ is the number of levels associated with it and $dat[\text{j}-1,i-1]$ is assumed to take the values $1$ to $levels[i-1]$, for $\text{j}=1,2,\dots ,n$.

yvidint

The column of $dat$ holding the dependent, $y$, variable.

cwidint

The column of $dat$ holding the case weights.

If $cwid=0$, no weights are used.

fvidint, array-like, shape $\left(\text{nfv}\right)$

The columns of the data matrix $dat$ holding the fixed independent variables with $fvid[i-1]$ holding the column number corresponding to the $i$th fixed variable.

fintint

Flag indicating whether a fixed intercept is included ($fint=1$).

rvidint, array-like, shape $\left(\text{nrv}\right)$

The columns of the data matrix $dat$ holding the random independent variables with $rvid[i-1]$ holding the column number corresponding to the $i$th random variable.

nvprint

If $rint=1$ and $svid\ne 0$, $nvpr$ is the number of variance components being $\text{estimated}-2$, ($g-1$), else $nvpr=g$.

If $\text{nrv}=0$, $nvpr$ is not referenced.

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

$vpr[i-1]$ 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=vpr[i-1]+1$ if $rint=1$ and $svid\ne 0$ and $j=vpr[i-1]$ otherwise. Random variables with the same value of $j$ are assumed to be taken from the same distribution.

rintint

Flag indicating whether a random intercept is included ($rint=1$).

If $svid=0$, $rint$ is not referenced.

svidint

The column of $dat$ holding the subject variable.

If $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\times Z_S$ denote the interaction between variables $Z_1$ and $Z_S$, fitting a model with $rint=0$, random-effects $Z_1+Z_2$ and subject variable $Z_S$ is equivalent to fitting a model with random-effects $Z_1\times Z_S+Z_2\times Z_S$ and no subject variable.

If $rint=1$ the model is equivalent to fitting $Z_S+Z_1\times Z_S+Z_2\times Z_S$ and no subject variable.

gammafloat, array-like, shape $(nvpr+2)$

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 ,g$. If $rint=1$ and $svid\ne 0$, $g=nvpr+1$, else $g=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.

lbint

The size of the array $b$.

maxitint, optional

The maximum number of iterations.

If $maxit<0$, the default value of $100$ is used.

If $maxit=0$, the parameter estimates $(\beta ,\nu )$ and corresponding standard errors are calculated based on the value of ${\gamma }_0$ supplied in $gamma$.

tolfloat, optional

The tolerance used to assess convergence.

If $tol\le 0.0$, the default value of ${ϵ}^{0.7}$ is used, where $ϵ$ is the machine precision.

io_managerFileObjManager, optional

Manager for I/O in this routine.

Returns

gammafloat, ndarray, shape $(nvpr+2)$

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

nffint

The number of fixed effects estimated (i.e., the number of columns, $p$, in the design matrix $X$).

nrfint

The number of random effects estimated (i.e., the number of columns, $q$, in the design matrix $Z$).

dfint

The degrees of freedom.

remlfloat

$-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 $\left(lb\right)$

The parameter estimates, $(\beta ,\nu )$, 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 $fint=1$, $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=levels\left[fvid\right[i-1]-1]$.

Define

if $fint=1$, $F_1=2$ else if $fint=0$, $F_1=1$;

$F_{i+1}=F_i+max(L_i-1,1)$, $i\ge 1$.

Then for $i=1,2,\dots ,\text{nfv}$:

if $L_i>1$, $b[F_i+\text{j}-2-1]$ contains the parameter estimate for the $\text{j}$th level of the $i$th fixed variable, for $\text{j}=2,3,\dots ,L_i$;

if $L_i\le 1$, $b[F_i-1]$ 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=levels\left[rvid\right[i-1]-1]$.

Define

if $rint=1$, $R_1=2$ else if $rint=0$, $R_1=1$;

$R_{i+1}=R_i+L_i$, $i\ge 1$.

Then for $i=1,2,\dots ,\text{nrv}$:

if $svid=0$,

if $L_i>1$, $b[nff+R_i+\text{j}-1-1]$ contains the parameter estimate for the $\text{j}$th level of the $i$th random variable, for $\text{j}=1,2,\dots ,L_i$;

if $L_i\le 1$, $b[nff+R_i-1]$ contains the parameter estimate for the $i$th random variable;

if $svid\ne 0$,

let $L_S$ denote the number of levels associated with the subject variable, that is $L_S=levels[svid-1]$;

if $L_i>1$, $b[nff+(\text{s}-1)L_S+R_i+\text{j}-1-1]$ contains the parameter estimate for the interaction between the $\text{s}$th level of the subject variable and the $\text{j}$th level of the $i$th random variable, for $\text{j}=1,2,\dots ,L_i$, for $\text{s}=1,2,\dots ,L_S$;

if $L_i\le 1$, $b[nff+(\text{s}-1)L_S+R_i-1]$ contains the parameter estimate for the interaction between the $\text{s}$th level of the subject variable and the $i$th random variable, for $\text{s}=1,2,\dots ,L_S$;

if $rint=1$, $b\left[nff\right]$ contains the estimate of the random intercept.

sefloat, ndarray, shape $\left(lb\right)$

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

warnint

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 $warn=1$, the negative estimate is set to zero and the estimation process allowed to continue.

Raises

NagValueError

(errno $1$)

On entry, $n=⟨value⟩$.

Constraint: $n\ge 1$.

(errno $1$)

On entry, $n=⟨value⟩$.

Constraint: number of observations with nonzero weights must be greater than one.

(errno $1$)

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

Constraint: $\text{ncol}\ge 1$.

(errno $1$)

On entry, $\text{nfv}=⟨value⟩$ and $\text{ncol}=⟨value⟩$.

Constraint: $0\le \text{nfv}<\text{ncol}$.

(errno $1$)

On entry, $\text{nrv}=⟨value⟩$ and $\text{ncol}=⟨value⟩$.

Constraint: $0\le \text{nrv}<\text{ncol}$ and $\text{nrv}+rint>0$.

(errno $1$)

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

Constraint: $0\le nvpr\le \text{nrv}$ and ($\text{nrv}\ne 0$ or $nvpr\ge 1$).

(errno $1$)

On entry, $yvid=⟨value⟩$ and $\text{ncol}=⟨value⟩$.

Constraint: $1\le yvid\le \text{ncol}$.

(errno $1$)

On entry, $svid=⟨value⟩$ and $\text{ncol}=⟨value⟩$.

Constraint: $0\le svid\le \text{ncol}$.

(errno $1$)

On entry, $cwid=⟨value⟩$ and $\text{ncol}=⟨value⟩$.

Constraint: $0\le cwid\le \text{ncol}$ and any supplied weights must be $\ge 0.0$.

(errno $1$)

On entry, $fint=⟨value⟩$.

Constraint: $fint=0$ or $1$.

(errno $1$)

On entry, $rint=⟨value⟩$.

Constraint: $rint=0$ or $1$.

(errno $1$)

On entry, $lb$ too small: $lb=⟨value⟩$.

(errno $2$)

On entry, $levels\left[\text{i}\right]<1$, for at least one $\text{i}$.

(errno $2$)

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

Constraint: $1\le fvid\left[i\right]\le \text{ncol}$, for all $i$.

(errno $2$)

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

Constraint: $1\le rvid\left[i\right]\le \text{ncol}$, for all $i$.

(errno $2$)

On entry, $nvpr=⟨value⟩$.

Constraint: $1\le vpr\left[i\right]\le nvpr$, for all $i$.

(errno $2$)

On entry, invalid data: categorical variable with value greater than that specified in $levels$.

(errno $2$)

On entry, $gamma\left[i\right]<0.0$, for at least one $i$.

(errno $3$)

Degrees of freedom $<1$: $df=⟨value⟩$.

(errno $4$)

Routine failed to converge to specified tolerance: $tol=⟨value⟩$.

(errno $4$)

Routine failed to converge in $maxit$ iterations: $maxit=⟨value⟩$.

Notes

mixeff_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

$$\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 expectations 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_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_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_reml fits the model using a quasi-Newton algorithm to maximize the restricted log-likelihood function:

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

where

$$V=ZG{Z}^{\prime }+R\text{,}r=y-Xb\text{and}b={\left(X^{\prime }{V}^{-1}X\right)}^{-1}{X}^{\prime }{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^{\prime }{V}^{-1}r\right)/(n-p)\text{.}$$

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