Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox: nag_correg_glm_estfunc (g02gn)

## Purpose

nag_correg_glm_estfunc (g02gn) gives the estimate of an estimable function along with its standard error from the results from fitting a generalized linear model.

## Syntax

[est, stat, sestat, z, ifail] = g02gn(irank, b, covar, v, f, tol, 'ip', ip)
[est, stat, sestat, z, ifail] = nag_correg_glm_estfunc(irank, b, covar, v, f, tol, 'ip', ip)

## Description

nag_correg_glm_estfunc (g02gn) computes the estimates of an estimable function for a generalized linear model which is not of full rank. It is intended for use after a call to nag_correg_glm_normal (g02ga), nag_correg_glm_binomial (g02gb), nag_correg_glm_poisson (g02gc) or nag_correg_glm_gamma (g02gd). An estimable function is a linear combination of the arguments such that it has a unique estimate. For a full rank model all linear combinations of arguments are estimable.
In the case of a model not of full rank the functions use a singular value decomposition (SVD) to find the parameter estimates, $\stackrel{^}{\beta }$, and their variance-covariance matrix. Given the upper triangular matrix $R$ obtained from the $QR$ decomposition of the independent variables the SVD gives
 $R=Q* D 0 0 0 PT,$
where $D$ is a $k$ by $k$ diagonal matrix with nonzero diagonal elements, $k$ being the rank of $R$, and ${Q}_{*}$ and $P$ are $p$ by $p$ orthogonal matrices. This leads to a solution:
 $β^=P1D-1Q*1Tc1,$
${P}_{1}$ being the first $k$ columns of $P$, i.e., $P=\left({P}_{1}{P}_{0}\right)$; ${Q}_{{*}_{1}}$ being the first $k$ columns of ${Q}_{*}$, and ${c}_{1}$ being the first $p$ elements of $c$.
Details of the SVD are made available in the form of the matrix ${P}^{*}$:
 $P*= D-1 P1T P0T$
as described by nag_correg_glm_normal (g02ga), nag_correg_glm_binomial (g02gb), nag_correg_glm_poisson (g02gc) and nag_correg_glm_gamma (g02gd).
A linear function of the arguments, $F={f}^{\mathrm{T}}\beta$, can be tested to see if it is estimable by computing $\zeta ={P}_{0}^{\mathrm{T}}f$. If $\zeta$ is zero, then the function is estimable, if not; the function is not estimable. In practice $\left|\zeta \right|$ is tested against some small quantity $\eta$.
Given that $F$ is estimable it can be estimated by ${f}^{\mathrm{T}}\stackrel{^}{\beta }$ and its standard error calculated from the variance-covariance matrix of $\stackrel{^}{\beta }$, ${C}_{\beta }$, as
 $seF=fTCβf.$
Also a $z$ statistic
 $z=fTβ^ seF ,$
can be computed. The distribution of $z$ will be approximately Normal.

## References

Golub G H and Van Loan C F (1996) Matrix Computations (3rd Edition) Johns Hopkins University Press, Baltimore
McCullagh P and Nelder J A (1983) Generalized Linear Models Chapman and Hall
Searle S R (1971) Linear Models Wiley

## Parameters

### Compulsory Input Parameters

1:     $\mathrm{irank}$int64int32nag_int scalar
$k$, the rank of the dependent variables.
Constraint: $1\le {\mathbf{irank}}\le {\mathbf{ip}}$.
2:     $\mathrm{b}\left({\mathbf{ip}}\right)$ – double array
The ip values of the estimates of the arguments of the model, $\stackrel{^}{\beta }$.
3:     $\mathrm{covar}\left({\mathbf{ip}}×\left({\mathbf{ip}}+1\right)/2\right)$ – double array
The upper triangular part of the variance-covariance matrix of the ip parameter estimates given in b. They are stored packed by column, i.e., the covariance between the parameter estimate given in ${\mathbf{b}}\left(i\right)$ and the parameter estimate given in ${\mathbf{b}}\left(j\right)$, $j\ge i$, is stored in ${\mathbf{covar}}\left(\left(j×\left(j-1\right)/2+i\right)\right)$.
4:     $\mathrm{v}\left(\mathit{ldv},{\mathbf{ip}}+7\right)$ – double array
ldv, the first dimension of the array, must satisfy the constraint $\mathit{ldv}\ge {\mathbf{ip}}$.
5:     $\mathrm{f}\left({\mathbf{ip}}\right)$ – double array
$f$, the linear function to be estimated.
6:     $\mathrm{tol}$ – double scalar
The tolerance value used in the check for estimability, $\eta$.
If ${\mathbf{tol}}\le 0.0$ then $\sqrt{\epsilon }$, where $\epsilon$ is the machine precision, is used instead.

### Optional Input Parameters

1:     $\mathrm{ip}$int64int32nag_int scalar
Default: the dimension of the arrays b, f and the first dimension of the array v. (An error is raised if these dimensions are not equal.)
$p$, the number of terms in the linear model.
Constraint: ${\mathbf{ip}}\ge 1$.

### Output Parameters

1:     $\mathrm{est}$ – logical scalar
Indicates if the function was estimable.
${\mathbf{est}}=\mathit{true}$
The function is estimable.
${\mathbf{est}}=\mathit{false}$
The function is not estimable and stat, sestat and z are not set.
2:     $\mathrm{stat}$ – double scalar
If ${\mathbf{est}}=\mathit{true}$, stat contains the estimate of the function, ${f}^{\mathrm{T}}\stackrel{^}{\beta }$
3:     $\mathrm{sestat}$ – double scalar
If ${\mathbf{est}}=\mathit{true}$, sestat contains the standard error of the estimate of the function, $\mathrm{se}\left(F\right)$.
4:     $\mathrm{z}$ – double scalar
If ${\mathbf{est}}=\mathit{true}$, z contains the $z$ statistic for the test of the function being equal to zero.
5:     $\mathrm{ifail}$int64int32nag_int scalar
${\mathbf{ifail}}={\mathbf{0}}$ unless the function detects an error (see Error Indicators and Warnings).

## Error Indicators and Warnings

Note: nag_correg_glm_estfunc (g02gn) may return useful information for one or more of the following detected errors or warnings.
Errors or warnings detected by the function:

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

${\mathbf{ifail}}=1$
 On entry, ${\mathbf{ip}}<1$, or ${\mathbf{irank}}<1$, or ${\mathbf{irank}}>{\mathbf{ip}}$, or $\mathit{ldv}<{\mathbf{ip}}$.
W  ${\mathbf{ifail}}=2$
 On entry, ${\mathbf{irank}}={\mathbf{ip}}$. In this case est is returned as true and all statistics are calculated.
W  ${\mathbf{ifail}}=3$
Standard error of statistic $\text{}=0.0$; this may be due to rounding errors if the standard error is very small or due to mis-specified inputs covar and f.
${\mathbf{ifail}}=-99$
${\mathbf{ifail}}=-399$
Your licence key may have expired or may not have been installed correctly.
${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.

## Accuracy

The computations are believed to be stable.

The value of estimable functions is independent of the solution chosen from the many possible solutions. While nag_correg_glm_estfunc (g02gn) may be used to estimate functions of the arguments of the model as computed by nag_correg_glm_constrain (g02gk), ${\beta }_{\mathrm{c}}$, these must be expressed in terms of the original arguments, $\beta$. The relation between the two sets of arguments may not be straightforward.

## Example

A loglinear model is fitted to a $3$ by $5$ contingency table by nag_correg_glm_poisson (g02gc). The model consists of terms for rows and columns. The table is:
 $141 67 114 79 39 131 66 143 72 35 36 14 38 28 16$
The number of functions to be tested is read in, then the linear functions themselves are read in and tested with nag_correg_glm_estfunc (g02gn). The results of nag_correg_glm_estfunc (g02gn) are printed.
```function g02gn_example

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

x = [
1  0  0  1  0  0  0  0;
1  0  0  0  1  0  0  0;
1  0  0  0  0  1  0  0;
1  0  0  0  0  0  1  0;
1  0  0  0  0  0  0  1;
0  1  0  1  0  0  0  0;
0  1  0  0  1  0  0  0;
0  1  0  0  0  1  0  0;
0  1  0  0  0  0  1  0;
0  1  0  0  0  0  0  1;
0  0  1  1  0  0  0  0;
0  0  1  0  1  0  0  0;
0  0  1  0  0  1  0  0;
0  0  1  0  0  0  1  0;
0  0  1  0  0  0  0  1];

y = [141  67 114  79  39 131  66 143  72  35  36  14  38  28  16];

[n,m] = size(x);
isx = ones(m,1,'int64');
ip = int64(m+1);

mean_p = 'M';
eps = 1e-6;
tol = 5e-5;

% Fit generalized linear model with Poisson errors
[dev, idf, b, irank, se, covar, v, ifail] = ...
g02gc( ...
link, mean_p, x, isx, ip, y, 'eps', eps, 'tol', tol);

% Display initial results
fprintf('Deviance           = %12.4e\n', dev);
fprintf('Degrees of freedom = %2d\n', idf);
fprintf('\nVariable   Parameter estimate   Standard error\n\n');
ivar = double([1:ip]');
fprintf('%6d%16.4f%20.4f\n',[ivar b se]');

f = [1  0  0;
1  1  1;
0 -1  0;
0  0  0;
1  0  0;
0  0  0;
0  0  0;
0  0  0;
0  0  0];
tol = 5e-05;

% Estimate the estimable functions
for j = 1:size(f,2)
[est, stat, sestat, z, ifail] = ...
g02gn( ...
irank, b, covar, v, f(:,j), tol, 'ip', ip);

% Display results
fprintf('\nFunction %2d\n\n', j);
fprintf('%6.1f', f(:,j)');
if est
fprintf('\n\nstat = %10.4f, se = %10.4f, z = %10.4f\n', stat, sestat, z);
else
fprintf('\n\nFunction not estimable\n');
end
end

```
```g02gn example results

Deviance           =   9.0379e+00
Degrees of freedom =  8

Variable   Parameter estimate   Standard error

1          2.5977              0.0258
2          1.2619              0.0438
3          1.2777              0.0436
4          0.0580              0.0668
5          1.0307              0.0551
6          0.2910              0.0732
7          0.9876              0.0559
8          0.4880              0.0675
9         -0.1996              0.0904

Function  1

1.0   1.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0

stat =     4.8903, se =     0.0674, z =    72.5934

Function  2

0.0   1.0  -1.0   0.0   0.0   0.0   0.0   0.0   0.0

stat =    -0.0158, se =     0.0672, z =    -0.2350

Function  3

0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0

Function not estimable
```