A interface functions to use Simon Wood's gam() and bam() functions within GAMLSS

Share:

Description

The ga() and ba() functions are a additive functions to be used within GAMLSS models. They are interfaces for the gam() and the bam() functions of package mgcv of Simon Wood. The functions gam() and the bam() allows the user to use all the available smoothers of the package mcgv() within gamlss. The great advantage of course come from fitting models outside the exponential family.

Usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
ga(formula, control = ga.control(...), ...)

ba(formula, control = ba.control(...), ...)

ga.control(offset = NULL, method = "REML", optimizer = c("outer", 
          "newton"), select = FALSE, knots = NULL,sp = NULL, 
          min.sp = NULL, H = NULL, gamma = 1, paraPen = NULL, 
          in.out = NULL, drop.unused.levels = TRUE,
          drop.intercept = NULL, ...)

ba.control(offset = NULL, method = "fREML", select = FALSE, 
          scale = 0, gamma = 1, knots = NULL, sp = NULL, 
          min.sp = NULL, paraPen = NULL, chunk.size = 10000, 
          rho = 0, AR.start = NULL, discrete = FALSE, 
          cluster = NULL, nthreads = NA, gc.level = 1, 
          use.chol = FALSE, samfrac = 1, drop.unused.levels = TRUE, 
          drop.intercept = NULL, ...)

Arguments

formula

A formula containing s() and te functions i.e. ~s(x1)+ te(x2,x3).

control

this allow to specify argument within the function gam() of mgcv

offset

the offset argument in gam() and bam()

method

the method argument in gam() and bam()

optimizer

the method optimizer in gam()

select

the select argument in gam() and bam()

knots

the knots argument in gam() and bam()

sp

the sp argument in gam() and bam()

min.sp

the min.sp argument in gam() and bam()

H

a user supplied fixed quadratic penalty on the parameters in gam()

gamma

the gamma argument in gam() and bam()

paraPen

the paraPen argument in gam() and bam()

in.out

the in.out argument in gam()

drop.unused.levels

by default unused levels are dropped from factors before fitting for gam() and bam()

drop.intercept

set to TRUE to force the model to really not have the a constant in the parametric model part for gam() and bam()

scale

scale parameter for bam() (should not used in gamlss)

chunk.size

the model matrix is created in chunks used in bam()

rho

an AR1 error model can be used for the residuals see the help in bam()

AR.start

see the help in bam() for explanation

discrete

with method="fREML" it is possible to discretize covariates see the help in bam()

cluster

whether QR decomposition in parallel see the help in bam()

nthreads

number of threads to use for non-cluster computation see the help in bam()

gc.level

to keep the memory footprint down see the help in bam()

use.chol

whether Choleski decomposition see the help in bam()

samfrac

sampling fraction see the help in bam()

...

extra options to pass to gam.control()

Details

Note that ga itself does no smoothing; it simply sets things up for the function gamlss() which in turn uses the function additive.fit() for back-fitting which in turn uses gamlss.ga()

Note that, in our (limited) experience, for normal errors or exponential family, the fitted models using gam() and ga() within gamlss() are identical or at least very similar. This is particularly true if the default values for gam() are used.

Value

the fitted values of the smoother is returned, endowed with a number of attributes. The smoother fitted values are used in the construction of the overall fitted values of the particular distribution parameter. The attributes can be use to obtain information about the individual fit. In particular the coefSmo within the parameters of the fitted model contains the final additive fit.

Warning

The function is experimental so please report any peculiar behaviour to the authors

Author(s)

Mikis Stasinopoulos

References

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. 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.

Wood S.N. (2006) Generalized Additive Models: An Introduction with R. Chapman and Hall/CRC Press.

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
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
library(mgcv)
data(rent)
#---------------------------------------------------------
## normal errors one x-variable
ga1 <- gam(R~s(Fl, bs="ps", k=20), data=rent, method="REML")
gn1 <- gamlss(R~ga(~s(Fl, bs="ps", k=20), method="REML"), data=rent) # additive
gb1 <- gamlss(R~pb(Fl), data=rent) # additive
AIC(ga1,gn1, gb1, k=0)
AIC(ga1,gn1, gb1)
#--------------------------------------------------------
## normal error additive in Fl and A
ga2 <- gam(R~s(Fl)+s(A), method="REML", data=rent)
gn2 <- gamlss(R~ga(~s(Fl)+s(A), method="REML"),  data=rent) # additive
gb2 <- gamlss(R~pb(Fl)+pb(A), data=rent) # additive
AIC(ga2,gn2, gb2, k=0)
AIC(ga2,gn2, gb2)
#---------------------------------------------------------
## Not run: 
## gamma error additive in Fl and A
ga3 <- gam(R~s(Fl)+s(A), method="REML", data=rent, family=Gamma(log))
gn3 <- gamlss(R~ga(~s(Fl)+s(A), method="REML"), data=rent, family=GA)# additive
gb3 <- gamlss(R~pb(Fl)+pb(A), data=rent, family=GA) # additive
AIC(ga3,gn3, gb3, k=0)
AIC(ga3,gn3, gb3)
#---------------------------------------------------------
## gamma error surface fitting
ga4 <-gam(R~s(Fl,A), method="REML", data=rent, family=Gamma(log))
gn4 <- gamlss(R~ga(~s(Fl,A), method="REML"), data=rent, family=GA) 
AIC(ga4,gn4, k=0)
AIC(ga4,gn4)
## plot the fitted surfaces
op<-par(mfrow=c(1,2))
vis.gam(ga4)
vis.gam(getSmo(gn4))
par(op)
## contour plot using mgcv's plot() function
plot(getSmo(gn4))
#---------------------------------------------------------
## predict
newrent <- data.frame(expand.grid(Fl=seq(30,120,5), A=seq(1890,1990,5 )))
newrent1 <-newrent2 <- newrent
newrent1$pred <- predict(ga4, newdata=newrent, type="response")
newrent2$pred <- predict(gn4, newdata=newrent, type="response")
library(lattice)
wf1<-wireframe(pred~Fl*A, newrent1, aspect=c(1,0.5), drape=TRUE, 
             colorkey=(list(space="right", height=0.6)), main="gam()")
wf2<-wireframe(pred~Fl*A, newrent2, aspect=c(1,0.5), drape=TRUE, 
          colorkey=(list(space="right", height=0.6)), main="gamlss()")
print(wf1, split=c(1,1,2,1), more=TRUE)
print(wf2, split=c(2,1,2,1))
#---------------------------------------------------------
##gamma error two variables te() function
ga5 <-  gam(R~te(Fl,A), data=rent, family=Gamma(log))
gn5 <- gamlss(R~ga(~te(Fl,A)), data=rent, family=GA) 
AIC(ga5,gn5)
AIC(ga5,gn5, k=0)
op<-par(mfrow=c(1,2))
vis.gam(ga5)
vis.gam(getSmo(gn5))
par(op)
#----------------------------------------------------------
## use of Markov random fields 
## example from package mgcv of Simon Wood
## Load Columbus Ohio crime data (see ?columbus for details and credits)
data(columb)       ## data frame
data(columb.polys) ## district shapes list
xt <- list(polys=columb.polys) ## neighbourhood structure info for MRF
## First a full rank MRF...
b <- gam(crime ~ s(district,bs="mrf",xt=xt),data=columb,method="REML")
bb <- gamlss(crime~ ga(~s(district,bs="mrf",xt=xt), method="REML"), data=columb)
AIC(b,bb, k=0)
op<-par(mfrow=c(2,2))
plot(b,scheme=1)
plot(bb$mu.coefSmo[[1]], scheme=1)
## Compare to reduced rank version...
b <- gam(crime ~ s(district,bs="mrf",k=20,xt=xt),data=columb,method="REML")
bb <- gamlss(crime~ ga(~s(district,bs="mrf",k=20,xt=xt), method="REML"), 
             data=columb)
AIC(b,bb, k=0)
plot(b,scheme=1)
plot(bb$mu.coefSmo[[1]], scheme=1)
par(op)
## An important covariate added...
b <- gam(crime ~ s(district,bs="mrf",k=20,xt=xt)+s(income),
         data=columb,method="REML")
## x in gam() 
bb <- gamlss(crime~ ga(~s(district,bs="mrf",k=20,xt=xt)+s(income), 
             method="REML"), data=columb)
## x in gamlss()
bbb <- gamlss(crime~ ga(~s(district,bs="mrf",k=20,xt=xt), 
             method="REML")+pb(income), data=columb)
AIC(b,bb,bbb)
## ploting the fitted models
op<-par(mfrow=c(2,2))
plot(b,scheme=c(0,1))
plot(getSmo(bb), scheme=c(0,1))
par(op)
plot(getSmo(bbb, which=2))
## plot fitted values by district
op<- par(mfrow=c(1,2))
fv <- fitted(b)
names(fv) <- as.character(columb$district)
fv1 <- fitted(bbb)
names(fv1) <- as.character(columb$district)
polys.plot(columb.polys,fv)
polys.plot(columb.polys,fv1)
par(op)
## End(Not run)
## bam