| cens {gamlss.cens} | R Documentation |
Function to Fit Censored Data Using a gamlss.family Distribution
Description
This function can be used to fit censored or interval response variables.
It takes as an argument an existing gamlss.family distribution
and generates
a new gamlss.family object which then can be used to fit
right, left or interval censored data.
Usage
cens(family = "NO", type = c("right", "left", "interval"), name = "cens",
local = TRUE, delta = NULL, ...)
Arguments
family |
a |
name |
the characters you want to add to the name of new functions, by default is |
type |
what type of censoring is required, |
local |
if TRUE the function will try to find the environment of |
delta |
the delta increment used in the numerical derivatives |
... |
for extra arguments |
Details
This function is created to help users to fit censored data using an existing
gamlss.family distribution.
It does this by taking an existing gamlss.family and changing
some of the components of the distribution to help the fitting process.
It particular it (i) creates a (d) function (for calculating the censored
likelihood) and a (p) function (for generating the quantile residuals)
within gamlss,
(ii) changes the global deviance function G.dev.incr,
the first derivative functions (see note below)
and other quantities from the original distribution.
Value
It returns a gamlss.family object which has all the components needed for fitting a distribution in gamlss.
Note
This function is experimental and could be changed in the future.
The function cens changes the first derivatives of the original gamlss family
d function to numerical derivatives for the new censored d function.
The default increment delta, for this numerical derivatives function,
is eps * pmax(abs(x), 1) where eps<-sqrt(.Machine$double.eps).
The default delta could be inappropriate for
specific applications and can be overwritten by using the argument delta.
Note that in order to get the correct standard errors you have to generate the "d" function by using
gen.cens().
Author(s)
Mikis Stasinopoulos d.stasinopoulos@londonmet.ac.uk and Bob Rigby r.rigby@londonmet.ac.uk
References
Rigby, R. A. and Stasinopoulos D. M. (2005). Generalized additive models for location, scale and shape,(with discussion), Appl. Statist., 54, part 3, 1-38.
Rigby, R. A., Stasinopoulos, D. M., Heller, G. Z., and De Bastiani, F. (2019) Distributions for modeling location, scale, and shape: Using GAMLSS in R, Chapman and Hall/CRC. An older version can be found in https://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, https://www.jstatsoft.org/v23/i07/.
Stasinopoulos D. M., Rigby R.A., Heller G., Voudouris V., and De Bastiani F., (2017) Flexible Regression and Smoothing: Using GAMLSS in R, Chapman and Hall/CRC.
(see also https://www.gamlss.com/).
See Also
Examples
# comparing output with the survreg() of package survival
library(gamlss.dist)
library(survival)
#--------------------------------------------------------------------
# right censoring example
# example from survreg()
# fitting the exponential distribution
mexp<-survreg(Surv(futime, fustat) ~ ecog.ps + rx, ovarian, dist='exponential')
gexp<-gamlss(Surv(futime, fustat) ~ ecog.ps + rx, data=ovarian,
family=cens(EXP), c.crit=0.00001)
if(abs(-2*mexp$loglik[2]-deviance(gexp))>0.001)
stop(paste("descrepancies in exponential models"))
if(sum(coef(mexp)-coef(gexp))>0.001)
warning(paste("descrepancies in coef in exponential models"))
summary(mexp)
gen.cens(EXP)
summary(gexp)
# fitting different distributions
# weibull
mwei <-survreg(Surv(futime, fustat) ~ ecog.ps + rx, ovarian, dist='weibull')
gwei<-gamlss(Surv(futime, fustat) ~ ecog.ps + rx, data=ovarian,
family=cens(WEI, delta=c(0.0001,0.0001)), c.crit=0.00001)
if(abs(-2*mwei$loglik[2]-deviance(gwei))>0.005)
stop(paste("descrepancies in deviance in WEI"))
scoef <- sum(coef(mwei)-coef(gwei))
if(abs(scoef)>0.005)
warning(cat("descrepancies in coef in WEI of ", scoef, "\n"))
# WEI3 is weibull parametrised with mu as the mean
gwei3 <- gamlss(Surv(futime, fustat) ~ ecog.ps + rx, data=ovarian,
family=cens(WEI3))
# log normal
mlogno <-survreg(Surv(futime, fustat) ~ ecog.ps + rx, ovarian,
dist='lognormal')
glogno<-gamlss(Surv(futime, fustat) ~ ecog.ps + rx, data=ovarian,
family=cens(LOGNO, delta=c(0.001,0.001)), c.cyc=0.00001)
if(abs(-2*mlogno$loglik[2]-deviance(glogno))>0.005)
stop(paste("descrepancies in deviance in LOGNO"))
coef(mlogno);coef(glogno)
#--------------------------------------------------------------------
# now interval response variable
data(lip)
with(lip, y)
mg1<-survreg(y ~ poly(Tem,2)+poly(pH,2)+poly(aw,2), data=lip, dist="weibull")
gg1<- gamlss(y ~ poly(Tem,2)+poly(pH,2)+poly(aw,2), data=lip,
family=cens(WEI,type="interval"), c.crit=0.00001, n.cyc=200, trace=FALSE)
summary(mg1)
gen.cens(WEI,type="interval")
summary(gg1)
#--------------------------------------------------------------------
# now fitting discretised continuous distribution to count data
# fitting discretised Gamma
data(species)
# first generate the distributions
gen.cens(GA, type="interval")
gen.cens(IG, type="interval")
mGA<-gamlss(Surv(fish,fish+1,type= "interval2")~log(lake)+I(log(lake)^2),
sigma.fo=~log(lake), data=species, family=GAic)
# fitting discretised inverse Gaussian
mIG<-gamlss(Surv(fish,fish+1,type= "interval2")~log(lake)+I(log(lake)^2),
sigma.fo=~log(lake), data=species, family=IGic)
AIC(mGA,mIG)
plot(fish~log(lake), data=species)
with(species, lines(log(lake)[order(lake)], fitted(mIG)[order(lake)]))
#--------------------------------------------------------------------