NAG Library Function Document

nag_regsn_ridge (g02kbc)

nag_regsn_ridge (Nag_OrderType order, Integer n, Integer m, const double x[], Integer pdx, const Integer isx[], Integer ip, const double y[], Integer lh, const double h[], double nep[], Nag_ParaOption wantb, double b[], Integer pdb, Nag_VIFOption wantvf, double vf[], Integer pdvf, Integer lpec, const Nag_PredictError pec[], double pe[], Integer pdpe, NagError *fail)

3 Description

A linear model has the form:

y = c + X β + ε,

where

$y$ is an $n$ by $1$ matrix of values of a dependent variable;
$c$ is a scalar intercept term;
$X$ is an $n$ by $m$ matrix of values of independent variables;
$β$ is a $m$ by $1$ matrix of unknown values of parameters;
$ε$ is an $n$ by $1$ matrix of unknown random errors such that variance of ${ε = σ}^{2} I$ .

Let

\tilde{X}

be the mean-centred

X

and

\tilde{y}

the mean-centred

y

. Furthermore,

\tilde{X}

is scaled such that the diagonal elements of the cross product matrix

{\tilde{X}}^{T} \tilde{X}

are one. The linear model now takes the form:

\tilde{y} = \tilde{X} \tilde{β} + ε .

Ridge regression estimates the parameters

\tilde{β}

in a penalised least squares sense by finding the

\tilde{b}

that minimizes

{‖\tilde{X} \tilde{b} - \tilde{y}‖}^{2} + h {‖\tilde{b}‖}^{2}, h > 0,

where

‖\cdot‖

denotes the

ℓ_{2}

-norm and

h

is a scalar regularization or ridge parameter. For a given value of

h

, the parameters estimates

\tilde{b}

are found by evaluating

\tilde{b} = {({\tilde{X}}^{T} \tilde{X} + h I)}^{- 1} {\tilde{X}}^{T} \tilde{y} .

Note that if

h = 0

the ridge regression solution is equivalent to the ordinary least squares solution.

Rather than calculate the inverse of (

{\tilde{X}}^{T} \tilde{X} + h I

) directly, nag_regsn_ridge (g02kbc) uses the singular value decomposition (SVD) of

\tilde{X}

. After decomposing

\tilde{X}

into

U D V^{T}

where

U

and

V

are orthogonal matrices and

D

is a diagonal matrix, the parameter estimates become

\tilde{b} = V {(D^{T} D + h I)}^{- 1} D U^{T} \tilde{y} .

A consequence of introducing the ridge parameter is that the effective number of parameters,

γ

, in the model is given by the sum of diagonal elements of

D^{T} D {(D^{T} D + h I)}^{- 1},

see Moody (1992) for details.

Any multi-collinearity in the design matrix

X

may be highlighted by calculating the variance inflation factors for the fitted model. The

j

th variance inflation factor,

v_{j}

, is a scaled version of the multiple correlation coefficient between independent variable

j

and the other independent variables,

R_{j}

, and is given by

v_{j} = \frac{1}{1 - R_{j}}, j = 1, 2, \dots, m .

The

m

variance inflation factors are calculated as the diagonal elements of the matrix:

{({\tilde{X}}^{T} \tilde{X} + h I)}^{- 1} {\tilde{X}}^{T} \tilde{X} {({\tilde{X}}^{T} \tilde{X} + h I)}^{- 1},

which, using the SVD of

\tilde{X}

, is equivalent to the diagonal elements of the matrix:

V {(D^{T} D + h I)}^{- 1} D^{T} D {(D^{T} D + h I)}^{- 1} V^{T} .

Given a value of

h

, any or all of the following prediction criteria are available:

(a)

Generalized cross-validation (GCV):

\frac{n s}{{(n - γ)}^{2}};

(b)

Unbiased estimate of variance (UEV):

\frac{s}{n - γ};

(c)

Future prediction error (FPE):

\frac{1}{n} (s + \frac{2 γ s}{n - γ});

(d)

Bayesian information criterion (BIC):

\frac{1}{n} (s + \frac{\log (n) γ s}{n - γ});

(e)

Leave-one-out cross-validation (LOOCV),

where

s

is the sum of squares of residuals.

Although parameter estimates

\tilde{b}

are calculated by using

\tilde{X}

, it is usual to report the parameter estimates

b

associated with

X

. These are calculated from

\tilde{b}

, and the means and scalings of

X

. Optionally, either

\tilde{b}

b

may be calculated.

4 References

Hastie T, Tibshirani R and Friedman J (2003) The Elements of Statistical Learning: Data Mining, Inference and Prediction Springer Series in Statistics

Moody J.E. (1992) The effective number of parameters: An analysis of generalisation and regularisation in nonlinear learning systems In: Neural Information Processing Systems (eds J E Moody, S J Hanson, and R P Lippmann) 4 847–854 Morgan Kaufmann San Mateo CA

5 Arguments

1: order – Nag_OrderTypeInput

On entry: the order argument specifies the two-dimensional storage scheme being used, i.e., row-major ordering or column-major ordering. C language defined storage is specified by

order = Nag_RowMajor

. See Section 3.2.1.3 in the Essential Introduction for a more detailed explanation of the use of this argument.

Constraint:

order = Nag_RowMajor

or Nag_ColMajor.

2: n – IntegerInput

On entry:

n

, the number of observations.

Constraint:

n \geq 1

3: m – IntegerInput

On entry: the number of independent variables available in the data matrix

X

Constraint:

m \leq n

4: x[ $\dim$ ] – const doubleInput

Note: the dimension, dim, of the array x must be at least

$\max (1, pdx \times m)$ when $order = Nag_ColMajor$ ;
$\max (1, n \times pdx)$ when $order = Nag_RowMajor$ .

The

(i, j)

th element of the matrix

X

is stored in

$x [(j - 1) \times pdx + i - 1]$ when $order = Nag_ColMajor$ ;
$x [(i - 1) \times pdx + j - 1]$ when $order = Nag_RowMajor$ .

On entry: the values of independent variables in the data matrix

X

5: pdx – IntegerInput

On entry: the stride separating row or column elements (depending on the value of order) in the array x.

Constraints:

if $order = Nag_ColMajor$ , $pdx \geq n$ ;
if $order = Nag_RowMajor$ , $pdx \geq m$ .

6: isx[m] – const IntegerInput

On entry: indicates which

m

independent variables are included in the model.

$isx [j - 1] = 1$: The $j$ th variable in x will be included in the model.
$isx [j - 1] = 0$: Variable $j$ is excluded.

Constraint:

isx [j - 1] = 0 ​ or ​ 1

, for

j = 1, 2, \dots, m

7: ip – IntegerInput

On entry:

m

, the number of independent variables in the model.

Constraints:

$1 \leq ip \leq m$ ;
exactly ip elements of isx must be equal to $1$ .

8: y[n] – const doubleInput

On entry: the

n

values of the dependent variable

y

9: lh – IntegerInput

On entry: the number of supplied ridge parameters.

Constraint:

lh > 0

10: h[lh] – const doubleInput

On entry:

h [j - 1]

is the value of the

j

th ridge parameter

h

Constraint:

h [j - 1] \geq 0.0

, for

j = 1, 2, \dots, lh

11: nep[lh] – doubleOutput

On exit:

nep [j - 1]

is the number of effective parameters,

γ

, in the

j

th model, for

j = 1, 2, \dots, lh

12: wantb – Nag_ParaOptionInput

On entry: defines the options for parameter estimates.

$wantb = Nag_NoPara$: Parameter estimates are not calculated and b is not referenced.
$wantb = Nag_OrigPara$: Parameter estimates $b$ are calculated for the original data.
$wantb = Nag_StandPara$: Parameter estimates $\tilde{b}$ are calculated for the standardized data.

Constraint:

wantb = Nag_NoPara

Nag_OrigPara

Nag_StandPara

13: b[ $\dim$ ] – doubleOutput

Note: the dimension, dim, of the array b must be at least

$pdb \times lh$ when $wantb \neq Nag_NoPara$ and $order = Nag_ColMajor$ ;
$\max (1, (ip + 1) \times pdb)$ when $wantb \neq Nag_NoPara$ and $order = Nag_RowMajor$ ;
$1$ otherwise.

Where

B (i, j)

appears in this document, it refers to the array element

$b [(j - 1) \times pdb + i - 1]$ when $order = Nag_ColMajor$ ;
$b [(i - 1) \times pdb + j - 1]$ when $order = Nag_RowMajor$ .

On exit: if

wantb \neq Nag_NoPara

, b contains the intercept and parameter estimates for the fitted ridge regression model in the order indicated by isx.

B (1, j)

, for

j = 1, 2, \dots, lh

, contains the estimate for the intercept;

B (i + 1, j)

contains the parameter estimate for the

i

th independent variable in the model fitted with ridge parameter

h [j - 1]

, for

i = 1, 2, \dots, ip

14: pdb – IntegerInput

On entry: the stride separating row or column elements (depending on the value of order) in the array b.

Constraints:

if order=Nag_ColMajor,
- if $wantb \neq Nag_NoPara$ , $pdb \geq ip + 1$ ;
- otherwise $pdb \geq 1$ ;
if order=Nag_RowMajor,
- if $wantb \neq Nag_NoPara$ , $pdb \geq lh$ ;
- otherwise $pdb \geq 1$ .

15: wantvf – Nag_VIFOptionInput

On entry: defines the options for variance inflation factors.

$wantvf = Nag_NoVIF$: Variance inflation factors are not calculated and the array vf is not referenced.
$wantvf = Nag_WantVIF$: Variance inflation factors are calculated.

Constraints:

$wantvf = Nag_NoVIF$ or $Nag_WantVIF$ ;
if $wantb = Nag_NoPara$ , $wantvf = Nag_WantVIF$ .

16: vf[ $\dim$ ] – doubleOutput

Note: the dimension, dim, of the array vf must be at least

$pdvf \times lh$ when $wantvf \neq Nag_NoVIF$ and $order = Nag_ColMajor$ ;
$\max (1, ip \times pdvf)$ when $wantvf \neq Nag_NoVIF$ and $order = Nag_RowMajor$ ;
$1$ otherwise.

Where

VF (i, j)

appears in this document, it refers to the array element

$vf [(j - 1) \times pdvf + i - 1]$ when $order = Nag_ColMajor$ ;
$vf [(i - 1) \times pdvf + j - 1]$ when $order = Nag_RowMajor$ .

On exit: if

wantvf = Nag_WantVIF

, the variance inflation factors. For the

i

th independent variable in a model fitted with ridge parameter

h [j - 1]

VF (i, j)

is the value of

v_{i}

, for

i = 1, 2, \dots, ip

17: pdvf – IntegerInput

On entry: the stride separating row or column elements (depending on the value of order) in the array vf.

Constraints:

if order=Nag_ColMajor,
- if $wantvf \neq Nag_NoVIF$ , $pdvf \geq ip$ ;
- otherwise $pdvf \geq 1$ ;
if order=Nag_RowMajor,
- if $wantvf \neq Nag_NoVIF$ , $pdvf \geq lh$ ;
- otherwise $pdvf \geq 1$ .

18: lpec – IntegerInput

On entry: the number of prediction error statistics to return; set

lpec \leq 0

for no prediction error estimates.

19: pec[ $lpec$ ] – const Nag_PredictErrorInput

On entry: if

lpec > 0

pec [j - 1]

defines the

j

th prediction error, for

j = 1, 2, \dots, lpec

; otherwise pec is not referenced.

$pec [j - 1] = Nag_BIC$: Bayesian information criterion (BIC).
$pec [j - 1] = Nag_FPE$: Future prediction error (FPE).
$pec [j - 1] = Nag_GCV$: Generalized cross-validation (GCV).
$pec [j - 1] = Nag_LOOCV$: Leave-one-out cross-validation (LOOCV).
$pec [j - 1] = Nag_EUV$: Unbiased estimate of variance (UEV).

Constraint: if

lpec > 0

pec [j - 1] = Nag_BIC

Nag_FPE

Nag_GCV

Nag_LOOCV

Nag_EUV

, for

j = 1, 2, \dots, lpec

20: pe[ $\dim$ ] – doubleOutput

Note: the dimension, dim, of the array pe must be at least

$pdpe \times lh$ when $lpec > 0$ and $order = Nag_ColMajor$ ;
$\max (1, lpec \times pdpe)$ when $lpec > 0$ and $order = Nag_RowMajor$ ;
$1$ otherwise.

Where

PE (i, j)

appears in this document, it refers to the array element

$pe [(j - 1) \times pdpe + i - 1]$ when $order = Nag_ColMajor$ ;
$pe [(i - 1) \times pdpe + j - 1]$ when $order = Nag_RowMajor$ .

On exit: if

lpec \leq 0

, pe is not referenced; otherwise

PE (i, j)

contains the prediction error of criterion

pec [i - 1]

for the model fitted with ridge parameter

h [j - 1]

, for

i = 1, 2, \dots, lpec

and

j = 1, 2, \dots, lh

21: pdpe – IntegerInput

On entry: the stride separating row or column elements (depending on the value of order) in the array pe.

Constraints:

if order=Nag_ColMajor,
- if $lpec > 0$ , $pdpe \geq lpec$ ;
- otherwise $pdpe \geq 1$ ;
if order=Nag_RowMajor,
- if $lpec > 0$ , $pdpe \geq lh$ ;
- otherwise $pdpe \geq 1$ .

22: fail – NagError *Input/Output

The NAG error argument (see Section 3.6 in the Essential Introduction).

6 Error Indicators and Warnings

NE_ALLOC_FAIL: Dynamic memory allocation failed.
NE_BAD_PARAM: On entry, argument $〈value〉$ had an illegal value.
NE_CONSTRAINT: On entry, $wantb = Nag_NoPara$ and $wantvf = Nag_NoVIF$ .
NE_ENUM_INT_2: On entry, $pdb = 〈value〉$ and $ip = 〈value〉$ .
Constraint: if $wantb \neq Nag_NoPara$ , $pdb \geq ip + 1$ .; On entry, $pdvf = 〈value〉$ and $ip = 〈value〉$ .
Constraint: if $wantvf \neq Nag_NoVIF$ , $pdvf \geq ip$ .; On entry, $wantb = 〈value〉$ , $pdb = 〈value〉$ , $lh = 〈value〉$ .
Constraint: if $wantb \neq Nag_NoPara$ , $pdb \geq lh$ ;
otherwise $pdb \geq 1$ .; On entry, $wantvf = 〈value〉$ , $pdvf = 〈value〉$ , $lh = 〈value〉$ .
Constraint: if $wantvf \neq Nag_NoVIF$ , $pdvf \geq lh$ ;
otherwise $pdvf \geq 1$ .
NE_INT: On entry, $lh = 〈value〉$ .
Constraint: $lh > 0$ .; On entry, $n = 〈value〉$ .
Constraint: $n \geq 1$ .
NE_INT_2: On entry, $m = 〈value〉$ and $n = 〈value〉$ .
Constraint: $m \leq n$ .; On entry, $pdpe = 〈value〉$ and $lpec = 〈value〉$ .
Constraint: $pdpe \geq lpec$ .; On entry, $pdx = 〈value〉$ and $m = 〈value〉$ .
Constraint: $pdx \geq m$ .; On entry, $pdx = 〈value〉$ and $n = 〈value〉$ .
Constraint: $pdx \geq n$ .
NE_INT_3: On entry, $pdpe = 〈value〉$ , $lpec = 〈value〉$ and $lh = 〈value〉$ .
Constraint: if $lpec > 0$ , $pdpe \geq lh$ ;
otherwise $pdpe \geq 1$ .
NE_INT_ARG_CONS: ip does not equal the sum of elements in isx.
NE_INT_ARRAY_VAL_1_OR_2: On entry, $isx [i - 1] \neq 0$ or $1$ for at least one $i$ .
NE_INTERNAL_ERROR: An internal error has occurred in this function. Check the function call and any array sizes. If the call is correct then please contact NAG for assistance.
NE_REAL_ARRAY_CONS: On entry, $h [i - 1] < 0$ for at least one $i$ .

7 Accuracy

The accuracy of nag_regsn_ridge (g02kbc) is closely related to that of the singular value decomposition.

8 Further Comments

nag_regsn_ridge (g02kbc) allocates internally

\max (5 \times (n - 1), 2 \times ip \times ip) + (n + 3) \times ip + n

elements of double precision storage.

9 Example

This example reads in data from an experiment to model body fat, and a selection of ridge regression models are calculated.

NAG Library Function Documentnag_regsn_ridge (g02kbc)

+− Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Further Comments

9 Example

9.1 Program Text

9.2 Program Data

9.3 Program Results

NAG Library Function Document

nag_regsn_ridge (g02kbc)