Fit Logistic-CoxPH Cure-Rate Model

This is a brief guide on package cureph illustrated with a simulated data set. First, you should load the package.

library(curephEM)

## Loading required package: survival

## Loading required package: Matrix

devtools::load_all()

## ℹ Loading curephEM

set.seed(531)

The simulated data

You may generate a simulated dataset using the build-in function cureph.simgen()

sim.cureph.data = cureph.simgen()

The true parameters are stored in the attributes true.coef and true.baseline.surv.

attr(sim.cureph.data, 'true.coef')

## $a
## [1]  1.00 -0.63  1.00  0.00  0.00  0.00  0.00
## 
## $b
## [1] -0.2  0.3  0.0  0.0  0.0  0.0

Fit the model

You may use the conventional syntax for coxph to fit a cureph, except for the use of a newly defined Surv.cure object instead of Surv.

If only one formula is provided, the set of covariates goes into both the logistic part and the cox part of the model.

fit=cureph(Surv.cure(time,time2,event,origin=0,end=20)~Z1+Z2+Z3+Z4,data=sim.cureph.data)

## Converge at step  55

Alternatively, you can provide two formulae—first formula for logistic part, then formula2 for cox part.

fit2=cureph(Surv.cure(time,time2,event,origin=0,end=20)~Z1+Z2+Z3+Z4,
  formula2 = ~ Z1+Z2,data=sim.cureph.data)

## Converge at step  62

Post Estimation

A detailed summary can be produced in a generic way. A multivariate Wald test table is activated if the two sets of covariates are detected to be the same. The null hypothesis is all the coefficients associated with the listed covariate are all zero.

summary(fit)

## Call:
## cureph(formula = Surv.cure(time, time2, event, origin = 0, end = 20) ~ 
##     Z1 + Z2 + Z3 + Z4, data = sim.cureph.data)
## 
## Logistic Model: 
## Surv.cure(time, time2, event, origin = 0, end = 20) ~ Z1 + Z2 + 
##     Z3 + Z4
## Cox Model: 
## Surv.cure(time, time2, event, origin = 0, end = 20) ~ Z1 + Z2 + 
##     Z3 + Z4
## 
##   n= 200, number of events= 32
##  
## Logistic:
##             coef    exp(coef) se(coef) z      Pr(>|z|)    
## (Intercept) 0.9459  2.5752    0.9336   1.013  0.311       
## Z1          -0.7521 0.4714    0.2282   -3.296 0.001    ***
## Z21         1.8655  6.4593    0.5319   3.507  5e-04    ***
## Z3B         -0.3145 0.7302    0.5646   -0.557 0.578       
## Z3C         0.214   1.2386    0.5917   0.362  0.718       
## Z41         -0.1346 0.874     0.516    -0.261 0.794       
## Z42         <NA>    <NA>      0        <NA>   <NA>        
## 
## Cox:
##     coef    exp(coef) se(coef) z      Pr(>|z|)  
## Z1  -0.1167 0.8899    0.3015   -0.387 0.699     
## Z21 0.7901  2.2037    0.5913   1.336  0.181     
## Z3B 0.2813  1.3248    0.5862   0.48   0.631     
## Z3C -0.132  0.8763    0.5723   -0.231 0.818     
## Z41 -0.0369 0.9637    0.5303   -0.07  0.944     
## Z42 <NA>    <NA>      0        <NA>   <NA>      
## 
## Combined Wald tests: 
##    Wald-chi.square df p-value    
## Z1 11.018          2  0.004   ** 
## Z2 15.46           2  4e-04   ***
## Z3 1.518           4  0.823      
## Z4 <NA>            4  <NA>       
## --- 
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Logistic:
##             exp(coef) exp(-coef) lower .95 upper .95
## (Intercept)    2.5752     0.3883    0.4132   16.0497
## Z1             0.4714     2.1214    0.3014    0.7372
## Z21            6.4593     0.1548    2.2773   18.3212
## Z3B            0.7302     1.3695    0.2414    2.2082
## Z3C            1.2386     0.8074    0.3884    3.9503
## Z41            0.8740     1.1441    0.3179    2.4031
## Z42                NA         NA        NA        NA
## 
## Cox:
##     exp(coef) exp(-coef) lower .95 upper .95
## Z1     0.8899     1.1237    0.4928     1.607
## Z21    2.2037     0.4538    0.6916     7.022
## Z3B    1.3248     0.7548    0.4199     4.179
## Z3C    0.8763     1.1411    0.2854     2.691
## Z41    0.9637     1.0376    0.3408     2.725
## Z42        NA         NA        NA        NA
## 
## Wald test = 43 on 11 df, p = 1.086e-05

Another function survpred combines conventional predict and survfit. It returns the linear predictors for both part, the estimated probabilities, the mean-baselevel survival in cox part and the marginal mean-baselevel survival.

mysurv = survpred(fit)

## Warning in prob0 * surv.cox: Recycling array of length 1 in array-vector arithmetic is deprecated.
##   Use c() or as.vector() instead.

## Warning in 1 - prob0 + prob0 * surv.cox: Recycling array of length 1 in array-vector arithmetic is deprecated.
##   Use c() or as.vector() instead.

The generic plot function produces the survival curves.

plot(mysurv)
plot(mysurv, pooled = F)

Left: Estimated marginal survival at mean-baselevel. Right: Baseline survival in cox part.