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_linregm_estfunc (g02dn)

Purpose

nag_correg_linregm_estfunc (g02dn) gives the estimate of an estimable function along with its standard error.

Syntax

[est, stat, sestat, t, ifail] = g02dn(irank, b, covar, p, f, tol, 'ip', ip)
[est, stat, sestat, t, ifail] = nag_correg_linregm_estfunc(irank, b, covar, p, f, tol, 'ip', ip)

Description

nag_correg_linregm_estfunc (g02dn) computes the estimates of an estimable function for a general linear regression model which is not of full rank. It is intended for use after a call to nag_correg_linregm_fit (g02da) or nag_correg_linregm_update (g02dd). 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 gives the 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 given by nag_correg_linregm_fit (g02da) and nag_correg_linregm_update (g02dd).
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 $t$-statistic,
 $t=fTβ^ seF ,$
can be computed. The $t$-statistic will have a Student's $t$-distribution with degrees of freedom as given by the degrees of freedom for the residual sum of squares for the model.

References

Golub G H and Van Loan C F (1996) Matrix Computations (3rd Edition) Johns Hopkins University Press, Baltimore
Hammarling S (1985) The singular value decomposition in multivariate statistics SIGNUM Newsl. 20(3) 2–25
Searle S R (1971) Linear Models Wiley

Parameters

Compulsory Input Parameters

1:     $\mathrm{irank}$int64int32nag_int scalar
$k$, the rank of the independent 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{p}\left({\mathbf{ip}}×{\mathbf{ip}}+2×{\mathbf{ip}}\right)$ – double array
5:     $\mathrm{f}\left({\mathbf{ip}}\right)$ – double array
$f$, the linear function to be estimated.
6:     $\mathrm{tol}$ – double scalar
$\eta$, the tolerance value used in the check for estimability.
${\mathbf{tol}}\le 0.0$
$\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. (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 t 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{t}$ – double scalar
If ${\mathbf{est}}=\mathit{true}$, t contains the $t$-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_linregm_estfunc (g02dn) 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}}$.
W  ${\mathbf{ifail}}=2$
 On entry, ${\mathbf{irank}}={\mathbf{ip}}$. In this case est is returned as true and all statistics are calculated.
${\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_linregm_estfunc (g02dn) may be used to estimate functions of the arguments of the model as computed by nag_correg_linregm_constrain (g02dk), ${\beta }_{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

Data from an experiment with four treatments and three observations per treatment are read in. A model, with a mean term, is fitted by nag_correg_linregm_fit (g02da). The number of functions to be tested is read in, then the linear functions themselves are read in and tested with nag_correg_linregm_estfunc (g02dn). The results of nag_correg_linregm_estfunc (g02dn) are printed.
```function g02dn_example

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

x = [1, 0, 0, 0;
0, 0, 0, 1;
0, 1, 0, 0;
0, 0, 1, 0;
0, 0, 0, 1;
0, 1, 0, 0;
0, 0, 0, 1;
1, 0, 0, 0;
0, 0, 1, 0;
1, 0, 0, 0;
0, 0, 1, 0;
0, 1, 0, 0];
y = [33.63;     39.62;     38.18;     41.46;     38.02;     35.83;
35.99;     36.58;     42.92;     37.80;     40.43;     37.89];

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

% Fit general linear regression model
[rss, idf, b, se, covar, res, h, q, svd, irank, p, wk, ifail] = ...
g02da(mean_p, x, isx, ip, y);

% Display initial results
fprintf('Estimates from g02da\n\n');
fprintf('Residual sum of squares = %12.4e\n', rss);
fprintf('Degrees of freedom      = %4d\n', idf);
fprintf('\nVariable   Parameter estimate   Standard error\n\n');
ivar = double([1:ip]');
fprintf('%6d%20.4e%20.4e\n',[ivar b se]');

% Estimable functions
f = [1 1  0 0 0;
0 1 -1 0 0;
0 1  0 0 0];
nf = size(f,1);
tol = 1e-05;

% Loop over estimable functions
for j = 1:nf
[est, stat, sestat, t, ifail] = ...
g02dn(irank, b, covar, p, f(j,:), tol);

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

```
```g02dn example results

Estimates from g02da

Residual sum of squares =   2.2227e+01
Degrees of freedom      =    8

Variable   Parameter estimate   Standard error

1          3.0557e+01          3.8494e-01
2          5.4467e+00          8.3896e-01
3          6.7433e+00          8.3896e-01
4          1.1047e+01          8.3896e-01
5          7.3200e+00          8.3896e-01

Function  1

1.00    1.00    0.00    0.00    0.00

stat =    36.0033, se =     0.9623, t =    37.4119

Function  2

0.00    1.00   -1.00    0.00    0.00

stat =    -1.2967, se =     1.3610, t =    -0.9528

Function  3

0.00    1.00    0.00    0.00    0.00

Function not estimable
```