mgcv: Multiple Smoothing Parameter Estimation by GCV or UBRE

Description

Function to efficiently estimate smoothing parameters in Generalized Ridge Regression Problem with multiple (quadratic) penalties, by GCV or UBRE. The function uses Newton's method in multi-dimensions.

Usage

mgcv(M)

Arguments

is the single argument to mgcv and should be a list with the elements listed below:

M$y

The response data vector.

M$w

A vector of weights for the data (often proportional to the reciprocal of the variance).

M$X

The design matrix for the problem, note that ncol(M$X) must give the number of model parameters, while nrow(M$X) should give the number of data.

M$C

Matrix containing any linear equality constraints on the problem (i.e. $\bf C$ in ${\bf Cp}={\bf 0}$).

M$S

A one dimensional array containing the non-zero elements of the penalty matrices. Let start<-sum(M$df[1:(k-1)]^2) or start<-0 if k==1. Then penalty matrix k has M$S[start+i+M$df[i]*(j-1)

M$off

Offset values locating the elements of M$S[] in the correct
               location within each penalty coefficient matrix.

M$fix

M$fix[i] is FALSE if the corresponding smoothing parameter 
is to be estimated, or TRUE if it is to be fixed at zero.

M$df

M$df[i] gives the number of rows and columns of
               M$S[i] that are non-zero.

M$sig2

This serves 2 purposes. If it is strictly positive then UBRE is 
used
                with sig2 taken as the known error variance. If it is 
non-positive 
                then GCV is used (and an estimate of the error variance 
returned).

M$sp

An array of initial smoothing parameter estimates if any of 
these is not strictly positive then automatic generation of initial values
is used.

M$conv.tol

M$max.half

Each Newton step of the GCV minimisation is halved if it fails to
  improve the GCV score. This is the maximum number of halvings to try
  at each step. 15 is usually enough.

`Value`

The function returns a list of the same form as the input list with the 
following
  elements added/modified:
M$pThe best fit parameters given the estimated smoothing parameters.
M$spThe estimated smoothing parameters 
($\lambda_i/\rho$)
M$sig2The estimated (GCV) or supplied (UBRE) error variance.
M$edfThe estimated degrees of freedom associated with each penalized term.
M$VpEstimated covariance matrix of model parameters.

`WARNING`

The method may not behave well with near column rank deficient ${\bf
X}$
especially in contexts where the weights vary wildly.

`Details`

This is documentation for the code implementing the method described in section 
4 of 
Wood (2000) . The method is a computationally efficient means of applying GCV to 
the 
problem of smoothing parameter selection in generalized ridge regression problems 
of 
the form:
$$minimise~ \| { \bf W}^{1/2} ({ \bf Xp - y} ) \|^2 \rho +  \sum_{i=1}^m
\lambda_i {\bf p^\prime S}_i{\bf p}$$
possibly subject to constraints ${\bf Cp}={\bf 0}$. 
${\bf X}$ is a design matrix, $\bf p$ a parameter vector, 
$\bf y$ a data vector, $\bf W$ a diagonal weight matrix,
${\bf S}_i$ a positive semi-definite matrix  of coefficients
defining the ith penalty and $\bf C$ a matrix of coefficients 
defining any linear equality constraints on the problem. The smoothing
parameters are the $\lambda_i$ but there is an overall
smoothing parameter $\rho$ as well. Note that ${\bf X}$
must be of full column rank, at least when projected  into the null space
of any equality constraints.  
The method operates by alternating very efficient direct searches for 
$\rho$
with Newton updates of the logs of the $\lambda_i$
The method is coded in C and is intended to be portable. It should be 
noted that seriously ill conditioned problems (i.e. with close to column rank 
deficiency in the design matrix) may cause problems, especially if weights vary 
wildly between observations.

`References`

Gu and Wahba (1991) Minimizing GCV/GML scores with multiple smoothing parameters via
the Newton method. SIAM J. Sci. Statist. Comput. 12:383-398

Wood (2000) Modelling and Smoothing Parameter Estimation 
with Multiple 
   Quadratic Penalties. JRSSB 62(2):413-428

http://www.ruwpa.st-and.ac.uk/simon.html

`See Also`

GAMsetup gam

`Examples`

Run this code# This example modified from routine SANtest()
    
    n<-100 # number of observations to simulate
    x <- runif(4 * n, 0, 1) # simulate covariates
    x <- array(x, dim = c(4, n)) # put into array for passing to GAMsetup
    pi <- asin(1) * 2  # begin simulating some data
    y <- 2 * sin(pi * x[1, ])
    y <- y + exp(2 * x[2, ]) - 3.75887
    y <- y + 0.2 * x[3, ]^11 * (10 * (1 - x[3, ]))^6 + 10 * (10 *
        x[3, ])^3 * (1 - x[3, ])^10 - 1.396
    sig2<- -1    # set magnitude of variance
    e <- rnorm(n, 0, sqrt(abs(sig2)))
    y <- y + e          # simulated data
    w <- matrix(1, n, 1) # weight matrix
    par(mfrow = c(2, 2)) # scatter plots of simulated data   
    plot(x[1, ], y)
    plot(x[2, ], y)
    plot(x[3, ], y)
    plot(x[4, ], y)
    G <- list(m = 4, n = n, nsdf = 0, df = c(15, 15, 15, 15),dim=c(1,1,1,1),s.type=c(0,0,0,0), 
       p.order=c(0,0,0,0), x = x) # creat list for passing to GAMsetup
    H <- GAMsetup(G)
    H$y <- y    # add data to H
    H$sig2 <- sig2  # add variance (signalling GCV use in this case) to H
    H$w <- w       # add weights to H
    H$sp<-array(-1,H$m)
    H$fix<-array(FALSE,H$m)
    H$conv.tol<-1e-6;H$max.half<-15
    H <- mgcv(H)  # select smoothing parameters and fit model
Run the code above in your browser using DataLab