0
$\begingroup$

I have done the following bit myself:

#Prediction using Mixture Cure Modelling
pd<-smcure(Surv(Line_Tenure_In_Days,Churn_Flag) ~ Months_On_Device,cureform = ~ Months_On_Device,data=df_acn4,model = "ph",nboot=100,Var=TRUE)
printsmcure(pd)

predm <- predictsmcure(pd,newX=c(0,1),newZ=c(0,1),model = "ph")
predm_predict<-predm$prediction
as.data.frame(predm_predict)
uncure_prob<-predm$newuncureprob
as.data.frame(uncure_prob)
plotpredictsmcure(predm,model="ph")

Now my aim is to obtain the individual survival probabilities at a fixed time instance using Mixture Cure model for which I am using this smcure() in R.

Questions: I am well aware of the mixture cure model see(https://deepblue.lib.umich.edu/bitstream/handle/2027.42/65901/j.0006-341X.2000.00227.x.pdf?sequence=1&isAllowed=y)

titled Estimation in a Cox Proportional Hazards Cure Model

I am well aware of the smcure() in R.

I also know the Breslow type estimator to calculate $ \Lambda_0(t|Y=1)=\Sigma_{i:t_i\le t}{\dfrac{d_i}{\Sigma_{l \in R_i}w_l^{m}e^{z_l^T.\beta}}}$

And finally $S_{pop}(t∣X, Z) = π(z)S(t∣X) + 1 - π(z)$

Now how do I obtain $\Lambda_0(t|Y=1)$ in R for the mixture cure model?

How do I obtain $S_{pop}(t∣X, Z) = π(z)S(t∣X) + 1 - π(z)$ for individual observations using the smcure() in R?(By this I mean the following: Is there an equivalent of predict(smcureobject,type="lp") that which we used to do find out the linear predictors in Cox PH model for individual observations?)

I give you how I obtained the individual survival probabilities using CoxPh in R

cox_Jan_Aug <-coxph(Surv(Line_Tenure_In_Days,Churn_Flag)~Months_On_Device+Subscriber_Activity_Price_Plan_Code+Months_On_Price_Plan+Price_Plan_Change_Flag+Total_MOU+Total_Calls+Market_Name+BYOD_Indicator1,method = "breslow")
summary(cox_Jan_Aug)
pred_Jan_Aug <- survfit(cox_Jan_Aug)
plot(pred_Jan_Aug,fun=function(y)log(y/(1-y)), ylab="Logit S(t)",col="Red")

plot(pred_Jan_Aug,
     main = 'Cox Regression Plot',col="Blue",width=1,height=1)
abline(h = 0.5,col="Red")
summary(pred_Jan_Aug)
pred_coef_times_vars4<- predict(cox_Jan_Aug,type = "lp")
basehaz_Jan_Aug <- basehaz(cox_Jan_Aug)

predicted_on_15thday4 <- as.data.frame((exp(-basehaz_Jan_Aug[16,1]))^(exp(pred_coef_times_vars4)))
colnames(predicted_on_15thday4)<- 'predicted survival probability on 15th day'

predicted_on_20thday4 <- as.data.frame((exp(-basehaz_Jan_Aug[21,1]))^(exp(pred_coef_times_vars4)))
colnames(predicted_on_20thday4)<- 'predicted survival probability on 20th day'

predicted_on_30thday4 <- as.data.frame((exp(-basehaz_Jan_Aug[31,1]))^(exp(pred_coef_times_vars4)))
colnames(predicted_on_30thday4)<- 'predicted survival probability on 30th day'

Indiv_Prob_Jan_Aug <- cbind(predicted_on_15thday4,predicted_on_20thday4$`predicted survival probability on 20th day`,predicted_on_30thday4$`predicted survival probability on 30th day`,Line_Tenure_In_Days,Total_MOU,Churn_Flag,Months_On_Device,Months_On_Price_Plan,Price_Plan_Change_Flag,Total_Calls,BYOD_Indicator1,Subscriber_Activity_Price_Plan_Code)
write.csv(Indiv_Prob_Jan_Aug,file = "Individual Probabilities on Jan through August.csv")

A few words:

  • Please try not to downvote this post for not posting the original data here
  • Please try to use any standard data like e1684 dataset.
$\endgroup$
2
  • 1
    $\begingroup$ smsurv looks like what you're after. NB: the "mixture cure model" is an abuse of terminology. Farwell and Sy appropriately referred to such models as mixture or cure models: Cox models where an unobserved proportion of individuals are not at risk for the event. NB also: "mixture models" have come to take a different meaning over time. If the smcure function gives you some posterior probability of being not-at-risk, you can classify based on thresholds, fit the normal Cox model only among those at risk, and use the basehaz function. $\endgroup$
    – AdamO
    Oct 18, 2018 at 14:00
  • 1
    $\begingroup$ instructing people to not downvote you seems like shady site etiquette. Are you not able to makea reproducible example yourself with "standard data like the e1684 dataset"? $\endgroup$
    – AdamO
    Oct 18, 2018 at 14:04

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.