gamlss: Generalised Additive Models for Location, Scale and Shape

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

View source: R/gamlss.r

Description

gamlss fits flexible univariate regression models with several continuous and discrete distributions, and types of covariate effects. The purpose of this function was only to provide, in some cases, starting values for the simultaneous models in the package, but it has now been made available in the form of a proper function should the user wish to fit univariate models using the general estimation approach of this package. The distributions implemented here have been parametrised according to Rigby and Stasinopoulos (2005).

Usage

1
2
3
4
5
6
gamlss(formula, data = list(), weights = NULL, subset = NULL,  
       margin = "N", surv = FALSE, cens = NULL,
       robust = FALSE, rc = 3, lB = NULL, uB = NULL, infl.fac = 1, 
       rinit = 1, rmax = 100, iterlimsp = 50, tolsp = 1e-07,
       gc.l = FALSE, parscale, extra.regI = "t", gev.par = -0.25,
       chunk.size = 10000)

Arguments

formula

List of equations. This should contain one or more equations.

data

An optional data frame, list or environment containing the variables in the model. If not found in data, the variables are taken from environment(formula), typically the environment from which gamlss is called.

weights

Optional vector of prior weights to be used in fitting.

subset

Optional vector specifying a subset of observations to be used in the fitting process.

margin

Possible distributions are normal ("N"), normal where sigma2 corresponds to the standard deviation instead of the variance ("N2"), log-normal ("LN"), Gumbel ("GU"), reverse Gumbel ("rGU"), logistic ("LO"), Weibull ("WEI"), inverse Gaussian ("iG"), gamma ("GA"), Dagum ("DAGUM"), Singh-Maddala ("SM"), beta ("BE"), Fisk ("FISK", also known as log-logistic distribution), Poisson ("PO"), zero truncated Poisson ("ZTP"), negative binomial - type I ("NBI"), negative binomial - type II ("NBII"), Poisson inverse Gaussian ("PIG"), generalised extreme value link function ("GEVlink", this is used for binary responses and is more stable and faster than the R package bgeva).

surv

If TRUE then a survival model is fitted. Here margin can be "PH" (generalised proportional hazards), "PO" (generalised proportional odds), "probit" (generalised probit).

cens

Binary censoring indicator. This is required when surv = TRUE.

robust

If TRUE then the robust version of the model is fitted.

rc

Robust constant.

lB, uB

Lower and upper bounds for integrals, when robust estimation is employed. If not provided then standard bounds are used.

infl.fac

Inflation factor for the model degrees of freedom in the approximate AIC. Smoother models can be obtained setting this parameter to a value greater than 1.

rinit

Starting trust region radius. The trust region radius is adjusted as the algorithm proceeds.

rmax

Maximum allowed trust region radius. This may be set very large. If set small, the algorithm traces a steepest descent path.

iterlimsp

A positive integer specifying the maximum number of loops to be performed before the smoothing parameter estimation step is terminated.

tolsp

Tolerance to use in judging convergence of the algorithm when automatic smoothing parameter estimation is used.

gc.l

This is relevant when working with big datasets. If TRUE then the garbage collector is called more often than it is usually done. This keeps the memory footprint down but it will slow down the routine.

parscale

The algorithm will operate as if optimizing objfun(x / parscale, ...) where parscale is a scalar. If missing then no rescaling is done. See the documentation of trust for more details.

extra.regI

If "t" then regularization as from trust is applied to the information matrix if needed. If different from "t" then extra regularization is applied via the options "pC" (pivoted Choleski - this will only work when the information matrix is semi-positive or positive definite) and "sED" (symmetric eigen-decomposition).

gev.par

GEV link parameter.

chunk.size

This is used for discrete robust models.

Details

The underlying algorithm is described in ?SemiParBIV.

There are many continuous/discrete distributions to choose from and we plan to include more options. Get in touch if you are interested in a particular distribution.

The "GEVlink" option is used for binary response additive models and is more stable and faster than the R package bgeva. This model has been incorporated into this package to take advantage of the richer set of smoother choices, and of the estimation approach. Details on the model can be found in Calabrese, Marra and Osmetti (2016).

Value

The function returns an object of class gamlss as described in gamlssObject.

WARNINGS

Convergence can be checked using conv.check which provides some information about the score and information matrix associated with the fitted model. The former should be close to 0 and the latter positive definite. gamlss() will produce some warnings if there is a convergence issue.

Convergence failure may sometimes occur. This is not necessarily a bad thing as it may indicate specific problems with a fitted model. In such a situation, the user may use some extra regularisation (see extra.regI) and/or rescaling (see parscale). However, the user should especially consider re-specifying/simplifying the model, and/or checking that the chosen distribution fits the response well. In our experience, we found that convergence failure typically occurs when the model has been misspecified and/or the sample size is low compared to the complexity of the model. It is also worth bearing in mind that the use of three parameter distributions requires the data to be more informative than a situation in which two parameter distributions are used instead.

Author(s)

Maintainer: Giampiero Marra [email protected]

References

Marra G. and Radice R. (2017), Bivariate Copula Additive Models for Location, Scale and Shape. Computational Statistics and Data Analysis, 112, 99-113.

Rigby R.A. and Stasinopoulos D.M. (2005). Generalized additive models for location, scale and shape (with discussion). Journal of the Royal Statistical Society, Series C, 54(3), 507-554.

Calabrese R., Marra G. and Osmetti SA (2016), Bankruptcy Prediction of Small and Medium Enterprises Using a Flexible Binary Generalized Extreme Value Model. Journal of the Operational Research Society, 67(4), 604-615.

See Also

JRM-package, gamlssObject, conv.check, summary.gamlss

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
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
## Not run:  

library(JRM)

set.seed(0)

n <- 400

x1 <- round(runif(n))
x2 <- runif(n)
x3 <- runif(n)
f1 <- function(x) cos(pi*2*x) + sin(pi*x)
y1 <- -1.55 + 2*x1 + f1(x2) + rnorm(n)

dataSim <- data.frame(y1, x1, x2, x3)
resp.check(y1, "N")

eq.mu <- y1 ~ x1 + s(x2) + s(x3)
eq.s2 <-    ~ s(x3)
fl    <- list(eq.mu, eq.s2)

out <- gamlss(fl, data = dataSim)

conv.check(out)
post.check(out)

plot(out, eq = 1, scale = 0, pages = 1, seWithMean = TRUE)
plot(out, eq = 2, seWithMean = TRUE)

summary(out)

AIC(out)
BIC(out)

################
# Robust example
################

eq.mu <- y1 ~ x1 + x2 + x3
fl    <- list(eq.mu)

out <- gamlss(fl, data = dataSim, margin = "N", robust = TRUE, 
                  rc = 3, lB = -Inf, uB = Inf)

conv.check(out)
summary(out)

##

eq.s2 <-    ~ x3
fl    <- list(eq.mu, eq.s2)

out <- gamlss(fl, data = dataSim, margin = "N", robust = TRUE)

conv.check(out)
summary(out)

##

eq.mu <- y1 ~ x1 + s(x2) + s(x3)
eq.s2 <-    ~ s(x3)
fl    <- list(eq.mu, eq.s2)

out1 <- gamlss(fl, data = dataSim, margin = "N", robust = TRUE)

conv.check(out1)
summary(out1)
AIC(out, out1)

plot(out1, eq = 1, all.terms = TRUE, pages = 1, seWithMean = TRUE)
plot(out1, eq = 2, seWithMean = TRUE)

##########################
## GEV link binary example
##########################
# this incorporates the bgeva
# model implemented in the bgeva package
# however this implementation is more general 
# stable and efficient

set.seed(0)

n <- 400

x1 <- round(runif(n)); x2 <- runif(n); x3 <- runif(n)

f1 <- function(x) cos(pi*2*x) + sin(pi*x)
f2 <- function(x) x+exp(-30*(x-0.5)^2)   

y  <- ifelse(-3.55 + 2*x1 + f1(x2) + rnorm(n) > 0, 1, 0)

dataSim <- data.frame(y, x1, x2, x3)

out1 <- gamlss(list(y ~ x1 + x2 + x3), margin = "GEVlink", data = dataSim)
out2 <- gamlss(list(y ~ x1 + s(x2) + s(x3)), margin = "GEVlink", data = dataSim)

conv.check(out1)
conv.check(out2)
summary(out1)
summary(out2)
AIC(out1, out2)
BIC(out1, out2)

plot(out2, eq = 1, all.terms = TRUE, pages = 1, seWithMean = TRUE)

##################
# prediction of Pr
##################

# Calculate eta (that is, X*model.coef)
# For a new data set the argument newdata should be used

eta <- predict(out2, eq = 1, type = "link")

# extract gev tail parameter

gev.par <- out2$gev.par

# multiply gev tail parameter by eta

gevpeta <- gev.par*eta 
  
# establish for which values the model is defined   

gevpetaIND <- ifelse(gevpeta < -1, FALSE, TRUE) 
gevpeta <- gevpeta[gevpetaIND]
    
# estimate probabilities  

pr <- exp(-(1 + gevpeta)^(-1/gev.par))

##################################
## Flexible survival model example
##################################

library(JRM)

########################################
## Simulate proportional hazards data ##
########################################

set.seed(0)
n  <- 2000
c  <- runif(n, 3, 8)
u  <- runif(n, 0, 1)
z1 <- rbinom(n, 1, 0.5)
z2 <- runif(n, 0, 1)
t <- rep(NA, n)

beta_0 <- -0.2357
beta_1 <- 1

f <- function(t, beta_0, beta_1, u, z1, z2){ 
  S_0 <- 0.7 * exp(-0.03*t^1.9) + 0.3*exp(-0.3*t^2.5)
  exp(-exp(log(-log(S_0))+beta_0*z1 + beta_1*z2))-u
}


for (i in 1:n){
   t[i] <- uniroot(f, c(0, 8), tol = .Machine$double.eps^0.5, 
                   beta_0 = beta_0, beta_1 = beta_1, u = u[i], 
                   z1 = z1[i], z2 = z2[i], extendInt = "yes" )$root
}

delta   <- ifelse(t < c, 1, 0)
u       <- apply(cbind(t, c), 1, min)
dataSim <- data.frame(u, delta, z1, z2)
1-mean(delta) # average censoring rate

out <- gamlss(list(u ~ z1 + s(z2) + s(u, bs = "mpi") ), data = dataSim, 
              surv = TRUE, margin = "PH", cens = delta)
post.check(out)
summary(out)
AIC(out)
BIC(out)
plot(out, eq = 1, scale = 0, pages = 1)
hazsurv.plot(out, newdata = data.frame(z1 = 0, z2 = 0), shade = TRUE, n.sim = 1000)
hazsurv.plot(out, type = "hazard", newdata = data.frame(z1 = 0, z2 = 0), 
             shade = TRUE, n.sim = 1000)

out1 <- gam(u ~ z1 + s(z2), family = cox.ph(), 
            data = dataSim, weights = delta)
summary(out1)
# estimates of z1 and s(z2) are
# nearly identical between out and out1 

# note that the Weibull is implemented as AFT
# as using the PH parametrisation makes
# computation unstable
out2 <- gamlss(list(u ~ z1 + s(z2) ), data = dataSim, surv = TRUE, 
               margin = "WEI", cens = delta)
 
#####################################
## Simulate proportional odds data ##
#####################################

set.seed(0)

n <- 2000
c <- runif(n, 4, 8)
u <- runif(n, 0, 1)
z <- rbinom(n, 1, 0.5)
beta_0 <- -1.05
t <- rep(NA, n)

f <- function(t, beta_0, u, z){ 
  S_0 <- 0.7 * exp(-0.03*t^1.9) + 0.3*exp(-0.3*t^2.5)
  1/(1 + exp(log((1-S_0)/S_0)+beta_0*z))-u
}



for (i in 1:n){
    t[i] <- uniroot(f, c(0, 8), tol = .Machine$double.eps^0.5, 
                    beta_0 = beta_0, u = u[i], z = z[i], 
                    extendInt="yes" )$root
}

delta   <- ifelse(t < c,1, 0)
u       <- apply(cbind(t, c), 1, min)
dataSim <- data.frame(u, delta, z)
1-mean(delta) # average censoring rate

out <- gamlss(list(u ~ z + s(u, bs = "mpi") ), data = dataSim, surv = TRUE, 
              margin = "PO", cens = delta)
post.check(out)
summary(out)
AIC(out)
BIC(out)
plot(out, eq = 1, scale = 0)
hazsurv.plot(out, newdata = data.frame(z =0), shade = TRUE, n.sim = 1000)
hazsurv.plot(out, type = "hazard", newdata = data.frame(z = 0), 
             shade = TRUE, n.sim = 1000)

# note that the Fisk is implemented as AFT
# as using the PH parametrisation makes
# computation unstable
out1 <- gamlss(list(u ~ z), data = dataSim, surv = TRUE, 
               margin = "FISK", cens = delta)


## End(Not run)

JRM documentation built on July 13, 2017, 5:03 p.m.