penReg: Function to fit penalised regression

Description

The function penReg() can be used to fit a P-spline. It can be used as demonstration of how the penalised B-splines can be fitted to one explanatory variable. For more that explanatory variable use the function pb() in gamlss.

Usage

penReg(y, x, w = rep(1, length(y)), df = NULL, lambda = NULL, start = 10, 
      inter = 20, order = 2, degree = 3,  plot = FALSE,
      method = c("ML", "ML-1", "GAIC", "GCV", "EM"), k = 2, ...)

Arguments

Value

Returns a fitted object of class penReg. The object contains 1) the fitted coefficients 2) the fitted.values 3) the response variable y, 4) the label of the response variable ylabel 5) the explanatory variable x, 6) the lebel of the explanatory variable 7) the smoothing parameter lambda, 8) the effective degrees of freedom df, 9) the estimete for sigma sigma, 10) the residual sum of squares rss, 11) the Akaike information criterion aic, 12) the Bayesian information criterion sbc and 13) the deviance

References

Eilers, P. H. C. and Marx, B. D. (1996). Flexible smoothing with B-splines and penalties (with comments and rejoinder). Statist. Sci, 11, 89-121. Rigby, R. A. and Stasinopoulos D. M. (2005). Generalized additive models for location, scale and shape,(with discussion), Appl. Statist., 54, part 3, pp 507-554. Stasinopoulos D. M. Rigby R.A. (2007) Generalized additive models for location scale and shape (GAMLSS) in R. Journal of Statistical Software, Vol. 23, Issue 7, Dec 2007, http://www.jstatsoft.org/v23/i07.

Examples

Run this code

set.seed(1234)
 x <- seq(0,10,length=200); y<-(yt<-1+2*x+.6*x^2-.1*x^3)+rnorm(200, 4)
 library(gamlss)
#------------------ 
# df fixed
 g1<-gamlss(y~pb(x, df=4))
 m1<-penReg(y,x, df=4) 
 cbind(g1$mu.coefSmo[[1]]$lambda, m1$lambda)
 cbind(g1$mu.df, m1$df)
 cbind(g1$aic, m1$aic)
 cbind(fitted(g1), fitted(m1))[1:10,]
# identical
#------------------
# estimate lambda using ML
g2<-gamlss(y~pb(x))
m2<-penReg(y,x) 
cbind(g2$mu.df, m2$df)
cbind(g2$mu.lambda, m2$lambda) 
cbind(g2$aic, m2$aic) # different lambda
cbind(fitted(g2), fitted(m2))[1:10,]
# identical
#------------------
#  estimate lambda using GCV
g3 <- gamlss(y~pb(x, method="GCV"))
m3 <- penReg(y,x, method="GCV") 
cbind(g3$mu.df, m3$df)
cbind(g3$mu.lambda, m3$lambda)
cbind(g3$aic, m3$aic)
cbind(fitted(g3), fitted(m3))[1:10,]
# almost identical
#------------------
# estimate lambda using EM
g4<-gamlss(y~pb(x, method="EM"))
m4<-penReg(y,x, method="EM") 
cbind(g4$mu.df, m4$df )
cbind(g4$mu.lambda, m4$lambda)
cbind(g4$aic, m4$aic)
cbind(fitted(g4), fitted(m4))[1:10,]
# almost identical
#------------------
#  estimate lambda using  GAIC(#=3)
g5<-gamlss(y~pb(x, method="GAIC", k=3))
m5<-penReg(y,x, method="GAIC", k=3) 
cbind(g5$mu.df, m5$df )
cbind(g5$mu.lambda, m5$lambda)
cbind(g5$aic, m5$aic)
cbind(g5$mu.df, m5$df)
cbind(g5$mu.lambda, m5$lambda)
cbind(fitted(g5), fitted(m5))[1:10,]
#-------------------
plot(y~x)
lines(fitted(m1)~x, col="green")
lines(fitted(m2)~x, col="red")
lines(fitted(m3)~x, col="blue")
lines(fitted(m4)~x, col="yellow")
lines(fitted(m4)~x, col="grey")