gamlss.add (version 4.0-0)

fk: A function to fit break points within GAMLSS

Description

The fk() function is a additive function to be used for GAMLSS models. It is an interface for the fitFreeKnots() function of package gamlss.util. The functions fitFreeKnots() was first based on the curfit.free.knot() function of package DierckxSpline of Sundar Dorai-Raj and Spencer Graves. The function fk() allows the user to use the free knots function fitFreeKnots() within gamlss. The great advantage of course comes from the fact GAMLSS models provide a variety of distributions and diagnostics.

Usage

fk(x, degree = 3, start = NULL, ...)

Arguments

x

the x-variable

degree

the degree of the spline function fitted

start

starting values for the breakpoints. If are set the number of break points is also determined by the length of start

...

for extra arguments

Value

• The gamlss object saved contains the last fitted object which can be accessed using obj$par.coefSmo where obj is the fitted gamlss object par is the relevant distribution parameter.

Details

Note that fk itself does no smoothing; it simply sets things up for the function gamlss() which in turn uses the function additive.fit() for backfitting which in turn uses gamlss.fk(). Note that, finding the break points is not a trivial problem and therefor multiple maximum points can occur. More details about the free knot splines can be found in Dierckx, (1991). The gamlss algorithm used a modified backfitting in this case, that is, it fits the linear part fist. Note that trying to predict outside the x-range can be dangerous as the example below shows.

References

Dierckx, P. (1991) Curve and Surface Fitting with Splines, Oxford Science Publications Stasinopoulos D. M., Rigby R.A. and Akantziliotou C. (2006) Instructions on how to use the GAMLSS package in R. Accompanying documentation in the current GAMLSS help files, (see also http://www.gamlss.com/). 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.

See Also

gamlss.fk

Examples

library(gamlss.util)
# creating  a linear + linear function
    x <- seq(0,10, length.out=201)
 knot <- 5
 set.seed(12543)
    mu <- ifelse(x<=knot,5+0.5*x,5+0.5*x+(x-knot))
    y <- rNO(201, mu=mu, sigma=.5)
# plot the data
 plot(y~x, xlim=c(-1,13), ylim=c(3,17))
# fit model using curfit
 m1 <- fitFreeKnots(x, y, knots=4, degree=1)
# fitted values
lines(fitted(m1)~x, col="red", lwd="3")
# predict
pm1<-predict(m1, newdata=-1:12)
 points(-1:12,pm1, col="red",pch = 21, bg="blue")
#------------------------------------------------
# now gamlss
#------------------------------------------------
# get the fk function
# fit model
g0 <- gamlss(y~fk(x, degree=1, start=c(4)))
# creating data frame
da <- data.frame(y,x)
g1 <- gamlss(y~fk(x, degree=1, start=c(4)), data=da)
nd<-data.frame(x=-1:12)
pg1<-predict(g1,  newdata=nd)
lines(fitted(g1)~x, col="purple", lwd="3")
points(pg1~I(-1:12), col="darkgreen", pch = 21, bg="green" )
# note the predicting values outside the x-range
# it predicts the global linear fit
#-------------------------------------------------
#-------------------------------------------------
# now negative binomial data 
eta1 <- ifelse(x<=knot,5+0.5*x,5+0.5*x+(x-knot))
 set.seed(143)
 y <- rNBI(201, mu=exp(eta1/7), sigma=.1)
 da <- data.frame(y=y,x=x)
#rm(y,x)
plot(y~x, data=da)
# fitting the model in gamlss using profile deviance
n1 <- quote(gamlss(y ~ x+I((x>this)*(x-this)),family=NBI,data=da))
prof.term(n1, min=1, max=9, step=.1, criterion="GD")
# now fit the model using fk
g1 <- gamlss(y~fk(x, degree=1, start=c(4)), data=da, family=NBI)
# get the breakpoint
knots(g1$mu.coefSmo[[1]])
lines(fitted(g1)~x, data=da)
#------------------------------------------------
# the aids data
# using fk()
data(aids)
a1<-gamlss(y~fk(x, degree=1, start=25)+qrt, data=aids, family=NBI)
knots(a1$mu.coefSmo[[1]])
# using profile deviance
aids.1 <- quote(gamlss(y ~ x+I((x>this)*(x-this))+qrt,family=NBI,data=aids))
prof.term(aids.1, min=16, max=21, step=.1, criterion="GD")
# using non-linear estimation
library(gamlss.nl)
naids <- nl.obj(~p1*(x-p2)*I((x>p2)), start=c(0.2, 20), data=aids)
mod1 <- gamlss(y~nl(naids)+x+qrt, data=aids , family=NBI) 
mod1$mu.coefSmo 
GAIC(mod1,a1) # almost identical
# plotting the fit
with(aids, plot(x,y,pch=21,bg=c("red","green3","blue","yellow")[unclass(qrt)]))
lines(fitted(a1)~aids$x)
#-------------------------------------------------