NAG Toolbox: nag_correg_linregm_fit (g02da)

Purpose

Syntax

Description

References

Parameters

Compulsory Input Parameters

1:     $\mathrm{mean\_p}$ – string (length ≥ 1)

Indicates if a mean term is to be included.

$\mathbf{mean\_p}=\text{'M'}$

A mean term, intercept, will be included in the model.

$\mathbf{mean\_p}=\text{'Z'}$

The model will pass through the origin, zero-point.

Constraint: $\mathbf{mean\_p}=\text{'M'}$ or $\text{'Z'}$.

2:     $\mathrm{x}\left(\mathit{ldx},\mathbf{m}\right)$ – double array

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

$\mathbf{x}\left(\mathit{i},\mathit{j}\right)$ must contain the $\mathit{i}$th observation for the $\mathit{j}$th independent variable, for $\mathit{i}=1,2,\dots ,\mathbf{n}$ and $\mathit{j}=1,2,\dots ,\mathbf{m}$.

3:     $\mathrm{isx}\left(\mathbf{m}\right)$ – int64int32nag_int array

Indicates which independent variables are to be included in the model.

$\mathbf{isx}\left(j\right)>0$

The variable contained in the $j$th column of x is included in the regression model.

Constraints:

• $\mathbf{isx}\left(\mathit{j}\right)\ge 0$, for $\mathit{j}=1,2,\dots ,m$;
• if $\mathbf{mean\_p}=\text{'M'}$, exactly $\mathbf{ip}-1$ values of isx must be $>0$;
• if $\mathbf{mean\_p}=\text{'Z'}$, exactly ip values of isx must be $>0$.

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

The number of independent variables in the model, including the mean or intercept if present.

Constraints:

• if $\mathbf{mean\_p}=\text{'M'}$, $1\le \mathbf{ip}\le \mathbf{m}+1$;
• if $\mathbf{mean\_p}=\text{'Z'}$, $1\le \mathbf{ip}\le \mathbf{m}$;
• otherwise $1\le \mathbf{ip}\le \mathbf{n}$.

5:     $\mathrm{y}\left(\mathbf{n}\right)$ – double array

$y$, the observations on the dependent variable.

Optional Input Parameters

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

Default: the dimension of the array y and the first dimension of the array x. (An error is raised if these dimensions are not equal.)

$n$, the number of observations.

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

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

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

$m$, the total number of independent variables in the dataset.

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

3:     $\mathrm{wt}\left(:\right)$ – double array

The dimension of the array wt must be at least $\mathbf{n}$ if $\mathit{weight}=\text{'W'}$, and at least $1$ otherwise

If If provided, wt must contain the weights to be used in the weighted regression.

If $\mathbf{wt}\left(i\right)=0.0$, the $i$th observation is not included in the model, in which case the effective number of observations is the number of observations with nonzero weights. The values of res and h will be set to zero for observations with zero weights.

If wt is not provided, the effective number of observations is $n$.

Constraint: if $\mathit{weight}=\text{'W'}$, $\mathbf{wt}\left(\mathit{i}\right)\ge 0.0$, for $\mathit{i}=1,2,\dots ,n$.

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

Default: $0.000001$

The value of tol is used to decide if the independent variables are of full rank and if not what is the rank of the independent variables. The smaller the value of tol the stricter the criterion for selecting the singular value decomposition. If $\mathbf{tol}=0.0$, the singular value decomposition will never be used; this may cause run time errors or inaccurate results if the independent variables are not of full rank.

Constraint: $\mathbf{tol}\ge 0.0$.

Output Parameters

1:     $\mathrm{rss}$ – double scalar

The residual sum of squares for the regression.

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

The degrees of freedom associated with the residual sum of squares.

3:     $\mathrm{b}\left(\mathbf{ip}\right)$ – double array

$\mathbf{b}\left(i\right)$, $i=1,2,\dots ,\mathbf{ip}$ contains the least squares estimates of the parameters of the regression model, $\hat{\beta }$.

If $\mathbf{mean\_p}=\text{'M'}$, $\mathbf{b}\left(1\right)$ will contain the estimate of the mean parameter and $\mathbf{b}\left(i+1\right)$ will contain the coefficient of the variable contained in column $j$ of x, where $\mathbf{isx}\left(j\right)$ is the $i$th positive value in the array isx.

If $\mathbf{mean\_p}=\text{'Z'}$, $\mathbf{b}\left(i\right)$ will contain the coefficient of the variable contained in column $j$ of x, where $\mathbf{isx}\left(j\right)$ is the $i$th positive value in the array isx.

4:     $\mathrm{se}\left(\mathbf{ip}\right)$ – double array

$\mathbf{se}\left(i\right)$, $i=1,2,\dots ,\mathbf{ip}$ contains the standard errors of the ip parameter estimates given in b.

5:     $\mathrm{covar}\left(\mathbf{ip}\times \left(\mathbf{ip}+1\right)/2\right)$ – double array

The first $\mathbf{ip}\times \left(\mathbf{ip}+1\right)/2$ elements of covar contain 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(j\times \left(j-1\right)/2+i\right)$.

6:     $\mathrm{res}\left(\mathbf{n}\right)$ – double array

The (weighted) residuals, $r_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,n$.

7:     $\mathrm{h}\left(\mathbf{n}\right)$ – double array

The diagonal elements of $H$, $h_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,n$.

8:     $\mathrm{q}\left(\mathit{ldq},\mathbf{ip}+1\right)$ – double array

The results of the $QR$ decomposition:

• the first column of q contains $c$;
• the upper triangular part of columns $2$ to $\mathbf{ip}+1$ contain the $R$ matrix;
• the strictly lower triangular part of columns $2$ to $\mathbf{ip}+1$ contain details of the $Q$ matrix.

9:     $\mathrm{svd}$ – logical scalar

If a singular value decomposition has been performed then svd will be true, otherwise svd will be false.

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

The rank of the independent variables.

If $\mathbf{svd}=\mathit{false}$, $\mathbf{irank}=\mathbf{ip}$.

If $\mathbf{svd}=\mathit{true}$, irank is an estimate of the rank of the independent variables.

irank is calculated as the number of singular values greater that $\mathbf{tol}\times$ (largest singular value). It is possible for the SVD to be carried out but irank to be returned as ip.

11:   $\mathrm{p}\left(2\times \mathbf{ip}+\mathbf{ip}\times \mathbf{ip}\right)$ – double array

Details of the $QR$ decomposition and SVD if used.

If $\mathbf{svd}=\mathit{false}$, only the first ip elements of p are used these will contain the zeta values for the $QR$ decomposition (see nag_lapack_dgeqrf (f08ae) for details).

If $\mathbf{svd}=\mathit{true}$, the first ip elements of p will contain the zeta values for the $QR$ decomposition (see nag_lapack_dgeqrf (f08ae) for details) and the next ip elements of p contain singular values. The following ip by ip elements contain the matrix $P^{*}$ stored by columns.

12:   $\mathrm{wk}\left(\max \left(2,5\times \left(\mathbf{ip}-1\right)+\mathbf{ip}\times \mathbf{ip}\right)\right)$ – double array

If on exit $\mathbf{svd}=\mathit{true}$, wk contains information which is needed by nag_correg_linregm_fit_newvar (g02dg); otherwise wk is used as workspace.

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

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

Error Indicators and Warnings

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{n}<2$,
or	$\mathbf{m}<1$,
or	$\mathit{ldx}<\mathbf{n}$,
or	$\mathit{ldq}<\mathbf{n}$,
or	$\mathbf{tol}<0.0$,
or	$\mathbf{ip}\le 0$,
or	$\mathbf{ip}>\mathbf{n}$.

   $\mathbf{ifail}=2$

On entry,	$\mathbf{mean\_p}\ne \text{'M'}$ or $\text{'Z'}$,
or	$\mathit{weight}\ne \text{'W'}$ or $\text{'U'}$.

   $\mathbf{ifail}=3$

On entry,

$\mathit{weight}=\text{'W'}$ and a value of $\mathbf{wt}<0.0$.

   $\mathbf{ifail}=4$

On entry,	a value of $\mathbf{isx}<0$,
or	the value of ip is incompatible with the values of mean_p and isx,
or	ip is greater than the effective number of observations.

W  $\mathbf{ifail}=5$

The degrees of freedom for the residuals are zero, i.e., the designated number of arguments is equal to the effective number of observations. In this case the parameter estimates will be returned along with the diagonal elements of $H$, but neither standard errors nor the variance-covariance matrix will be calculated.

   $\mathbf{ifail}=6$

The singular value decomposition has failed to converge, see nag_eigen_real_triang_svd (f02wu). This is an unlikely error.

   $\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: g02da_example

function g02da_example


fprintf('g02da 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);

% Calculate total sums of squares about mean
[sw, wmean, c, ifail] = g02bu(y);

% Effective number of observations
en = double(idf + irank);
% Calculate R-squared, corrected R-squared and AIC
rsq = 1 - rss/c(1);
mult = (en-1)/double(idf);
arsq = 1 - mult*(1-rsq);
aic = en*log(rss/en) + 2*double(irank);

% Display results
if svd
  fprintf('Model not of full rank, rank = %4d\n\n', irank);
end
fprintf('Residual sum of squares = %12.4e\n', rss);
fprintf('Degrees of freedom      = %4d\n', idf);
fprintf('R-squared               = %12.4e\n', rsq);
fprintf('Adjusted R-squared      = %12.4e\n', arsq);
fprintf('AIC                     = %12.4e\n', aic);
fprintf('\nVariable   Parameter estimate   Standard error\n\n');
ivar = double([1:ip]');
fprintf('%6d%20.4e%20.4e\n',[ivar b se]');
fprintf('\n   Obs          Residuals              H\n\n');
ivar = double([1:n]');
fprintf('%6d%20.4e%20.4e\n',[ivar res h]');

g02da example results

Model not of full rank, rank =    4

Residual sum of squares =   2.2227e+01
Degrees of freedom      =    8
R-squared               =   7.0042e-01
Adjusted R-squared      =   5.8808e-01
AIC                     =   1.5397e+01

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

   Obs          Residuals              H

     1         -2.3733e+00          3.3333e-01
     2          1.7433e+00          3.3333e-01
     3          8.8000e-01          3.3333e-01
     4         -1.4333e-01          3.3333e-01
     5          1.4333e-01          3.3333e-01
     6         -1.4700e+00          3.3333e-01
     7         -1.8867e+00          3.3333e-01
     8          5.7667e-01          3.3333e-01
     9          1.3167e+00          3.3333e-01
    10          1.7967e+00          3.3333e-01
    11         -1.1733e+00          3.3333e-01
    12          5.9000e-01          3.3333e-01