NAG Toolbox: nag_correg_mixeff_ml (g02jb)

Purpose

Syntax

Description

References

Parameters

Compulsory Input Parameters

1:     $\mathrm{nvpr}$ – int64int32nag_int scalar

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

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

Constraint: if $\mathbf{nrv}\ne 0$, $1\le \mathbf{nvpr}\le \mathbf{nrv}$.

2:     $\mathrm{levels}\left(\mathbf{ncol}\right)$ – int64int32nag_int array

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

Constraint: $\mathbf{levels}\left(\mathit{i}\right)\ge 1$, for $\mathit{i}=1,2,\dots ,\mathbf{ncol}$.

3:     $\mathrm{yvid}$ – int64int32nag_int scalar

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

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

4:     $\mathrm{fvid}\left(\mathbf{nfv}\right)$ – int64int32nag_int array

The columns of the data matrix dat holding the fixed independent variables with $\mathbf{fvid}\left(i\right)$ holding the column number corresponding to the $i$th fixed variable.

Constraint: $1\le \mathbf{fvid}\left(\mathit{i}\right)\le \mathbf{ncol}$, for $\mathit{i}=1,2,\dots ,\mathbf{nfv}$.

5:     $\mathrm{rvid}\left(\mathbf{nrv}\right)$ – int64int32nag_int array

The columns of the data matrix $\mathbf{dat}$ holding the random independent variables with $\mathbf{rvid}\left(i\right)$ holding the column number corresponding to the $i$th random variable.

Constraint: $1\le \mathbf{rvid}\left(\mathit{i}\right)\le \mathbf{ncol}$, for $\mathit{i}=1,2,\dots ,\mathbf{nrv}$.

6:     $\mathrm{svid}$ – int64int32nag_int scalar

The column of dat holding the subject variable.

If $\mathbf{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 $\mathbf{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 $\mathbf{rint}=1$ the model is equivalent to fitting $Z_S+Z_1\times Z_S+Z_2\times Z_S$ and no subject variable.

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

7:     $\mathrm{cwid}$ – int64int32nag_int scalar

The column of dat holding the case weights.

If $\mathbf{cwid}=0$, no weights are used.

Constraint: $0\le \mathbf{cwid}\le \mathbf{ncol}$.

8:     $\mathrm{vpr}\left(\mathbf{nrv}\right)$ – int64int32nag_int array

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

Constraint: $1\le \mathbf{vpr}\left(\mathit{i}\right)\le \mathbf{nvpr}$, for $\mathit{i}=1,2,\dots ,\mathbf{nrv}$.

9:     $\mathrm{dat}\left(\mathit{lddat},\mathbf{ncol}\right)$ – double array

lddat, the first dimension of the array, must satisfy the constraint $\mathit{lddat}\ge \mathbf{n}$.

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

• $\mathbf{dat}\left(i,\mathbf{yvid}\right)$ holds the dependent variable, $y$;
• if $\mathbf{cwid}\ne 0$, $\mathbf{dat}\left(i,\mathbf{cwid}\right)$ holds the case weights;
• if $\mathbf{svid}\ne 0$, $\mathbf{dat}\left(i,\mathbf{svid}\right)$ holds the subject variable.

The remaining columns hold the values of the independent variables.

Constraints:

• if $\mathbf{cwid}\ne 0$, $\mathbf{dat}\left(i,\mathbf{cwid}\right)\ge 0.0$;
• if $\mathbf{levels}\left(j\right)\ne 1$, $1\le \mathbf{dat}\left(i,j\right)\le \mathbf{levels}\left(j\right)$.

10:   $\mathrm{fint}$ – int64int32nag_int scalar

Flag indicating whether a fixed intercept is included ($\mathbf{fint}=1$).

Constraint: $\mathbf{fint}=0$ or $1$.

11:   $\mathrm{rint}$ – int64int32nag_int scalar

Flag indicating whether a random intercept is included ($\mathbf{rint}=1$).

If $\mathbf{svid}=0$, rint is not referenced.

Constraint: $\mathbf{rint}=0$ or $1$.

12:   $\mathrm{lb}$ – int64int32nag_int scalar

The size of the array b.

Constraint: $\mathbf{lb}\ge \mathbf{fint}+{\displaystyle \sum _{i=1}^{\mathbf{nfv}}}\max \left(\mathbf{levels}\left(\mathbf{fvid}\left(i\right)\right)-1,1\right)+L_S\times \left(\mathbf{rint}+{\displaystyle \sum _{i=1}^{\mathbf{nrv}}}\mathbf{levels}\left(\mathbf{rvid}\left(i\right)\right)\right)$ where $L_S=\mathbf{levels}\left(\mathbf{svid}\right)$ if $\mathbf{svid}\ne 0$ and $1$ otherwise.

13:   $\mathrm{gamma}\left(\mathbf{nvpr}+2\right)$ – double array

Holds the initial values of the variance components, ${\gamma }_0$, with $\mathbf{gamma}\left(\mathit{i}\right)$ the initial value for ${\sigma }_{\mathit{i}}^2/{\sigma }_R^2$, for $\mathit{i}=1,2,\dots ,g$. If $\mathbf{rint}=1$ and $\mathbf{svid}\ne 0$, $g=\mathbf{nvpr}+1$, else $g=\mathbf{nvpr}$.

If $\mathbf{gamma}\left(1\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.

Constraint: $\mathbf{gamma}\left(1\right)=-1.0\text{ or }\mathbf{gamma}\left(\mathit{i}\right)\ge 0.0$, for $\mathit{i}=1,2,\dots ,g$.

Optional Input Parameters

1:     $\mathrm{n}$ – int64int32nag_int scalar

Default: the first dimension of the array dat.

$n$, the number of observations.

Constraint: $\mathbf{n}\ge 1$.

2:     $\mathrm{ncol}$ – int64int32nag_int scalar

Default: the dimension of the array levels and the second dimension of the array dat. (An error is raised if these dimensions are not equal.)

The number of columns in the data matrix, dat.

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

3:     $\mathrm{nfv}$ – int64int32nag_int scalar

Default: the dimension of the array fvid.

The number of independent variables in the model which are to be treated as being fixed.

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

4:     $\mathrm{nrv}$ – int64int32nag_int scalar

Default: the dimension of the arrays rvid, vpr. (An error is raised if these dimensions are not equal.)

The number of independent variables in the model which are to be treated as being random.

Constraints:

• $0\le \mathbf{nrv}<\mathbf{ncol}$;
• $\mathbf{nrv}+\mathbf{rint}>0$.

5:     $\mathrm{maxit}$ – int64int32nag_int scalar

Default: $-1$

The maximum number of iterations.

If $\mathbf{maxit}<0$, the default value of $100$ is used.

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

6:     $\mathrm{tol}$ – double scalar

Default: $0$

The tolerance used to assess convergence.

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

Output Parameters

1:     $\mathrm{nff}$ – int64int32nag_int scalar

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

2:     $\mathrm{nrf}$ – int64int32nag_int scalar

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

3:     $\mathrm{df}$ – int64int32nag_int scalar

The degrees of freedom.

4:     $\mathrm{ml}$ – double scalar

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

5:     $\mathrm{b}\left(\mathbf{lb}\right)$ – double array

The parameter estimates, $\left(\beta ,\nu \right)$, 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 $\mathbf{fint}=1$, $\mathbf{b}\left(1\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=\mathbf{levels}\left(\mathbf{fvid}\left(i\right)\right)$. Define

• if $\mathbf{fint}=1$, $F_1=2$ else if $\mathbf{fint}=0$, $F_1=1$;
• $F_{i+1}=F_i+\max \left(L_i-1,1\right)$, $i\ge 1$.

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

• if $L_i>1$, $\mathbf{b}\left(F_i+\mathit{j}-2\right)$ contains the parameter estimate for the $\mathit{j}$th level of the $i$th fixed variable, for $\mathit{j}=2,3,\dots ,L_i$;
• if $L_i\le 1$, $\mathbf{b}\left(F_i\right)$ 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=\mathbf{levels}\left(\mathbf{rvid}\left(i\right)\right)$. Define

• if $\mathbf{rint}=1$, $R_1=2$ else if $\mathbf{rint}=0$, $R_1=1$;
  $R_{i+1}=R_i+L_i$, $i\ge 1$.

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

• if $\mathbf{svid}=0$,
  • if $L_i>1$, $\mathbf{b}\left(\mathbf{nff}+R_i+\mathit{j}-1\right)$ contains the parameter estimate for the $\mathit{j}$th level of the $i$th random variable, for $\mathit{j}=1,2,\dots ,L_i$;
  • if $L_i\le 1$, $\mathbf{b}\left(\mathbf{nff}+R_i\right)$ contains the parameter estimate for the $i$th random variable;
• if $\mathbf{svid}\ne 0$,
  • let $L_S$ denote the number of levels associated with the subject variable, that is $L_S=\mathbf{levels}\left(\mathbf{svid}\right)$;
  • if $L_i>1$, $\mathbf{b}\left(\mathbf{nff}+\left(\mathit{s}-1\right)L_S+R_i+\mathit{j}-1\right)$ contains the parameter estimate for the interaction between the $\mathit{s}$th level of the subject variable and the $\mathit{j}$th level of the $i$th random variable, for $\mathit{s}=1,2,\dots ,L_S$ and $\mathit{j}=1,2,\dots ,L_i$;
  • if $L_i\le 1$, $\mathbf{b}\left(\mathbf{nff}+\left(\mathit{s}-1\right)L_S+R_i\right)$ contains the parameter estimate for the interaction between the $\mathit{s}$th level of the subject variable and the $i$th random variable, for $\mathit{s}=1,2,\dots ,L_S$;
  • if $\mathbf{rint}=1$, $\mathbf{b}\left(\mathbf{nff}+1\right)$ contains the estimate of the random intercept.

6:     $\mathrm{se}\left(\mathbf{lb}\right)$ – double array

The standard errors of the parameter estimates given in b.

7:     $\mathrm{gamma}\left(\mathbf{nvpr}+2\right)$ – double array

$\mathbf{gamma}\left(\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,g$, holds the final estimate of ${\sigma }_{\mathit{i}}^2$ and $\mathbf{gamma}\left(g+1\right)$ holds the final estimate for ${\sigma }_R^2$.

8:     $\mathrm{warn}$ – int64int32nag_int scalar

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

9:     $\mathrm{ifail}$ – int64int32nag_int scalar

$\mathbf{ifail}=\mathbf{0}$ unless the function detects an error (see Error Indicators and Warnings).

Error Indicators and Warnings

   $\mathbf{ifail}=1$

On entry,	$\mathbf{n}<2$,
or	$\mathbf{ncol}<1$,
or	$\mathit{lddat}<\mathbf{n}$,
or	$\mathbf{yvid}<1$ or $\mathbf{yvid}>\mathbf{ncol}$,
or	$\mathbf{cwid}<0$ or $\mathbf{cwid}>\mathbf{ncol}$,
or	$\mathbf{nfv}<0$ or $\mathbf{nfv}\ge \mathbf{ncol}$,
or	$\mathbf{fint}\ne 0$ and $\mathbf{fint}\ne 1$,
or	$\mathbf{nrv}<0$ or $\mathbf{nrv}>\mathbf{ncol}$ or $\mathbf{nrv}+\mathbf{rint}<1$,
or	$\mathbf{nvpr}<0$ or $\mathbf{nvpr}>\mathbf{nrv}$,
or	$\mathbf{rint}\ne 0$ and $\mathbf{rint}\ne 1$,
or	$\mathbf{svid}<0$ or $\mathbf{svid}>\mathbf{ncol}$,
or	lb is too small.

   $\mathbf{ifail}=2$

On entry,	$\mathbf{levels}\left(i\right)<1$, for at least one $i$,
or	$\mathbf{fvid}\left(i\right)<1$, or $\mathbf{fvid}\left(i\right)>\mathbf{ncol}$, for at least one $i$,
or	$\mathbf{rvid}\left(i\right)<1$, or $\mathbf{rvid}\left(i\right)>\mathbf{ncol}$, for at least one $i$,
or	$\mathbf{vpr}\left(i\right)<1$ or $\mathbf{vpr}\left(i\right)>\mathbf{nvpr}$, for at least one $i$,
or	at least one discrete variable in array dat has a value greater than that specified in levels,
or	$\mathbf{gamma}\left(i\right)<0$, for at least one $i$, and $\mathbf{gamma}\left(1\right)\ne -1$.

   $\mathbf{ifail}=3$

Degrees of freedom $<1$. The number of arguments exceed the effective number of observations.

   $\mathbf{ifail}=4$

The function failed to converge to the specified tolerance in maxit iterations. See Further Comments for advice.

   $\mathbf{ifail}=-99$

An unexpected error has been triggered by this routine. Please contact NAG.

   $\mathbf{ifail}=-399$

Your licence key may have expired or may not have been installed correctly.

   $\mathbf{ifail}=-999$

Dynamic memory allocation failed.

Accuracy

Further Comments

Example

Open in the MATLAB editor: g02jb_example

function g02jb_example


fprintf('g02jb example results\n\n');

dat = [56, 1, 1, 1, 1;
       50, 1, 2, 1, 1;
       39, 1, 3, 1, 1;
       30, 2, 1, 1, 1;
       36, 2, 2, 1, 1;
       33, 2, 3, 1, 1;
       32, 3, 1, 1, 1;
       31, 3, 2, 1, 1;
       15, 3, 3, 1, 1;
       30, 4, 1, 1, 1;
       35, 4, 2, 1, 1;
       17, 4, 3, 1, 1;
       41, 1, 1, 2, 1;
       36, 1, 2, 2, 2;
       35, 1, 3, 2, 3;
       25, 2, 1, 2, 1;
       28, 2, 2, 2, 2;
       30, 2, 3, 2, 3;
       24, 3, 1, 2, 1;
       27, 3, 2, 2, 2;
       19, 3, 3, 2, 3;
       25, 4, 1, 2, 1;
       30, 4, 2, 2, 2;
       18, 4, 3, 2, 3];
[n,ncol] = size(dat);

% Number of levels in each variable
levels = [int64(1);4;3;2;3];

% Model information
yvid =  int64(1);
fvid = [int64(3);  4;  5];
rvid = [int64(3)];
svid =  int64(2);
cwid =  int64(0);
fint =  int64(1);
rint =  int64(1);

% Calculate lb
lb = (rint + sum(levels(rvid)))*prod(levels(svid)) + ...
     fint + sum(levels(fvid)) - numel(fvid);

% Variance component
vpr = [int64(1)];
nvpr = int64(numel(vpr));

% Initial gamma
gamma = [1; 1; 0];
% Fit the linear mixed effects regresion model
[nff, nrf, df, ml, b, se, gamma, warn, ifail] = ...
g02jb( ...
       nvpr, levels, yvid, fvid, rvid, svid, cwid, vpr, ...
       dat, fint, rint, lb, gamma);

% Display results
if warn
   fprintf(['Warning: At least one variance component was ', ...
            'estimated to be negative and then reset to zero\n\n']);
end
fprintf('Fixed effects (Estimate and Standard Deviation)\n\n');
k = 1;
if fint==1
  fprintf('Intercept%15s%10.4f%10.4f\n', ' ',b(k), se(k));
  k = k + 1;
end
for i = 1:numel(fvid)
  for j = 1:levels(fvid(i))
    if levels(fvid(i))==1 || j>1
      fprintf('Variable %3d Level %3d: %10.4f%10.4f\n', i, j, b(k), se(k));
      k = k + 1;
    end
  end
end
fprintf('\nRandom Effects (Estimate and Standard Deviation\n\n');
if svid==0
  for i = 1:numel(rvid)
    for j = 1:levels(rvid(i))
      fprintf('Variable %4d Level %4d: %10.4f %10.4f\n', i, j, b(k), se(k));
      k = k + 1;
    end
  end
else
  for l = 1:levels(svid)
    if (rint==1)
      fprintf('Intercept for Subject Level %4d:%12s%10.4f%10.4f\n', ...
              l, ' ', b(k), se(k));
      k = k + 1;
    end
    for i = 1:numel(rvid)
      for j = 1:levels(rvid(i))
        fprintf('Subject Level %4d Variable %4d Level %4d: %10.4f%10.4f\n', ...
               l, i, j, b(k), se(k));
        k = k + 1;
      end
    end
  end
end
fprintf('\n Variance Components\n');
for i = 1:nvpr+rint
  fprintf('%4d%10.4f\n',i,gamma(i));
end
fprintf('\nsigma^2           = %10.4f\n', gamma(nvpr+rint+1));
fprintf('-2log likelihood  = %10.4f\n', ml);
fprintf('DF                = %16d\n', df);

g02jb example results

Fixed effects (Estimate and Standard Deviation)

Intercept                  37.0000    4.0421
Variable   1 Level   2:     1.0000    3.0461
Variable   1 Level   3:   -11.0000    3.0461
Variable   2 Level   2:    -8.2500    1.8736
Variable   3 Level   2:     0.5000    2.6497
Variable   3 Level   3:     7.7500    2.6497

Random Effects (Estimate and Standard Deviation

Intercept for Subject Level    1:               10.7631    3.8855
Subject Level    1 Variable    1 Level    1:     3.7276    2.6268
Subject Level    1 Variable    1 Level    2:    -1.4476    2.6268
Subject Level    1 Variable    1 Level    3:     0.3733    2.6268
Intercept for Subject Level    2:               -0.5269    3.8855
Subject Level    2 Variable    1 Level    1:    -3.7171    2.6268
Subject Level    2 Variable    1 Level    2:    -1.2253    2.6268
Subject Level    2 Variable    1 Level    3:     4.8125    2.6268
Intercept for Subject Level    3:               -5.6450    3.8855
Subject Level    3 Variable    1 Level    1:     0.5903    2.6268
Subject Level    3 Variable    1 Level    2:     0.3987    2.6268
Subject Level    3 Variable    1 Level    3:    -2.3806    2.6268
Intercept for Subject Level    4:               -4.5912    3.8855
Subject Level    4 Variable    1 Level    1:    -0.6009    2.6268
Subject Level    4 Variable    1 Level    2:     2.2742    2.6268
Subject Level    4 Variable    1 Level    3:    -2.8052    2.6268

 Variance Components
   1   46.7969
   2   11.5365

sigma^2           =     7.0208
-2log likelihood  =   141.6877
DF                =               16