cens: Function to Fit Censored Data Using a gamlss.family...

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

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

1
2
cens(family = "NO", type = c("right", "left", "interval"), name = "cens", 
       local = TRUE, delta = NULL, ...)

Arguments

family

a gamlss.family object, which is used to define the distribution and the link functions of the various parameters. The distribution families supported by gamlss() can be found in gamlss.family and in the package gamlss.dist.

name

the characters you want to add to the name of new functions, by default is cens

type

what type of censoring is required, right, left or interval.

local

if TRUE the function will try to find the environment of gamlss to generate the d and p functions required for the fitting, if FALSE the functions will be generated in the global environment

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.

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

cens.d, cens.p, gen.cens

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
# 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)
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)
summary(gg1)
#--------------------------------------------------------------------
# now fitting discretised continuous distribution to count data
# fitting discretised Gamma
data(species)
 mGA<-gamlss(Surv(fish,fish+1,type= "interval2")~log(lake)+I(log(lake)^2), sigma.fo=~log(lake), 
              data=species, family=cens(GA, type="interval"))
# 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=cens(IG, type="interval"))
AIC(mGA,mIG)
plot(fish~log(lake), data=species)
with(species, lines(log(lake)[order(lake)], fitted(mIG)[order(lake)]))             
#--------------------------------------------------------------------

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: survival
GAMLSS-RS iteration 1: Global Deviance = 194.3967 
GAMLSS-RS iteration 2: Global Deviance = 194.3963 
GAMLSS-RS iteration 3: Global Deviance = 194.3962 
GAMLSS-RS iteration 4: Global Deviance = 194.3961 
GAMLSS-RS iteration 5: Global Deviance = 194.3961 
GAMLSS-RS iteration 6: Global Deviance = 194.3961 
GAMLSS-RS iteration 7: Global Deviance = 194.3961 

Call:
survreg(formula = Surv(futime, fustat) ~ ecog.ps + rx, data = ovarian, 
    dist = "exponential")
             Value Std. Error     z       p
(Intercept)  6.962      1.322  5.27 1.4e-07
ecog.ps     -0.433      0.587 -0.74    0.46
rx           0.582      0.587  0.99    0.32

Scale fixed at 1 

Exponential distribution
Loglik(model)= -97.2   Loglik(intercept only)= -98
	Chisq= 1.67 on 2 degrees of freedom, p= 0.43 
Number of Newton-Raphson Iterations: 4 
n= 26 

******************************************************************
Family:  c("EXPrc", "right censored Exponential") 

Call:  
gamlss(formula = Surv(futime, fustat) ~ ecog.ps + rx, family = cens(EXP),  
    data = ovarian, c.crit = 1e-05) 

Fitting method: RS() 

------------------------------------------------------------------
Mu link function:  log
Mu Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   6.9608     0.8457   8.231 2.62e-08 ***
ecog.ps      -0.4324     0.3934  -1.099    0.283    
rx            0.5812     0.3922   1.482    0.152    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

------------------------------------------------------------------
No. of observations in the fit:  26 
Degrees of Freedom for the fit:  3
      Residual Deg. of Freedom:  23 
                      at cycle:  7 
 
Global Deviance:     194.3961 
            AIC:     200.3961 
            SBC:     204.1704 
******************************************************************
Warning message:
In summary.gamlss(gexp) :
  summary: vcov has failed, option qr is used instead

GAMLSS-RS iteration 1: Global Deviance = 194.6214 
GAMLSS-RS iteration 2: Global Deviance = 194.1882 
GAMLSS-RS iteration 3: Global Deviance = 194.1719 
GAMLSS-RS iteration 4: Global Deviance = 194.1707 
GAMLSS-RS iteration 5: Global Deviance = 194.17 
GAMLSS-RS iteration 6: Global Deviance = 194.1696 
GAMLSS-RS iteration 7: Global Deviance = 194.1693 
GAMLSS-RS iteration 8: Global Deviance = 194.1692 
GAMLSS-RS iteration 9: Global Deviance = 194.1691 
GAMLSS-RS iteration 10: Global Deviance = 194.1691 
GAMLSS-RS iteration 11: Global Deviance = 194.169 
GAMLSS-RS iteration 12: Global Deviance = 194.169 
GAMLSS-RS iteration 13: Global Deviance = 194.169 
GAMLSS-RS iteration 14: Global Deviance = 194.169 
GAMLSS-RS iteration 1: Global Deviance = 194.88 
GAMLSS-RS iteration 2: Global Deviance = 194.2411 
GAMLSS-RS iteration 3: Global Deviance = 194.1826 
GAMLSS-RS iteration 4: Global Deviance = 194.1731 
GAMLSS-RS iteration 5: Global Deviance = 194.1709 
GAMLSS-RS iteration 6: Global Deviance = 194.1702 
GAMLSS-RS iteration 1: Global Deviance = 192.6928 
GAMLSS-RS iteration 2: Global Deviance = 191.9227 
GAMLSS-RS iteration 3: Global Deviance = 191.8977 
GAMLSS-RS iteration 4: Global Deviance = 191.8955 
GAMLSS-RS iteration 5: Global Deviance = 191.8951 
(Intercept)     ecog.ps          rx 
  5.8784403  -0.2292750   0.8133591 
(Intercept)     ecog.ps          rx 
  5.8729054  -0.2259548   0.8104691 
  [1]  1-       1-       1-       1-      [11, 18]  1-       1-       1-     
  [9]  1-      [ 2,  4]  1-      [ 1,  2]  1-      [ 2,  4] 39+       1-     
 [17]  1-      [18, 25] 39+      39+       1-       1-      [ 1,  2] [ 4, 11]
 [25] [18, 25]  1-      [ 1,  2] [ 1,  2] [ 4, 11] [18, 25]  1-       1-     
 [33] [ 2,  4] 39+      39+      [ 2,  4] [ 2,  4] 39+      39+      39+     
 [41]  1-      [ 4, 11] [ 4, 11] 39+      39+       1-      [ 4, 11] [ 4, 11]
 [49] 39+      39+      [ 2,  4] 39+      39+      39+      39+      [ 4, 11]
 [57] 39+      39+      39+      39+       1-      [ 4, 11] 39+      39+     
 [65] 39+      [ 4, 11] 39+      [ 4, 11] 39+      39+      [ 4, 11] 39+     
 [73] 39+      39+      39+      39+      39+      39+      39+      39+     
 [81] [ 1,  2] 39+      39+      39+      39+      [ 4, 11] 39+      39+     
 [89] 39+      39+      [25, 32] 39+      39+      39+      39+      39+     
 [97] 39+      39+      39+      39+      [ 2,  4] 39+      39+      39+     
[105] 39+      [ 4, 11] 39+      39+      39+      39+      39+      39+     
[113] 39+      39+      39+      39+      39+      39+      39+      39+     

Call:
survreg(formula = y ~ poly(Tem, 2) + poly(pH, 2) + poly(aw, 2), 
    data = lip, dist = "weibull")
                Value Std. Error      z       p
(Intercept)     4.686      0.290  16.17 < 2e-16
poly(Tem, 2)1 -30.798      3.024 -10.19 < 2e-16
poly(Tem, 2)2   1.357      1.674   0.81   0.417
poly(pH, 2)1  -18.078      2.361  -7.66 1.9e-14
poly(pH, 2)2    3.894      1.685   2.31   0.021
poly(aw, 2)1  -27.596      2.760 -10.00 < 2e-16
poly(aw, 2)2   -1.072      1.554  -0.69   0.491
Log(scale)     -0.208      0.153  -1.36   0.173

Scale= 0.812 

Weibull distribution
Loglik(model)= -69   Loglik(intercept only)= -165.4
	Chisq= 192.77 on 6 degrees of freedom, p= 6.6e-39 
Number of Newton-Raphson Iterations: 13 
n= 120 

******************************************************************
Family:  c("WEIic", "interval censored Weibull") 

Call:  gamlss(formula = y ~ poly(Tem, 2) + poly(pH, 2) + poly(aw, 2),  
    family = cens(WEI, type = "interval"), data = lip, c.crit = 1e-05,  
    n.cyc = 200, trace = FALSE) 

Fitting method: RS() 

------------------------------------------------------------------
Mu link function:  log
Mu Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)     4.68997    0.07424  63.170  < 2e-16 ***
poly(Tem, 2)1 -30.83741    0.81330 -37.917  < 2e-16 ***
poly(Tem, 2)2   1.36136    0.81330   1.674   0.0969 .  
poly(pH, 2)1  -18.10342    0.81330 -22.259  < 2e-16 ***
poly(pH, 2)2    3.89797    0.81330   4.793 5.04e-06 ***
poly(aw, 2)1  -27.63020    0.81330 -33.973  < 2e-16 ***
poly(aw, 2)2   -1.07515    0.81330  -1.322   0.1888    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

------------------------------------------------------------------
Sigma link function:  log
Sigma Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   0.2067     0.0676   3.058  0.00276 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

------------------------------------------------------------------
No. of observations in the fit:  120 
Degrees of Freedom for the fit:  8
      Residual Deg. of Freedom:  112 
                      at cycle:  106 
 
Global Deviance:     138.0923 
            AIC:     154.0923 
            SBC:     176.3922 
******************************************************************
Warning message:
In summary.gamlss(gg1) :
  summary: vcov has failed, option qr is used instead

GAMLSS-RS iteration 1: Global Deviance = 604.3054 
GAMLSS-RS iteration 2: Global Deviance = 604.3006 
GAMLSS-RS iteration 3: Global Deviance = 604.3006 
GAMLSS-RS iteration 1: Global Deviance = 602.8633 
GAMLSS-RS iteration 2: Global Deviance = 602.8628 
    df      AIC
mIG  5 612.8628
mGA  5 614.3006

gamlss.cens documentation built on May 2, 2019, 6:49 p.m.