fk: A function to fit break points within GAMLSS

Description Usage Arguments Details Value Author(s) References See Also Examples

View source: R/fk-gamlss.R

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

1
2
fk(x, start=NULL, control=fk.control(...), ...) 
fk.control(degree = 1, all.fixed = FALSE, fixed = NULL, base = c("trun", "Bbase"))

Arguments

x

the x-variable

start

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

control

the degree of the spline function fitted

...

for extra arguments

degree

the degree of the based function

all.fixed

whether to fix all parameter

fixed

the fixed break points

base

Which base should be used

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 therefore multiple maximum points can occur. More details about the free knot splines can be found in package 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.

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.

Author(s)

Mikis Stasinopoulos mikis.stasinopoulos@gamlss.org, Bob Rigby

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.org/).

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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
## 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+1.5*(x-knot))
y <- rNO(201, mu=mu, sigma=.5)
## plot the data
plot(y~x, xlim=c(-1,13), ylim=c(3,23))
## fit model using curfit
m1 <- fitFreeKnots(y, x, knots=3, degree=1)
knots(m1)
## 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
#------------------------------------------------
## now negative binomial data 
knot=4
eta1 <- ifelse(x<=knot,0.8+0.08*x,.8+0.08*x+.3*(x-knot))
plot(eta1~x)
set.seed(143)
y <- rNBI(201, mu=exp(eta1), sigma=.1)
da <- data.frame(y=y,x=x)
plot(y~x, data=da)
## getting the break point using profile deviance
n1 <- quote(gamlss(y ~ x+I((x>this)*(x-this)), family=NBI, data=da))
prof.term(n1, min=1, max=9, criterion="GD", start.prev=FALSE)
## now fit the model using fk
g1 <- gamlss(y~fk(x, degree=1, start=c(4)), data=da, family=NBI)
## get the breakpoint
knots(getSmo(g1))
## summary of the gamlss object FreeBreakPointsReg object
getSmo(g1)
## plot fitted model
plot(y~x, data=da)
lines(fitted(g1)~x, data=da, col="red")
#------------------------------------------------
## the aids data as example where things can go wrong
## using fk()
data(aids)
a1<-gamlss(y~x+fk(x, degree=1, start=25)+qrt, data=aids, family=NBI)
knots(getSmo(a1))
# 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,  start.prev=FALSE)
## The Maximum Likelihood estimator is  18.33231 not 17.37064 
## plotting the fit
with(aids, plot(x,y,pch=21,bg=c("red","green3","blue","yellow")[unclass(qrt)]))
lines(fitted(a1)~aids$x)
#-------------------------------------------------

Example output

Loading required package: gamlss.dist
Loading required package: MASS
Loading required package: gamlss
Loading required package: splines
Loading required package: gamlss.data
Loading required package: nlme
Loading required package: parallel
 **********   GAMLSS Version 5.1-2  ********** 
For more on GAMLSS look at http://www.gamlss.org/
Type gamlssNews() to see new features/changes/bug fixes.

Loading required package: mgcv
This is mgcv 1.8-25. For overview type 'help("mgcv-package")'.
Loading required package: nnet

Attaching package: 'nnet'

The following object is masked from 'package:mgcv':

    multinom

Loading required package: rpart
     BP1 
5.011936 
GAMLSS-RS iteration 1: Global Deviance = 1032.275 
GAMLSS-RS iteration 2: Global Deviance = 1030.958 
GAMLSS-RS iteration 3: Global Deviance = 1030.957 
GAMLSS-RS iteration 4: Global Deviance = 1030.957 
GAMLSS-RS iteration 1: Global Deviance = 1021.241 
GAMLSS-RS iteration 2: Global Deviance = 1020.682 
GAMLSS-RS iteration 3: Global Deviance = 1020.682 
GAMLSS-RS iteration 1: Global Deviance = 1012.962 
GAMLSS-RS iteration 2: Global Deviance = 1012.919 
GAMLSS-RS iteration 3: Global Deviance = 1012.919 
GAMLSS-RS iteration 1: Global Deviance = 1013.683 
GAMLSS-RS iteration 2: Global Deviance = 1013.649 
GAMLSS-RS iteration 3: Global Deviance = 1013.649 
GAMLSS-RS iteration 1: Global Deviance = 1020.127 
GAMLSS-RS iteration 2: Global Deviance = 1019.822 
GAMLSS-RS iteration 3: Global Deviance = 1019.822 
GAMLSS-RS iteration 1: Global Deviance = 1027.483 
GAMLSS-RS iteration 2: Global Deviance = 1026.827 
GAMLSS-RS iteration 3: Global Deviance = 1026.827 
GAMLSS-RS iteration 1: Global Deviance = 1034.174 
GAMLSS-RS iteration 2: Global Deviance = 1032.899 
GAMLSS-RS iteration 3: Global Deviance = 1032.898 
GAMLSS-RS iteration 4: Global Deviance = 1032.898 
****************************************************************** 
The Maximum Likelihood estimator is  4.182408 
with a Global Deviance equal to  1012.211 
A  95 % Confidence interval is: ( 2.990368 , 5.595213 ) 
****************************************************************** 
GAMLSS-RS iteration 1: Global Deviance = 1012.078 
GAMLSS-RS iteration 2: Global Deviance = 1012.077 
GAMLSS-RS iteration 3: Global Deviance = 1012.077 
     BP1 
4.150024 

Call:  fitFreeKnots(y = y, x = xvar, weights = w, degree = degree,  
    knots = lambda, fixed = control$fixed, base = control$base) 

Coefficients:
(Intercept)            x       XatBP1  
   -0.90843      0.07212      0.29733  
Estimated Knots: 
 BP1  
4.15  

GAMLSS-RS iteration 1: Global Deviance = 381.5559 
GAMLSS-RS iteration 2: Global Deviance = 380.1881 
GAMLSS-RS iteration 3: Global Deviance = 380.1878 
     BP1 
17.37064 
GAMLSS-RS iteration 1: Global Deviance = 394.1451 
GAMLSS-RS iteration 2: Global Deviance = 392.7581 
GAMLSS-RS iteration 3: Global Deviance = 392.7562 
GAMLSS-RS iteration 4: Global Deviance = 392.7557 
GAMLSS-RS iteration 1: Global Deviance = 392.9952 
GAMLSS-RS iteration 2: Global Deviance = 391.5406 
GAMLSS-RS iteration 3: Global Deviance = 391.539 
GAMLSS-RS iteration 4: Global Deviance = 391.5384 
GAMLSS-RS iteration 1: Global Deviance = 391.8902 
GAMLSS-RS iteration 2: Global Deviance = 390.3677 
GAMLSS-RS iteration 3: Global Deviance = 390.3655 
GAMLSS-RS iteration 4: Global Deviance = 390.3653 
GAMLSS-RS iteration 1: Global Deviance = 390.8327 
GAMLSS-RS iteration 2: Global Deviance = 389.2424 
GAMLSS-RS iteration 3: Global Deviance = 389.24 
GAMLSS-RS iteration 4: Global Deviance = 389.2398 
GAMLSS-RS iteration 1: Global Deviance = 389.8248 
GAMLSS-RS iteration 2: Global Deviance = 388.1668 
GAMLSS-RS iteration 3: Global Deviance = 388.1642 
GAMLSS-RS iteration 4: Global Deviance = 388.1639 
GAMLSS-RS iteration 1: Global Deviance = 388.8692 
GAMLSS-RS iteration 2: Global Deviance = 387.1436 
GAMLSS-RS iteration 3: Global Deviance = 387.1403 
GAMLSS-RS iteration 4: Global Deviance = 387.14 
GAMLSS-RS iteration 1: Global Deviance = 387.9669 
GAMLSS-RS iteration 2: Global Deviance = 386.1731 
GAMLSS-RS iteration 3: Global Deviance = 386.1698 
GAMLSS-RS iteration 4: Global Deviance = 386.1696 
GAMLSS-RS iteration 1: Global Deviance = 387.12 
GAMLSS-RS iteration 2: Global Deviance = 385.2575 
GAMLSS-RS iteration 3: Global Deviance = 385.2544 
GAMLSS-RS iteration 4: Global Deviance = 385.2543 
GAMLSS-RS iteration 1: Global Deviance = 386.3296 
GAMLSS-RS iteration 2: Global Deviance = 384.3982 
GAMLSS-RS iteration 3: Global Deviance = 384.3952 
GAMLSS-RS iteration 4: Global Deviance = 384.3951 
GAMLSS-RS iteration 1: Global Deviance = 385.5971 
GAMLSS-RS iteration 2: Global Deviance = 383.5956 
GAMLSS-RS iteration 3: Global Deviance = 383.5935 
GAMLSS-RS iteration 4: Global Deviance = 383.5932 
GAMLSS-RS iteration 1: Global Deviance = 384.9225 
GAMLSS-RS iteration 2: Global Deviance = 382.8511 
GAMLSS-RS iteration 3: Global Deviance = 382.8483 
GAMLSS-RS iteration 4: Global Deviance = 382.8481 
GAMLSS-RS iteration 1: Global Deviance = 384.194 
GAMLSS-RS iteration 2: Global Deviance = 382.0384 
GAMLSS-RS iteration 3: Global Deviance = 382.0356 
GAMLSS-RS iteration 4: Global Deviance = 382.0355 
GAMLSS-RS iteration 1: Global Deviance = 383.5353 
GAMLSS-RS iteration 2: Global Deviance = 381.2952 
GAMLSS-RS iteration 3: Global Deviance = 381.2926 
GAMLSS-RS iteration 4: Global Deviance = 381.2925 
GAMLSS-RS iteration 1: Global Deviance = 382.9474 
GAMLSS-RS iteration 2: Global Deviance = 380.6226 
GAMLSS-RS iteration 3: Global Deviance = 380.62 
GAMLSS-RS iteration 4: Global Deviance = 380.6199 
GAMLSS-RS iteration 1: Global Deviance = 382.4312 
GAMLSS-RS iteration 2: Global Deviance = 380.0215 
GAMLSS-RS iteration 3: Global Deviance = 380.0186 
GAMLSS-RS iteration 4: Global Deviance = 380.0184 
GAMLSS-RS iteration 1: Global Deviance = 381.9868 
GAMLSS-RS iteration 2: Global Deviance = 379.4907 
GAMLSS-RS iteration 3: Global Deviance = 379.4879 
GAMLSS-RS iteration 4: Global Deviance = 379.4877 
GAMLSS-RS iteration 1: Global Deviance = 381.6142 
GAMLSS-RS iteration 2: Global Deviance = 379.0304 
GAMLSS-RS iteration 3: Global Deviance = 379.0277 
GAMLSS-RS iteration 4: Global Deviance = 379.0275 
GAMLSS-RS iteration 1: Global Deviance = 381.3128 
GAMLSS-RS iteration 2: Global Deviance = 378.6399 
GAMLSS-RS iteration 3: Global Deviance = 378.6372 
GAMLSS-RS iteration 4: Global Deviance = 378.6371 
GAMLSS-RS iteration 1: Global Deviance = 381.0821 
GAMLSS-RS iteration 2: Global Deviance = 378.3177 
GAMLSS-RS iteration 3: Global Deviance = 378.3153 
GAMLSS-RS iteration 4: Global Deviance = 378.3152 
GAMLSS-RS iteration 1: Global Deviance = 380.9202 
GAMLSS-RS iteration 2: Global Deviance = 378.0629 
GAMLSS-RS iteration 3: Global Deviance = 378.0604 
GAMLSS-RS iteration 4: Global Deviance = 378.0603 
GAMLSS-RS iteration 1: Global Deviance = 380.8259 
GAMLSS-RS iteration 2: Global Deviance = 377.8733 
GAMLSS-RS iteration 3: Global Deviance = 377.8707 
GAMLSS-RS iteration 4: Global Deviance = 377.8706 
GAMLSS-RS iteration 1: Global Deviance = 380.7194 
GAMLSS-RS iteration 2: Global Deviance = 377.6608 
GAMLSS-RS iteration 3: Global Deviance = 377.6579 
GAMLSS-RS iteration 4: Global Deviance = 377.6578 
GAMLSS-RS iteration 1: Global Deviance = 380.6919 
GAMLSS-RS iteration 2: Global Deviance = 377.5254 
GAMLSS-RS iteration 3: Global Deviance = 377.5222 
GAMLSS-RS iteration 4: Global Deviance = 377.5222 
GAMLSS-RS iteration 1: Global Deviance = 380.7412 
GAMLSS-RS iteration 2: Global Deviance = 377.4651 
GAMLSS-RS iteration 3: Global Deviance = 377.4618 
GAMLSS-RS iteration 4: Global Deviance = 377.4618 
GAMLSS-RS iteration 1: Global Deviance = 380.8655 
GAMLSS-RS iteration 2: Global Deviance = 377.4779 
GAMLSS-RS iteration 3: Global Deviance = 377.4746 
GAMLSS-RS iteration 4: Global Deviance = 377.4745 
GAMLSS-RS iteration 1: Global Deviance = 381.0621 
GAMLSS-RS iteration 2: Global Deviance = 377.5616 
GAMLSS-RS iteration 3: Global Deviance = 377.558 
GAMLSS-RS iteration 4: Global Deviance = 377.558 
GAMLSS-RS iteration 1: Global Deviance = 381.3283 
GAMLSS-RS iteration 2: Global Deviance = 377.7133 
GAMLSS-RS iteration 3: Global Deviance = 377.7095 
GAMLSS-RS iteration 4: Global Deviance = 377.7095 
GAMLSS-RS iteration 1: Global Deviance = 381.6609 
GAMLSS-RS iteration 2: Global Deviance = 377.9302 
GAMLSS-RS iteration 3: Global Deviance = 377.9262 
GAMLSS-RS iteration 4: Global Deviance = 377.9261 
GAMLSS-RS iteration 1: Global Deviance = 382.0564 
GAMLSS-RS iteration 2: Global Deviance = 378.2087 
GAMLSS-RS iteration 3: Global Deviance = 378.2049 
GAMLSS-RS iteration 4: Global Deviance = 378.2049 
GAMLSS-RS iteration 1: Global Deviance = 382.5117 
GAMLSS-RS iteration 2: Global Deviance = 378.5468 
GAMLSS-RS iteration 3: Global Deviance = 378.5426 
GAMLSS-RS iteration 4: Global Deviance = 378.5426 
GAMLSS-RS iteration 1: Global Deviance = 383.0228 
GAMLSS-RS iteration 2: Global Deviance = 378.9407 
GAMLSS-RS iteration 3: Global Deviance = 378.9361 
GAMLSS-RS iteration 4: Global Deviance = 378.936 
GAMLSS-RS iteration 1: Global Deviance = 383.5725 
GAMLSS-RS iteration 2: Global Deviance = 379.4074 
GAMLSS-RS iteration 3: Global Deviance = 379.4019 
GAMLSS-RS iteration 4: Global Deviance = 379.4018 
GAMLSS-RS iteration 1: Global Deviance = 384.184 
GAMLSS-RS iteration 2: Global Deviance = 379.9371 
GAMLSS-RS iteration 3: Global Deviance = 379.931 
GAMLSS-RS iteration 4: Global Deviance = 379.9308 
GAMLSS-RS iteration 1: Global Deviance = 384.8531 
GAMLSS-RS iteration 2: Global Deviance = 380.5264 
GAMLSS-RS iteration 3: Global Deviance = 380.5194 
GAMLSS-RS iteration 4: Global Deviance = 380.5192 
GAMLSS-RS iteration 1: Global Deviance = 385.5752 
GAMLSS-RS iteration 2: Global Deviance = 381.1712 
GAMLSS-RS iteration 3: Global Deviance = 381.1633 
GAMLSS-RS iteration 4: Global Deviance = 381.163 
GAMLSS-RS iteration 1: Global Deviance = 386.3461 
GAMLSS-RS iteration 2: Global Deviance = 381.8673 
GAMLSS-RS iteration 3: Global Deviance = 381.8582 
GAMLSS-RS iteration 4: Global Deviance = 381.8582 
GAMLSS-RS iteration 1: Global Deviance = 387.1614 
GAMLSS-RS iteration 2: Global Deviance = 382.6114 
GAMLSS-RS iteration 3: Global Deviance = 382.601 
GAMLSS-RS iteration 4: Global Deviance = 382.601 
GAMLSS-RS iteration 1: Global Deviance = 388.0167 
GAMLSS-RS iteration 2: Global Deviance = 383.3994 
GAMLSS-RS iteration 3: Global Deviance = 383.3874 
GAMLSS-RS iteration 4: Global Deviance = 383.3874 
GAMLSS-RS iteration 1: Global Deviance = 388.9078 
GAMLSS-RS iteration 2: Global Deviance = 384.2275 
GAMLSS-RS iteration 3: Global Deviance = 384.2138 
GAMLSS-RS iteration 4: Global Deviance = 384.2138 
GAMLSS-RS iteration 1: Global Deviance = 389.8308 
GAMLSS-RS iteration 2: Global Deviance = 385.092 
GAMLSS-RS iteration 3: Global Deviance = 385.0764 
GAMLSS-RS iteration 4: Global Deviance = 385.0764 
GAMLSS-RS iteration 1: Global Deviance = 390.7817 
GAMLSS-RS iteration 2: Global Deviance = 385.9893 
GAMLSS-RS iteration 3: Global Deviance = 385.9718 
GAMLSS-RS iteration 4: Global Deviance = 385.9718 
GAMLSS-RS iteration 1: Global Deviance = 391.5393 
GAMLSS-RS iteration 2: Global Deviance = 386.6922 
GAMLSS-RS iteration 3: Global Deviance = 386.6738 
GAMLSS-RS iteration 4: Global Deviance = 386.6738 
GAMLSS-RS iteration 1: Global Deviance = 392.3358 
GAMLSS-RS iteration 2: Global Deviance = 387.4411 
GAMLSS-RS iteration 3: Global Deviance = 387.4204 
GAMLSS-RS iteration 4: Global Deviance = 387.4204 
GAMLSS-RS iteration 1: Global Deviance = 393.1679 
GAMLSS-RS iteration 2: Global Deviance = 388.2305 
GAMLSS-RS iteration 3: Global Deviance = 388.2081 
GAMLSS-RS iteration 4: Global Deviance = 388.208 
GAMLSS-RS iteration 1: Global Deviance = 394.0319 
GAMLSS-RS iteration 2: Global Deviance = 389.0582 
GAMLSS-RS iteration 3: Global Deviance = 389.0337 
GAMLSS-RS iteration 4: Global Deviance = 389.0335 
GAMLSS-RS iteration 1: Global Deviance = 394.9245 
GAMLSS-RS iteration 2: Global Deviance = 389.9202 
GAMLSS-RS iteration 3: Global Deviance = 389.8935 
GAMLSS-RS iteration 4: Global Deviance = 389.8933 
GAMLSS-RS iteration 1: Global Deviance = 395.8423 
GAMLSS-RS iteration 2: Global Deviance = 390.8145 
GAMLSS-RS iteration 3: Global Deviance = 390.7846 
GAMLSS-RS iteration 4: Global Deviance = 390.7844 
GAMLSS-RS iteration 1: Global Deviance = 396.7822 
GAMLSS-RS iteration 2: Global Deviance = 391.7366 
GAMLSS-RS iteration 3: Global Deviance = 391.7039 
GAMLSS-RS iteration 4: Global Deviance = 391.7037 
GAMLSS-RS iteration 1: Global Deviance = 397.7413 
GAMLSS-RS iteration 2: Global Deviance = 392.683 
GAMLSS-RS iteration 3: Global Deviance = 392.6485 
GAMLSS-RS iteration 4: Global Deviance = 392.6483 
GAMLSS-RS iteration 1: Global Deviance = 398.7167 
GAMLSS-RS iteration 2: Global Deviance = 393.6527 
GAMLSS-RS iteration 3: Global Deviance = 393.6155 
GAMLSS-RS iteration 4: Global Deviance = 393.6152 
GAMLSS-RS iteration 1: Global Deviance = 399.7058 
GAMLSS-RS iteration 2: Global Deviance = 394.6422 
GAMLSS-RS iteration 3: Global Deviance = 394.6018 
GAMLSS-RS iteration 4: Global Deviance = 394.6017 
****************************************************************** 
The Maximum Likelihood estimator is  18.33231 
with a Global Deviance equal to  377.4579 
A  95 % Confidence interval is: ( 17.19903 , 19.42018 ) 
****************************************************************** 

gamlss.add documentation built on Feb. 4, 2020, 9:08 a.m.