g02 Chapter Contents
g02 Chapter Introduction
NAG Library Manual

NAG Library Function Documentnag_ml_hier_mixed_regsn (g02jec)

1  Purpose

nag_ml_hier_mixed_regsn (g02jec) fits a multi-level linear mixed effects regression model using maximum likelihood (ML). Prior to calling nag_ml_hier_mixed_regsn (g02jec) the initialization function nag_hier_mixed_init (g02jcc) must be called.

2  Specification

 #include #include
 void nag_ml_hier_mixed_regsn (Integer lvpr, const Integer vpr[], Integer nvpr, double gamma[], Integer *effn, Integer *rnkx, Integer *ncov, double *lnlike, Integer lb, Integer id[], Integer pdid, double b[], double se[], double czz[], Integer pdczz, double cxx[], Integer pdcxx, double cxz[], Integer pdcxz, const double rcomm[], const Integer icomm[], const Integer iopt[], Integer liopt, const double ropt[], Integer lropt, NagError *fail)

3  Description

nag_ml_hier_mixed_regsn (g02jec) 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$ by $p$ design matrix for the fixed independent variables, $\beta$ is a vector of length $p$ of unknown fixed effects, $Z$ is a known $n$ by $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$, nag_ml_hier_mixed_regsn (g02jec) 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. nag_ml_hier_mixed_regsn (g02jec) 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).
nag_ml_hier_mixed_regsn (g02jec) fits the model by maximizing the log-likelihood function:
 $-2 l R = log V + n log rT V-1 r + log 2 π / n$
where
 $V = ZG ZT + R, r=y-Xb and b = XT V-1 X -1 XT V-1 y .$
Once the final estimates for ${\gamma }^{*}$ have been obtained, the value of ${\sigma }_{R}^{2}$ is given by
 $σR2 = rT V-1 r / 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).

4  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

5  Arguments

Note: prior to calling nag_ml_hier_mixed_regsn (g02jec) the initialization function nag_hier_mixed_init (g02jcc) must be called, therefore this documention should be read in conjunction with the document for nag_hier_mixed_init (g02jcc).
In particular some argument names and conventions described in that document are also relevant here, but their definition has not been repeated. Specifically, RNDM, wtnnffnrfnlsvlevelsfixed, DAT, licomm and lrcomm should be interpreted identically in both functions.
1:     lvprIntegerInput
On entry: the sum of the number of random parameters and the random intercept flags specified in the call to nag_hier_mixed_init (g02jcc).
Constraint: ${\mathbf{lvpr}}={\sum }_{i}{\mathbf{RNDM}}\left(1,i\right)+{\mathbf{RNDM}}\left(2,i\right)$.
2:     vpr[lvpr]const IntegerInput
On entry: a vector of flags indicating the mapping between the random variables specified in rndm and the variance components, ${\sigma }_{i}^{2}$. See Section 9 for more details.
Constraint: $1\le {\mathbf{vpr}}\left[\mathit{i}-1\right]\le {\mathbf{nvpr}}$, for $\mathit{i}=1,2,\dots ,{\mathbf{lvpr}}$.
3:     nvprIntegerInput
On entry: $g$, the number of variance components being estimated (excluding the overall variance, ${\sigma }_{R}^{2}$).
Constraint: $1\le {\mathbf{nvpr}}\le {\mathbf{lvpr}}$.
4:     gamma[${\mathbf{nvpr}}+1$]doubleInput/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}$.
Constraint: ${\mathbf{gamma}}\left[0\right]=-1.0\text{​ or ​}{\mathbf{gamma}}\left[\mathit{i}-1\right]\ge 0.0$, for $\mathit{i}=1,2,\dots ,g$.
5:     effnInteger *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, then ${\mathbf{effn}}={\mathbf{n}}$.
6:     rnkxInteger *Output
On exit: the rank of the design matrix, $X$, for the fixed effects.
7:     ncovInteger *Output
On exit: number of variance components not estimated to be zero. If none of the variance components are estimated to be zero, then ${\mathbf{ncov}}={\mathbf{nvpr}}$.
8:     lnlikedouble *Output
On exit: $-2{l}_{R}\left(\stackrel{^}{\gamma }\right)$ where ${l}_{R}$ is the log of the maximum likelihood calculated at $\stackrel{^}{\gamma }$, the estimated variance components returned in gamma.
9:     lbIntegerInput
On entry: the dimension of the arrays b and se.
Constraint: ${\mathbf{lb}}\ge {\mathbf{nff}}+{\mathbf{nrf}}×{\mathbf{nlsv}}$.
10:   id[${\mathbf{pdid}}×{\mathbf{lb}}$]IntegerOutput
Note: where ${\mathbf{ID}}\left(i,j\right)$ appears in this document, it refers to the array element ${\mathbf{id}}\left[\left(j-1\right)×{\mathbf{pdid}}+i-1\right]$.
On exit: an array describing the parameter estimates returned in b. The first ${\mathbf{nlsv}}×{\mathbf{nrf}}$ columns of ID describe the parameter estimates for the random effects and the last nff columns the parameter estimates for the fixed effects.
A print function for decoding the parameter estimates given in b using information from id is supplied with the example program for this function.
For fixed effects:
• for $l={\mathbf{nrf}}×{\mathbf{nlsv}}+1,\dots ,{\mathbf{nrf}}×{\mathbf{nlsv}}+{\mathbf{nff}}$
• if ${\mathbf{b}}\left[l-1\right]$ contains the parameter estimate for the intercept then
 $ID1,l = ID2,l = ID3,l = 0 ;$
• if ${\mathbf{b}}\left[l-1\right]$ 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 DAT when ${\mathbf{fixed}}\left[j+1\right]=k$ then
 $ID1,l=0, ID2,l=j, ID3,l=i;$
• if the $j$th variable is continuous or binary, that is ${\mathbf{levels}}\left[{\mathbf{fixed}}\left[j+1\right]-1\right]=1$, then ${\mathbf{ID}}\left(3,l\right)=0$;
• any remaining rows of the $l$th column of ID 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}}={\mathbf{RNDM}}\left(1,b\right)$;
• ${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 DAT when ${\mathbf{RNDM}}\left(2+j,b\right)=k$;
• ${N}_{{S}_{b}}$ denote the number of subject variables in the $b$th random statement, that is ${N}_{{S}_{b}}={\mathbf{RNDM}}\left(3+{N}_{{R}_{b}},b\right)$;
• ${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 DAT when ${\mathbf{RNDM}}\left(3+{N}_{{R}_{b}}+j,b\right)=k$;
• $L\left({S}_{jb}\right)$ denote the number of levels for ${S}_{jb}$, that is $L\left({S}_{jb}\right)={\mathbf{levels}}\left[{\mathbf{RNDM}}\left(3+{N}_{{R}_{b}}+j,b\right)-1\right]$;
• then
• for $l=1,2,\dots {\mathbf{nrf}}×{\mathbf{nlsv}}$, if ${\mathbf{b}}\left[l-1\right]$ contains the parameter estimate for the $i$th level of ${R}_{jb}$ when ${S}_{\mathit{k}b}={s}_{\mathit{k}}$, for $\mathit{k}=1,2,\dots ,{N}_{{S}_{b}}$ and $1\le {s}_{\mathit{k}}\le L\left({S}_{jb}\right)$, i.e., ${s}_{k}$ is a valid value for the $k$th subject variable, then
 $ID1,l=b, ID2,l=j, ID3,l=i, ID3+k,l=sk, ​k=1,2,…,NSb;$
• if the parameter being estimated is for the intercept then ${\mathbf{ID}}\left(2,l\right)={\mathbf{ID}}\left(3,l\right)=0$;
• if the $j$th variable is continuous, or binary, that is $L\left({S}_{jb}\right)=1$, then ${\mathbf{ID}}\left(3,l\right)=0$;
• the remaining rows of the $l$th column of ID 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 ID are set to $-999$.
11:   pdidIntegerInput
On entry: the stride separating matrix row elements in the array id.
Constraint: ${\mathbf{pdid}}\ge 3+\underset{j}{\mathrm{max}}\phantom{\rule{0.25em}{0ex}}\left({\mathbf{RNDM}}\left(3+{\mathbf{RNDM}}\left(1,j\right),j\right)\right)$, i.e., $3+\text{}$ maximum number of subject variables (see nag_hier_mixed_init (g02jcc)).
12:   b[lb]doubleOutput
On exit: the parameter estimates, with the first ${\mathbf{nrf}}×{\mathbf{nlsv}}$ elements of ${\mathbf{b}}$ containing the parameter estimates for the random effects, $\nu$, and the remaining nff elements containing the parameter estimates for the fixed effects, $\beta$. The order of these estimates are described by the id argument.
13:   se[lb]doubleOutput
On exit: the standard errors of the parameter estimates given in b.
14:   czz[$\mathit{dim}$]doubleOutput
Note: the dimension, dim, of the array czz must be at least ${\mathbf{pdczz}}×{\mathbf{nrf}}×{\mathbf{nlsv}}$.
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{nlsv}}=1$, then 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{nlsv}}>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 subject variable, the next nrf columns the part corresponding to the second level of the overall subject variable etc.
15:   pdczzIntegerInput
On entry: the stride separating matrix row elements in the array czz.
Constraint: ${\mathbf{pdczz}}\ge {\mathbf{nrf}}$.
16:   cxx[$\mathit{dim}$]doubleOutput
Note: the dimension, dim, of the array cxx must be at least ${\mathbf{pdcxx}}×{\mathbf{nff}}$.
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: CXX holds the lower triangular portion of the matrix $\left(1/{\sigma }^{2}\right){X}^{\mathrm{T}}{\stackrel{^}{V}}^{-1}X$, where $\stackrel{^}{V}$ is the estimated value of $V$.
17:   pdcxxIntegerInput
On entry: the stride separating matrix row elements in the array cxx.
Constraint: ${\mathbf{pdcxx}}\ge {\mathbf{nff}}$.
18:   cxz[$\mathit{dim}$]doubleOutput
Note: the dimension, dim, of the array cxz must be at least ${\mathbf{pdcxz}}×{\mathbf{nlsv}}×{\mathbf{nrf}}$.
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: if ${\mathbf{nlsv}}=1$, then 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{nlsv}}>1$ then CXZ 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 subject variable, the next nrf columns the part corresponding to the second level of the overall subject variable etc.
19:   pdcxzIntegerInput
On entry: the stride separating matrix row elements in the array cxz.
Constraint: ${\mathbf{pdcxz}}\ge {\mathbf{nff}}$.
20:   rcomm[$\mathit{dim}$]const doubleCommunication Array
Note: the dimension, dim, of the array rcomm must be at least ${\mathbf{lrcomm}}$.
On entry: communication array initialized by a call to nag_hier_mixed_init (g02jcc).
21:   icomm[$\mathit{dim}$]const IntegerCommunication Array
Note: the dimension, dim, of the array icomm must be at least ${\mathbf{licomm}}$.
On entry: communication array initialized by a call to nag_hier_mixed_init (g02jcc).
22:   iopt[liopt]const IntegerInput
On entry: optional arguments passed to the optimization function.
By default nag_ml_hier_mixed_regsn (g02jec) fits the specified model using a modified Newton optimization algorithm as implemented in the NAG Fortran Library routine E04LBF. 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 the NAG Fortran Library routine E04UCF can be chosen by setting ${\mathbf{iopt}}\left[4\right]=1$. If ${\mathbf{liopt}}<4$ or ${\mathbf{iopt}}\left[4\right]\ne 1$ then E04LBF will be used.
Different optional arguments 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 liopt values of iopt are used, so for example, if only the first element of iopt needs changing and default values for all other optional arguments are sufficient liopt can be set to $1$.
NAG Fortran Library routine E04LBF is being used
 $i$ Description Equivalent E04LBF argument Default Value $0$ Number of iterations MAXCAL $1000$ $1$ Unit number for monitoring information. See nag_open_file (x04acc) for details on how to assign a file to a unit number. n/a Output sent to stdout $2$ Print optional arguments ($1=$ print) n/a $-1$ (no printing performed) $3$ Frequency that monitoring information is printed IPRINT $-1$ $4$ Optimizer used n/a n/a
If requested, monitoring information is displayed in a similar format to that given by E04LBF.
NAG Fortran Library routine E04UCF is being used
 $i$ Description Equivalent E04UCF argument Default Value $0$ Number of iterations Major Iteration Limit $\mathrm{max}\left(50,3×{\mathbf{nvpr}}\right)$ $1$ Unit number for monitoring information. See nag_open_file (x04acc) for details on how to assign a file to a unit number. n/a Output sent to stdout $2$ Print optional arguments ($1=\text{}$ 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 $\mathrm{max}\left(50,3×{\mathbf{nvpr}}\right)$ $6$ Frequency that additional monitoring information is printed Minor Print Level $0$
If ${\mathbf{liopt}}\le 0$ then default values are used for all optional arguments and iopt may be set to NULL.
23:   lioptIntegerInput
On entry: length of the options array iopt.
24:   ropt[lropt]const doubleInput
On entry: optional arguments passed to the optimization function.
Different optional arguments 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 lropt values of ropt are used, so for example, if only the first element of ropt needs changing and default values for all other optional arguments are sufficient lropt can be set to $1$.
NAG Fortran Library routine E04LBF is being used
 $i$ Description Equivalent E04LBF argument Default Value $0$ Sweep tolerance n/a $\mathrm{max}\left(\sqrt{\mathrm{eps}},\sqrt{\mathrm{eps}}×\underset{i}{\mathrm{max}}\phantom{\rule{0.25em}{0ex}}\left({\mathrm{zz}}_{ii}\right)\right)$ $1$ Accuracy of linear minimizations ETA $0.9$ $2$ Accuracy to which solution is required XTOL $0.0$ $3$ Initial distance from solution STEPMX $100000.0$
NAG Fortran Library routine E04UCF is being used
 $i$ Description Equivalent E04UCF argument Default Value $0$ Sweep tolerance n/a $\mathrm{max}\left(\sqrt{\mathrm{eps}},\sqrt{\mathrm{eps}}×\underset{i}{\mathrm{max}}\phantom{\rule{0.25em}{0ex}}\left({\mathrm{zz}}_{ii}\right)\right)$ $1$ Lower bound for ${\gamma }^{*}$ n/a $\mathrm{eps}/100$ $2$ Upper bound for ${\gamma }^{*}$ n/a ${10}^{20}$ $3$ Line search tolerance Line Search Tolerance $0.9$ $4$ Optimality tolerance Optimality Tolerance ${\mathrm{eps}}^{0.72}$
where $\mathrm{eps}$ is the machine precision returned by nag_machine_precision (X02AJC) and ${\mathrm{zz}}_{ii}$ denotes the $i$ diagonal element of ${Z}^{\mathrm{T}}Z$.
If ${\mathbf{lropt}}\le 0$ then default values are used for all optional arguments and ropt and may be set to NULL.
25:   lroptIntegerInput
On entry: length of the options array ropt.
26:   failNagError *Input/Output
The NAG error argument (see Section 3.6 in the Essential Introduction).

6  Error Indicators and Warnings

NE_ALLOC_FAIL
Dynamic memory allocation failed.
On entry, argument $⟨\mathit{\text{value}}⟩$ had an illegal value.
NE_INT
On entry, ${\mathbf{lb}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{lb}}\ge ⟨\mathit{\text{value}}⟩$.
On entry, ${\mathbf{lvpr}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{lvpr}}\ge ⟨\mathit{\text{value}}⟩$.
On entry, ${\mathbf{nvpr}}=⟨\mathit{\text{value}}⟩$.
Constraint: $1\le {\mathbf{nvpr}}\le ⟨\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, ${\mathbf{pdid}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{pdid}}\ge ⟨\mathit{\text{value}}⟩$.
NE_INT_ARRAY
On entry, at least one value of i, for $\mathit{i}=1,2,\dots ,{\mathbf{nvpr}}$, does not appear in vpr.
On entry, icomm has not been initialized correctly.
On entry, ${\mathbf{vpr}}\left[⟨\mathit{\text{value}}⟩\right]=⟨\mathit{\text{value}}⟩$ and ${\mathbf{nvpr}}=⟨\mathit{\text{value}}⟩$.
Constraint: $1\le {\mathbf{vpr}}\left[i-1\right]\le {\mathbf{nvpr}}$.
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.
NE_NEG_ELEMENT
At least one negative estimate for gamma was obtained. All negative estimates have been set to zero.
NE_REAL_ARRAY
On entry, ${\mathbf{gamma}}\left[⟨\mathit{\text{value}}⟩\right]=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{gamma}}\left[i-1\right]\ge 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.

8  Parallelism and Performance

nag_ml_hier_mixed_regsn (g02jec) is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
nag_ml_hier_mixed_regsn (g02jec) 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.

The argument vpr gives the mapping between the random variables and the variance components. In most cases ${\mathbf{vpr}}\left[\mathit{i}-1\right]=\mathit{i}$, for $\mathit{i}=1,2,\dots ,{\sum }_{\mathit{i}}{\mathbf{RNDM}}\left(1,\mathit{i}\right)+{\mathbf{RNDM}}\left(2,\mathit{i}\right)$. However, in some cases it might be necessary to associate more than one random variable with a single variance component, for example, when the columns of DAT hold dummy variables.
Consider a dataset with three variables:
 $DAT= 113.6 214.5 311.1 128.3 227.2 326.1$
where the first column corresponds to a categorical variable with three levels, the next to a categorical variable with two levels and the last column to a continuous variable. So in a call to nag_hier_mixed_init (g02jcc)
 $levels=321$
also assume a model with no fixed effects, no random intercept, no nesting and all three variables being included as random effects, then
 $fixed=00; RNDM=30123T.$
Each of the three columns in DAT therefore correspond to a single variable and hence there are three variance components, one for each random variable included in the model, so
 $vpr=123.$
This is the recommended way of supplying the data to nag_ml_hier_mixed_regsn (g02jec), however it is possible to reformat the above dataset by replacing each of the categorical variables with a series of dummy variables, one for each level. The dataset then becomes
 $DAT= 100103.6 010104.5 001101.1 100018.3 010017.2 001016.1$
where each column only has one level
 $levels= 111111 .$
Again a model with no fixed effects, no random intercept, no nesting and all variables being included as random effects is required, so
 $fixed=00 ; RNDM= 60123456T .$
With the data entered in this manner, the first three columns of DAT correspond to a single variable (the first column of the original dataset) as do the next two columns (the second column of the original dataset). Therefore vpr must reflect this
 $vpr= 111223 .$
In most situations it is more efficient to supply the data to nag_hier_mixed_init (g02jcc) in terms of categorical variables rather than transform them into dummy variables.

10  Example

This example fits a random effects model with three levels of nesting to a simulated dataset with $90$ observations and $12$ variables.

10.1  Program Text

Program Text (g02jece.c)

10.2  Program Data

Program Data (g02jece.d)

10.3  Program Results

Program Results (g02jece.r)