# NAG CL Interfaceg02jhc (lmm_​fit)

Settings help

CL Name Style:

## 1Purpose

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

## 2Specification

 #include
 void g02jhc (void *hlmm, Integer nvpr, double gamma[], Integer *effn, Integer *rnkx, Integer *ncov, double *lnlike, Integer lb, double b[], double se[], double czz[], Integer pdczz, double cxx[], Integer pdcxx, double cxz[], Integer pdcxz, double rcomm[], Integer icomm[], NagError *fail)
The function may be called by the names: g02jhc or nag_correg_lmm_fit.

## 3Description

g02jhc 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×p$ design matrix for the fixed independent variables, $\beta$ is a vector of length $p$ of unknown fixed effects, $Z$ is a known $n×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 variance/covariance matrix defined by
 $Var[ ν ε ] = [ G 0 0 R ]$
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 expectation 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$, g02jhc 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. g02jhc 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).
g02jhc fits the model by maximizing the restricted log-likelihood function:
 $−2 l R = log(|V|) + (n-p) log( r TV-1r) + log| X TV-1X| + (n-p) (1+log(2π/(n-p)))$
or the log-likelihood function:
 $−2 l R = log(|V|) + n log( r TV-1r) + log(2π/n)$
where
 $V = ZG ZT + R, r=y-Xb and b = (XTV-1X)-1 XT V-1 y .$
By default the restricted log-likelihood function is used, the log-likelihood function can be chosen through the optional parameter ${\mathbf{Solver}}$ as detailed in the documentation for g02jfc.
Once the final estimates for ${\gamma }^{*}$ have been obtained, the value of ${\sigma }_{R}^{2}$ is given by
 $σR2 = (rTV-1r) / (n-p) .$
Case weights, ${W}_{c}$, can be incorporated into the model by replacing ${X}^{\mathrm{T}}X$ and ${Z}^{\mathrm{T}}Z$ with ${X}^{\mathrm{T}}{W}_{c}X$ and ${Z}^{\mathrm{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).

## 4References

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

## 5Arguments

Note: prior to calling g02jhc the initialization function g02jfc must be called, therefore, this documentation should be read in conjunction with the document for g02jfc. In particular some argument names and conventions described in that document are also relevant here, but their definition has not been repeated. Specifically, hlmm should be interpreted identically in both functions.
1: $\mathbf{hlmm}$void * Input
On entry: a G22 handle to the internal data structure containing a description of the required model as returned in hlmm by g02jfc.
2: $\mathbf{nvpr}$Integer Input
On entry: $g$, the number of variance components being estimated (excluding the overall variance, ${\sigma }_{R}^{2}$).
This should be set to the value of nvpr returned by g02jfc.
3: $\mathbf{gamma}\left[{\mathbf{nvpr}}+1\right]$double Input/Output
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 ,{\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 ,{\mathbf{nvpr}}$, holds the final estimate of ${\sigma }_{\mathit{i}}^{2}$ and ${\mathbf{gamma}}\left[{\mathbf{nvpr}}\right]$ holds the final estimate for ${\sigma }_{R}^{2}$.
Labels for the variance components can be obtained using g22ydc.
Constraint: ${\mathbf{gamma}}\left[0\right]=-1.0$ or ${\mathbf{gamma}}\left[\mathit{i}-1\right]\ge 0.0$, for $\mathit{i}=1,2,\dots ,g$.
4: $\mathbf{effn}$Integer * Output
On exit: effective number of observations. If there are no weights (i.e., ${\mathbf{wt}}\phantom{\rule{0.25em}{0ex}}\text{is}\phantom{\rule{0.25em}{0ex}}\mathbf{NULL}$), or all weights are nonzero, ${\mathbf{effn}}={\mathbf{n}}$.
5: $\mathbf{rnkx}$Integer * Output
On exit: the rank of the design matrix, $X$, for the fixed effects.
6: $\mathbf{ncov}$Integer * Output
On exit: number of variance components not estimated to be zero. If none of the variance components are estimated to be zero, ${\mathbf{ncov}}={\mathbf{nvpr}}$.
7: $\mathbf{lnlike}$double * Output
On exit: $-2{l}_{R}\left(\stackrel{^}{\gamma }\right)$ where ${l}_{R}$ is the log of the restricted maximum likelihood calculated at $\stackrel{^}{\gamma }$, the estimated variance components returned in gamma.
8: $\mathbf{lb}$Integer Input
On entry: the dimension of the arrays b and se.
Constraint: ${\mathbf{lb}}\ge {\mathbf{fnlsv}}×{\mathbf{nff}}+{\mathbf{rnlsv}}×{\mathbf{nrf}}$.
9: $\mathbf{b}\left[{\mathbf{lb}}\right]$double Output
On exit: the parameter estimates, with the first ${\mathbf{nrf}}×{\mathbf{rnlsv}}$ elements of ${\mathbf{b}}$ containing the parameter estimates for the random effects, $\nu$, and the remaining ${\mathbf{nff}}×{\mathbf{fnlsv}}$ elements containing the parameter estimates for the fixed effects, $\beta$.
Labels for the parameter estimates can be obtained using g22ydc.
10: $\mathbf{se}\left[{\mathbf{lb}}\right]$double Output
On exit: the standard errors of the parameter estimates given in b.
11: $\mathbf{czz}\left[{\mathbf{pdczz}}×{\mathbf{rnlsv}}×{\mathbf{nrf}}\right]$double Output
Note: where ${\mathbf{CZZ}}\left(i,j\right)$ appears in this document, it refers to the array element ${\mathbf{czz}}\left[\left(j-1\right)×{\mathbf{pdczz}}+i-1\right]$.
On exit: if ${\mathbf{rnlsv}}=1$, CZZ holds the lower triangular portion of the matrix $\left(1/{\sigma }^{2}\right)\left({Z}^{\mathrm{T}}{\stackrel{^}{R}}^{-1}Z+{\stackrel{^}{G}}^{-1}\right)$, where $\stackrel{^}{R}$ and $\stackrel{^}{G}$ are the estimates of $R$ and $G$ respectively.
If ${\mathbf{rnlsv}}>1$, then CZZ holds this matrix in compressed form, with the first nrf columns holding the part of the matrix corresponding to the first level of the overall random subject variable, the next nrf columns holding the part of the matrix corresponding to the second level of the overall random subject variable etc.
If ${\mathbf{pdczz}}=0$, czz is not referenced and may be NULL.
12: $\mathbf{pdczz}$Integer Input
On entry: the stride separating matrix row elements in the array czz.
Constraints:
• ${\mathbf{pdczz}}=0$ or ${\mathbf{pdczz}}\ge {\mathbf{nrf}}$;
• if ${\mathbf{pdczz}}=0$, ${\mathbf{pdcxx}}={\mathbf{pdcxz}}=0$.
(see g02jfc)
13: $\mathbf{cxx}\left[{\mathbf{pdcxx}}×{\mathbf{fnlsv}}×{\mathbf{nff}}\right]$double Output
Note: where ${\mathbf{CXX}}\left(i,j\right)$ appears in this document, it refers to the array element ${\mathbf{cxx}}\left[\left(j-1\right)×{\mathbf{pdcxx}}+i-1\right]$.
On exit: if ${\mathbf{fnlsv}}=1$, CXX holds the lower triangular portion of the matrix $\left(1/{\sigma }^{2}\right)\left({X}^{\mathrm{T}}{\stackrel{^}{V}}^{-1}X\right)$, where $\stackrel{^}{V}$ is the estimated value of $V$.
If ${\mathbf{fnlsv}}>1$, then CXX holds this matrix in compressed form, with the first nff columns holding the part of the matrix corresponding to the first level of the overall fixed subject variable, the next nff columns holding the part of the matrix corresponding to the second level of the overall fixed subject variable, etc.
If ${\mathbf{pdcxx}}=0$, cxx is not referenced and may be NULL.
14: $\mathbf{pdcxx}$Integer Input
On entry: the stride separating matrix row elements in the array cxx.
Constraint: ${\mathbf{pdcxx}}=0$ or ${\mathbf{pdcxx}}\ge {\mathbf{nff}}$.  (see g02jfc)
15: $\mathbf{cxz}\left[{\mathbf{pdcxz}}×{\mathbf{fnlsv}}×{\mathbf{nrf}}\right]$double Output
Note: where ${\mathbf{CXZ}}\left(i,j\right)$ appears in this document, it refers to the array element ${\mathbf{cxz}}\left[\left(j-1\right)×{\mathbf{pdcxz}}+i-1\right]$.
On exit: CXZ holds the matrix $\left(1/{\sigma }^{2}\right)\left({X}^{\mathrm{T}}{\stackrel{^}{V}}^{-1}Z\right)\stackrel{^}{G}$, where $\stackrel{^}{V}$ and $\stackrel{^}{G}$ are the estimates of $V$ and $G$ respectively.
If ${\mathbf{pdcxz}}=0$, cxz is not referenced and may be NULL.
16: $\mathbf{pdcxz}$Integer Input
On entry: the stride separating matrix row elements in the array cxz.
Constraint: ${\mathbf{pdcxz}}=0$ or ${\mathbf{pdcxz}}\ge {\mathbf{fnlsv}}×{\mathbf{nff}}$.  (see )
17: $\mathbf{rcomm}\left[\mathit{dim}\right]$double Communication Array
On entry: communication array initialized by a call to g02jfc.
18: $\mathbf{icomm}\left[\mathit{dim}\right]$Integer Communication Array
On entry: communication array initialized by a call to g02jfc.
19: $\mathbf{fail}$NagError * Input/Output
The NAG error argument (see Section 7 in the Introduction to the NAG Library CL Interface).

## 6Error Indicators and Warnings

NE_ALLOC_FAIL
Dynamic memory allocation failed.
See Section 3.1.2 in the Introduction to the NAG Library CL Interface for further information.
NE_ARRAY_SIZE
On entry, ${\mathbf{lb}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{lb}}\ge ⟨\mathit{\text{value}}⟩$.
On entry, ${\mathbf{pdcxx}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{pdcxx}}\ge ⟨\mathit{\text{value}}⟩$.
On entry, ${\mathbf{pdcxz}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{pdcxz}}\ge ⟨\mathit{\text{value}}⟩$.
On entry, ${\mathbf{pdczz}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{pdczz}}\ge ⟨\mathit{\text{value}}⟩$.
On entry, argument $⟨\mathit{\text{value}}⟩$ had an illegal value.
NE_HANDLE
hlmm has not been initialized or is corrupt.
hlmm is not a G22 handle as generated by g02jfc.
NE_ILLEGAL_COMM
On entry, the communication arrays, icomm and rcomm, have not been initialized correctly.
NE_INT
On entry, ${\mathbf{nvpr}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{nvpr}}\ge ⟨\mathit{\text{value}}⟩$.
NE_INTERNAL_ERROR
An internal error has occurred in this function. Check the function call and any array sizes. If the call is correct then please contact NAG for assistance.
See Section 7.5 in the Introduction to the NAG Library CL Interface for further information.
NE_NEG_ELEMENT
At least one negative estimate for gamma was obtained. All negative estimates have been set to zero.
NE_NO_LICENCE
Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library CL Interface for further information.
NE_REAL_ARRAY
On entry, ${\mathbf{gamma}}\left[⟨\mathit{\text{value}}⟩\right]=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{gamma}}\left[0\right]=-1.0$ or ${\mathbf{gamma}}\left[i-1\right]\ge 0.0$.
NW_KT_CONDITIONS
Current point cannot be improved upon.
NW_NOT_CONVERGED
Optimal solution found, but requested accuracy not achieved.
NW_TOO_MANY_ITER
Too many major iterations.

Not applicable.

## 8Parallelism and Performance

g02jhc is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
g02jhc makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.
Please consult the X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this function. Please also consult the Users' Note for your implementation for any additional implementation-specific information.

If g02jgc has been called, then rather than use the values of fnlsv, nff, rnlsv, nrf and nvpr as returned by g02jfc they should be obtained by calling g02zlc. See Section 9 in g02jgc for more details.

## 10Example

This example fits a random effects model with three random submodels and two fixed effects to a simulated dataset with $90$ observations and $12$ variables. The model is fit using restricted maximum likelihood (REML). Custom labels for the parameter estimates and variance components are then constructed from information returned by g22ydc. See g02jfc for an example that uses the standard parameter labels directly.

### 10.1Program Text

Program Text (g02jhce.c)

### 10.2Program Data

Program Data (g02jhce.d)

### 10.3Program Results

Program Results (g02jhce.r)