sreg(x, y, lambda = NA, df = NA, offset = 0,
weights = rep(1, length(x)), cost = 1,
nstep.cv = 80, tol=1e-5,find.diagA = TRUE, trmin = 2.01,
trmax = NA, lammin = NA,
lammax = NA, verbose = FALSE,
do.cv = TRUE, method = "GCV", rmse = NA,
na.rm = TRUE)ESTIMATE: A smoothing spline is a locally weighted average of the y's based on the relative locations of the x values. Formally the estimate is the curve that minimizes the criterion:
(1/n) sum(k=1,n) w.k( Y.k - f( X.k))**2 + lambda R(f)
where R(f) is the integral of the squared second derivative of f over the range of the X values. Because of the inclusion of the (1/n) in the sum of squares the lambda parameter in sreg corresponds to the a value of lambda*n in the Tps function and in the Krig function. The solution to this minimization is a piecewise cubic polynomial with the join points at the unique set of X values. The polynomial segments are constructed so that the entire curve has continuous first and second derivatives and the second and third derivatives are zero at the boundaries. The smoothing has the range [0,infinity]. Lambda equal to zero gives a cubic spline interpolation of the data. As lambda diverges to infinity ( e.g lambda =1e20) the estimate will converge to the straight line estimated by least squares.
The values of the estimated function at the data points can be expressed in the matrix form: predicted values= A(lambda)Y
where A is an nXn symmetric matrix that does NOT depend on Y. The diagonal elements are the leverage values for the estimate and the sum of these (trace(A(lambda)) can be interpreted as the effective number of parameters that are used to define the spline function. IF there are replicate points the A matrix is the result of finding group averages and applying a weighted spline to the means. The A matrix is also used to find "Bayesian" confidence intervals for the estimate, see the example below.
CROSS-VALIDATION:The GCV criterion with no replicate points for a fixed value of lambda is
(1/n)(Residual sum of squares)/((1-(tr(A)-offset)*cost + offset)/n)**2,
Usually offset =0 and cost =1. Variations on GCV with replicate points are described in the documentation help file for Krig. With an appropriate choice for the smoothing parameter, the estimate of sigma**2 is found by (Residual sum of squares)/tr(A).
COMPUTATIONS: The computations for 1-d splines exploit the banded structure of the matrices needed to solve for the spline coefficients. Banded structure also makes it possible to get the diagonal elements of A quickly. This approach is different from the algorithms in Tps and tremendously more efficient for larger numbers of unique x values ( say > 200). The advantage of Tps is getting "Bayesian" standard errors at predictions different from the observed x values. This function is similar to the S-Plus smooth.spline. The main advantages are more information and control over the choice of lambda and also the FORTRAN source code is available (css.f).
See also the function splint which is designed to be a bare bones
but fast smoothing spline.
# fit a GCV spline to
# control group of rats.
fit<- sreg(rat.diet$t,rat.diet$con)
summary( fit)
set.panel(2,2)
plot(fit) # four diagnostic plots of fit
set.panel()
predict( fit) # predicted values at data points
xg<- seq(0,110,,50)
sm<-predict( fit, xg) # spline fit at 50 equally spaced points
der.sm<- predict( fit, xg, deriv=1) # derivative of spline fit
set.panel( 2,1)
plot( fit$x, fit$y) # the data
lines( xg, sm) # the spline
plot( xg,der.sm, type="l") # plot of estimated derivative
set.panel() # reset panel to 1 plot
# the same fit using the thin plate spline numerical algorithms
# sreg does not scale the obs so instruct Tps not to sacel either
# this will make lambda comparable within factor of n.
fit.tps<-Tps( rat.diet$t,rat.diet$con, scale="unscaled")
summary( fit.tps)
# compare sreg and Tps results to show the adjustment to lambda.
predict( fit)-> look
predict( fit.tps, lambda=fit$lambda*fit$N)-> look2
test.for.zero( look, look2) # silence means it checks to 1e-8
# finding approximate standard errors at observations
SE<- fit$shat.GCV*sqrt(fit$diagA)
# compare to predict.se( fit.tps) differences are due to
# slightly different lambda values and using shat.MLE instad of shat.GCV
#
# 95% pointwise prediction intervals
Zvalue<- qnorm(.0975)
upper<- fit$fitted.values + Zvalue* SE
lower<- fit$fitted.values - Zvalue* SE
#
# conservative, simultaneous Bonferroni bounds
#
ZBvalue<- qnorm(1- .025/fit$N)
upperB<- fit$fitted.values + ZBvalue* SE
lowerB<- fit$fitted.values - ZBvalue* SE
#
# take a look
plot( fit$x, fit$y)
lines( fit$predicted, lwd=2)
matlines( fit$x,
cbind( lower, upper, lowerB, upperB), type="l", col=c( 2,2,4,4), lty=1)
title( "95 pct pointwise and simultaneous intervals")
# or try the more visually honest:
plot( fit$x, fit$y)
lines( fit$predicted, lwd=2)
segments( fit$x, lowerB, fit$x, upperB, col=4)
segments( fit$x, lower, fit$x, upper, col=2, lwd=2)
title( "95 pct pointwise and simultaneous intervals")
set.panel( 1,1)Run the code above in your browser using DataLab