find_best_sp {lazy.stat}

R Documentation

Find the Best Subset of Spline Regression

Description

This function uses "leaps::regsubset" function of "leaps" package.

Usage

find_best_sp(
  object,
  maxnp = 10,
  force.in = NULL,
  eps = 1e-09,
  title = NULL,
  print = 1,
  plot = 1
)

Arguments

object	A fitted model object which contains y The criterion varaible X The regressor matrix yhat The predicted value rmse The rmse = sqrt(RSS/n) x Original univariate x
maxnp	Max size (# of variables) to be searched
force.in	A list of variable numbers to be forced in.
eps	Criterion to judge too small beta coeff.
title	A string to be used as plot title
print	= 0 to surpress the output
plot	= 1 to plot the univariate result

Details

First, using leaps::regsubset, the best subset of regressors of size 1 through maxnp will be found.
Then, using BIC, the best model will be chosen.
Prior to the application of regsubset, those columns of X whose beta values are too small will be removed.

Value

A list of
locsub The location of the colums selected
betasub The regression coefficient vector
Xsub The subset of X matrix: X[,locsub]
ysub The predicted value: X %*% betasub
maxnp Max size specified
np # of original variables
np1 # of selected varaibles: ncol(Xsub) = length(locsub) rmse0 RMSE of the full model rmse1 RMSE of the subset

Examples

# Generate data
set.seed(1701)
n <- 100
errstd <- 0.2
x <- seq( 0, 3, length=n )
yhat <- sin(2*pi*x^(1/2))^3
y <- yhat+errstd*rnorm(n)

# Spline regression of y on x using 11 knots:  b-spline
res <- spreg( y, x, n=11, ylim=c(-1.5,1.5), plot=1, print=1 )
# Find the best subset of size up to 11 variables
resbest <- find_best_sp( res,11 )

# Finding the best subset of knots using p-spline
res <- spreg( y, x, n=11, ylim=c(-1.5,1.5), plot=1, print=1, type="p" )
knots=res$knots

# Find the best subset of size up to 11 variables
resbest <- find_best_sp( res,11, force.in=1:4 )
knots1 <- knots[resbest$locsub[-(1:4)]-4] # assuming force.in=1:4

par(new=TRUE)
text(1.5,-0.9
 , paste("The Original Knots = ", paste(knots, collapse=", ", sep="")
 , sep=""), cex=.8)
par(new=TRUE)
text(1.5,-1.1
   , paste("The Best Knots = ", paste(knots1, collapse=", ", sep="")
   , sep=""), cex=.8)
par(new=FALSE)

---