naginterfaces.library.correg.linregm_​fit

naginterfaces.library.correg.linregm_fit(x, isx, y, mean='M', wt=None, tol=1e-06)

linregm_fit performs a general multiple linear regression when the independent variables may be linearly dependent. Parameter estimates, standard errors, residuals and influence statistics are computed. linregm_fit may be used to perform a weighted regression.

For full information please refer to the NAG Library document for g02da

https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/g02/g02daf.html

Parameters

xfloat, array-like, shape $(n,m)$

$x[\text{i}-1,\text{j}-1]$ must contain the $\text{i}$th observation for the $\text{j}$th independent variable, for $\text{j}=1,2,\dots ,m$, for $\text{i}=1,2,\dots ,n$.

isxint, array-like, shape $\left(m\right)$

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

$isx[j-1]>0$

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

yfloat, array-like, shape $\left(n\right)$

$y$, the observations on the dependent variable.

meanstr, length 1, optional

Indicates if a mean term is to be included.

$mean=\text{'M'}$

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

$mean=\text{'Z'}$

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

wtNone or float, array-like, shape $\left(n\right)$, optional

If provided $wt$ must contain the weights to be used with the model.

If $wt[i-1]=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$.

tolfloat, optional

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 $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.

Returns

rssfloat

The residual sum of squares for the regression.

idfint

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

bfloat, ndarray, shape $\left(\text{ip}\right)$

$b[i-1]$, $i=1,2,\dots ,\text{ip}$ contains the least squares estimates of the parameters of the regression model, $\hat{\beta }$.

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

If $mean=\text{'Z'}$, $b[i-1]$ will contain the coefficient of the variable contained in column $j$ of $x$, where $isx[j-1]$ is the $i$th positive value in the array $isx$.

sefloat, ndarray, shape $\left(\text{ip}\right)$

$se[i-1]$, $i=1,2,\dots ,\text{ip}$ contains the standard errors of the $\text{ip}$ parameter estimates given in $b$.

covfloat, ndarray, shape $(\text{ip}\times (\text{ip}+1)/2)$

The first $\text{ip}\times (\text{ip}+1)/2$ elements of $cov$ contain the upper triangular part of the variance-covariance matrix of the $\text{ip}$ parameter estimates given in $b$. They are stored packed by column, i.e., the covariance between the parameter estimate given in $b[i-1]$ and the parameter estimate given in $b[j-1]$, $j\ge i$, is stored in $cov[j\times (j-1)/2+i-1]$.

resfloat, ndarray, shape $\left(n\right)$

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

hfloat, ndarray, shape $\left(n\right)$

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

qfloat, ndarray, shape $(n,\text{ip}+1)$

The results of the $QR$ decomposition:

the first column of $q$ contains $c$;

the upper triangular part of columns $2$ to $\text{ip}+1$ contain the $R$ matrix;

the strictly lower triangular part of columns $2$ to $\text{ip}+1$ contain details of the $Q$ matrix.

svdbool

If a singular value decomposition has been performed then $svd$ will be $True$, otherwise $svd$ will be $False$.

irankint

The rank of the independent variables.

If $svd=False$, $irank=\text{ip}$.

If $svd=True$, $irank$ is an estimate of the rank of the independent variables.

$irank$ is calculated as the number of singular values greater that $tol\times$ (largest singular value).

It is possible for the SVD to be carried out but $irank$ to be returned as $\text{ip}$.

pfloat, ndarray, shape $(2\times \text{ip}+\text{ip}\times \text{ip})$

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

If $svd=False$, only the first $\text{ip}$ elements of $p$ are used these will contain the zeta values for the $QR$ decomposition (see lapackeig.dgeqrf for details).

If $svd=True$, the first $\text{ip}$ elements of $p$ will contain the zeta values for the $QR$ decomposition (see lapackeig.dgeqrf for details) and the next $\text{ip}$ elements of $p$ contain singular values.

The following $\text{ip}$ by $\text{ip}$ elements contain the matrix $P^{\ast }$ stored by columns.

wkfloat, ndarray, shape $(max(2,5\times (\text{ip}-1)+\text{ip}\times \text{ip}))$

If on exit $svd=True$, $wk$ contains information which is needed by linregm_fit_newvar(); otherwise $wk$ is used as workspace.

Raises

NagValueError

(errno $1$)

On entry, $\text{ip}=⟨value⟩$.

Constraint: $\text{ip}\ge 1$.

(errno $1$)

On entry, $tol=⟨value⟩$.

Constraint: $tol\ge 0.0$.

(errno $1$)

On entry, $\text{ip}=⟨value⟩$ and $n=⟨value⟩$.

Constraint: $\text{ip}\le n$.

(errno $1$)

On entry, $m=⟨value⟩$.

Constraint: $m\ge 1$.

(errno $1$)

On entry, $n=⟨value⟩$.

Constraint: $n\ge 2$.

(errno $2$)

On entry, $\text{weight}=⟨value⟩$.

Constraint: $\text{weight}=\text{'U'}$ or $\text{'W'}$.

(errno $2$)

On entry, $mean=⟨value⟩$.

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

(errno $3$)

On entry, $wt[⟨value⟩]<0.0$.

Constraint: $wt[i-1]\ge 0.0$, for $i=1,2,\dots ,n$.

(errno $4$)

On entry, $isx[⟨value⟩]<0$.

Constraint: $isx[i-1]\ge 0.0$, for $i=1,2,\dots ,m$.

(errno $4$)

On entry, $\text{ip}=⟨value⟩$.

Constraint: $\text{ip}$ must be compatible with the number of nonzero elements in $isx$.

(errno $6$)

SVD solution failed to converge.

Warns

NagAlgorithmicWarning

(errno $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.

Notes

In the NAG Library the traditional C interface for this routine uses a different algorithmic base. Please contact NAG if you have any questions about compatibility.

The general linear regression model is defined by

$$y=X\beta +ϵ\text{,}$$

where

$y$ is a vector of $n$ observations on the dependent variable,

$X$ is an $n\times p$ matrix of the independent variables of column rank $k$,

$\beta$ is a vector of length $p$ of unknown parameters, and

$ϵ$ is a vector of length $n$ of unknown random errors such that $var\left(ϵ\right)=V{\sigma }^2$, where $V$ is a known diagonal matrix.

If $V=I$, the identity matrix, then least squares estimation is used. If $V\ne I$, then for a given weight matrix $W\propto V^{-1}$, weighted least squares estimation is used.

The least squares estimates $\hat{\beta }$ of the parameters $\beta$ minimize ${(y-X\beta )}^T(y-X\beta )$ while the weighted least squares estimates minimize ${(y-X\beta )}^{T}W(y-X\beta )$.

linregm_fit finds a $QR$ decomposition of $X$ (or $W^{1/2}X$ in weighted case), i.e.,

$$X=Q{R}^{\ast }(\text{or}{W}^{1/2}X=Q{R}^{\ast })\text{,}$$

where $R^{\ast }=\left(\begin{array}{c}R\\ 0\end{array}\right)$ and $R$ is a $p\times p$ upper triangular matrix and $Q$ is an $n\times n$ orthogonal matrix. If $R$ is of full rank, then $\hat{\beta }$ is the solution to

$$R\hat{\beta }=c_1\text{,}$$

where $c=Q^{T}y$ (or $Q^T{W}^{1/2}y$) and $c_1$ is the first $p$ elements of $c$. If $R$ is not of full rank a solution is obtained by means of a singular value decomposition (SVD) of $R$,

$$\begin{array}{c}R=Q_{\ast }\left(\begin{array}{cc}D& 0\\ 0& 0\end{array}\right)P^T\text{,}\end{array}$$

where $D$ is a $k\times k$ diagonal matrix with nonzero diagonal elements, $k$ being the rank of $R$, and $Q_{\ast }$ and $P$ are $p\times p$ orthogonal matrices. This gives the solution

$$\hat{\beta }=P_1{D}^{-1}{Q}_{{\ast }_1}^T{c}_1\text{,}$$

$P_1$ being the first $k$ columns of $P$, i.e., $P=\left(\begin{array}{cc}{P}_1& P_0\end{array}\right)$, and $Q_{{\ast }_1}$ being the first $k$ columns of $Q_{\ast }$.

Details of the SVD, are made available, in the form of the matrix $P^{\ast }$:

$$\begin{array}{c}{P}^{\ast }=\left(\begin{array}{c}{D}^{-1}{P}_1^T\\ P_0^T\end{array}\right)\text{.}\end{array}$$

This will be only one of the possible solutions. Other estimates may be obtained by applying constraints to the parameters. These solutions can be obtained by using linregm_constrain() after using linregm_fit. Only certain linear combinations of the parameters will have unique estimates; these are known as estimable functions.

The fit of the model can be examined by considering the residuals, $r_i=y_i-\hat{y}$, where $\hat{y}=X\hat{\beta }$ are the fitted values. The fitted values can be written as $Hy$ for an $n\times n$ matrix $H$. The $i$th diagonal elements of $H$, $h_i$, give a measure of the influence of the $i$th values of the independent variables on the fitted regression model. The values $h_i$ are sometimes known as leverages. Both $r_i$ and $h_i$ are provided by linregm_fit.

The output of linregm_fit also includes $\hat{\beta }$, the residual sum of squares and associated degrees of freedom, $(n-k)$, the standard errors of the parameter estimates and the variance-covariance matrix of the parameter estimates.

In many linear regression models the first term is taken as a mean term or an intercept, i.e., $X_{i,1}=1$, for $i=1,2,\dots ,n$. This is provided as an option. Also only some of the possible independent variables are required to be included in a model, a facility to select variables to be included in the model is provided.

Details of the $QR$ decomposition and, if used, the SVD, are made available. These allow the regression to be updated by adding or deleting an observation using linregm_obs_edit(), adding or deleting a variable using linregm_var_add() and linregm_var_del() or estimating and testing an estimable function using linregm_estfunc(). For the same matrix of independent variables, a new set of parameter estimates can be quickly calculated from a new vector of dependent variables using linregm_fit_newvar(). The details of the factorizations held in $\text{q}$, $\text{p}$ and $\text{wk}$ are only for use by this suite of functions and cannot be used by other functions that use such factorizations, e.g., lapackeig.dormqr since these will expect a different storage scheme for the input factorization.

References

Cook, R D and Weisberg, S, 1982, Residuals and Influence in Regression, Chapman and Hall

Draper, N R and Smith, H, 1985, Applied Regression Analysis, (2nd Edition), Wiley

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

McCullagh, P and Nelder, J A, 1983, Generalized Linear Models, Chapman and Hall

Searle, S R, 1971, Linear Models, Wiley

See also

naginterfaces.library.examples.correg.linregm_fit_ex.main()