library(copula) library(survival) library(splines2) library(survivalMPLdc) ##-- Copula types copula1 <- 'clayton' copula2 <- 'gumbel' copula3 <- 'frank' copula4 <- 'independent' ##-- Marginal distribution for T, C, and A a <- 2 lambda <- 2 cons7 <- 0.2 cons9 <- 10 tau <- 0.8 betas <- c(-0.5, 0.1) phis <- c(0.3, 0.2) distr.ev <- 'weibull' distr.ce <- 'exponential' ##-- Sample size n <- 200 ##-- One sample Monte Carlo dataset cova <- cbind(rbinom(n, 1, 0.5), runif(n, min=-10, max=10)) surv <- surv_data_dc(n, a, cova, lambda, betas, phis, cons7, cons9, tau, copula3, distr.ev, distr.ce) n <- nrow(cova) p <- ncol(cova) ##-- event and dependent censoring proportions colSums(surv)[c(2,3)]/n X <- surv[,1] # Observed time del<-surv[,2] # failure status eta<-surv[,3] # dependent censoring status ##-- Selecting bin sample or number of knots using AIC vs binCounts plot binCounts <- seq(10, 100, 10) bn <- length(binCounts) aics<-rep(0, bn) for (j in 1:bn) { control=coxph_mpl_dc.control(ordSp = 4, binCount = binCounts[j], tie = 'No', tau = 0.8, copula = copula3, pent = 'penalty_mspl', smpart = 0, penc = 'penalty_mspl', smparc = 0, maxit2 = 100, maxit = 5000, mid = 1, asy = 0, ac = 1, cv = 0, ac.theta = 1e-5, ac.gamma = 1e-5, ac.Utheta = -1e-2, ac.Ugamma = -1e-2, min.theta = 1e-7, min.gamma = 1e-7, min.ht = 1e-7, min.hc = 1e-7, min.St = 1e-7, min.Sc = 1e-7, min.C = 1e-7, min.dC = 1e-7, eps = 1e-5, tol.thga = 1e-5, tol.bph = 1e-5, tol.smpar = 1e-2 ) aics[j]<-coxph_mpl_dc(surv, cova, control)$mpl_aic print(j) } binCount <- 10 + ( which.min( aics ) - 1 ) * 10 binCount plot(binCounts, aics) control <- coxph_mpl_dc.control(ordSp = 4, binCount = binCount, tie = 'No', tau = 0.8, copula = copula3, pent = 'penalty_mspl', smpart = 'REML', penc = 'penalty_mspl', smparc = 'REML', maxit2 = 100, maxit = 50000, mid = 1, asy = 1, ac = 1, cv = 1, ac.theta = 1e-5, ac.gamma = 1e-5, ac.Utheta = -1e-2, ac.Ugamma = -1e-2, min.theta = 1e-7, min.gamma = 1e-7, min.ht = 1e-7, min.hc = 1e-7, min.St = 1e-7, min.Sc = 1e-7, min.C = 1e-7, min.dC = 1e-7, eps = 1e-5, tol.thga = 1e-5, tol.bph = 1e-5, tol.smpar = 1e-2 ) control coxMPLests5 <- coxph_mpl_dc(surv=surv, cova=cova, control=control, ) mpl_theta5 <- coxMPLests5$mpl_theta mpl_gamma5 <- coxMPLests5$mpl_gamma mpl_beta5 <- coxMPLests5$mpl_beta mpl_phi5 <- coxMPLests5$mpl_phi mpl_beta_sd5 <- coxMPLests5$mpl_beta_sd mpl_phi_sd5 <- coxMPLests5$mpl_phi_sd mpl_h0t_sd5 <- coxMPLests5$mpl_h0t_sd mpl_h0c_sd5 <- coxMPLests5$mpl_h0c_sd mpl_beta_phi_zp5 <- coxMPLests5$mpl_beta_phi_zp mpl_h0t5 <- coxMPLests5$mpl_h0t mpl_beta_phi_zp5 coxMPLests5$smootht coxMPLests5$smoothc coxMPLests5$iter mpl_h0Ti5 <- approx( X, mpl_h0t5, xout = seq(0, 5.4, 0.01), method="constant", rule = 2, ties = mean)$y ##-- Real marginal baseline hazard for T ht0b <- a * (seq(0, 5.4, 0.01) ^ (a - 1)) / (lambda ^ a) ##-- Selectiong a smoothing parameter for T using CV method #given zero for the smoothing parameter of C smparts <- seq(10, 100, 10) smn <- length(smparts) smparcs <- rep(0, smn) cvls<-rep(0, smn) for (j in 1:smn){ control <- coxph_mpl_dc.control(ordSp = 4, binCount = binCount, tie = 'No', tau = 0.8, copula = copula3, pent = 'penalty_mspl', smpart = smparts[j], penc = 'penalty_mspl', smparc = smparcs[j], maxit2 = 100, maxit = 10000, mid = 1, asy = 0, ac = 0, cv = 1, ac.theta=1e-5, ac.gamma=1e-5, ac.Utheta=-1e-2, ac.Ugamma=-1e-2, min.theta=1e-7, min.gamma=1e-7, min.ht=1e-7, min.hc=1e-7, min.St=1e-7, min.Sc=1e-7, min.C=1e-7, min.dC=1e-7, eps=1e-5, tol.thga=1e-5, tol.bph=1e-5, tol.smpar=1e-2) cvls[j]<-coxph_mpl_dc(surv, cova, control)$mpl_cvl print(j) } cvls smpart <- 10 + ( which.max( cvls ) - 1 ) * 10 smpart plot(smparts, cvls) coxMPLests4 <- coxph_mpl_dc(surv, cova, ordSp=4, binCount=binCount, tie='No', tau=0.8, copula=copula3, pent='penalty_mspl', smpart=smpart, penc='penalty_mspl', smparc=0, maxit2=100, maxit=100000, mid=1, asy=1, ac=1, cv=1, ac.theta=1e-5, ac.gamma=1e-5, ac.Utheta=-1e-2, ac.Ugamma=-1e-2, min.theta=1e-7, min.gamma=1e-7, min.ht=1e-7, min.hc=1e-7, min.St=1e-7, min.Sc=1e-7, min.C=1e-7, min.dC=1e-7, eps=1e-5, tol.thga=1e-5, tol.bph=1e-5, tol.smpar=1e-2) mpl_theta4 <- coxMPLests4$mpl_theta mpl_gamma4 <- coxMPLests4$mpl_gamma mpl_beta4 <- coxMPLests4$mpl_beta mpl_phi4 <- coxMPLests4$mpl_phi mpl_beta_sd4 <- coxMPLests4$mpl_beta_sd mpl_phi_sd4 <- coxMPLests4$mpl_phi_sd mpl_h0t_sd4 <- coxMPLests4$mpl_h0t_sd mpl_h0c_sd4 <- coxMPLests4$mpl_h0c_sd mpl_beta_phi_zp4 <- coxMPLests4$mpl_beta_phi_zp mpl_beta_phi_zp4 mpl_h0t4 <- coxMPLests4$mpl_h0t coxMPLests4$smootht coxMPLests4$smoothc coxMPLests4$iter mpl_h0Ti4 <- approx( X, mpl_h0t4, xout = seq(0, 5.4, 0.01), method="constant", rule = 2, ties = mean)$y ##-- Fitting cox ph hazard model for T using MPL and an correct copula #with zero smoothing parameters coxMPLests3 <- coxph_mpl_dc(surv, cova, ordSp=4, binCount=binCount, tie='No', tau=0.8, copula=copula3, pent='penalty_mspl', smpart=0, penc='penalty_mspl', smparc=0, maxit2=100, maxit=100000, mid=1, asy=1, ac=1, cv=1, ac.theta=1e-5, ac.gamma=1e-5, ac.Utheta=-1e-2, ac.Ugamma=-1e-2, min.theta=1e-7, min.gamma=1e-7, min.ht=1e-7, min.hc=1e-7, min.St=1e-7, min.Sc=1e-7, min.C=1e-7, min.dC=1e-7, eps=1e-5, tol.thga=1e-5, tol.bph=1e-5, tol.smpar=1e-2) mpl_theta3 <- coxMPLests3$mpl_theta mpl_gamma3 <- coxMPLests3$mpl_gamma mpl_beta3 <- coxMPLests3$mpl_beta mpl_phi3 <- coxMPLests3$mpl_phi mpl_beta_sd3 <- coxMPLests3$mpl_beta_sd mpl_phi_sd3 <- coxMPLests3$mpl_phi_sd mpl_h0t_sd3 <- coxMPLests3$mpl_h0t_sd mpl_h0c_sd3 <- coxMPLests3$mpl_h0c_sd mpl_beta_phi_zp3 <- coxMPLests3$mpl_beta_phi_zp mpl_beta_phi_zp3 mpl_h0t3 <- coxMPLests3$mpl_h0t coxMPLests3$smootht coxMPLests3$smoothc coxMPLests3$iter mpl_h0Ti3 <- approx( X, mpl_h0t3, xout = seq(0, 5.4, 0.01), method="constant", rule = 2, ties = mean)$y ##-- Plot the true and estimated baseline hazards for T t_up <- 3.5 y_uplim <- 2 plot(seq(0, 5.4, 0.01)[seq(0, 5.4, 0.01)<=t_up], mpl_h0Ti5[seq(0, 5.4, 0.01)<=t_up], type="l", col="grey", lty=4, lwd=3, cex.axis=1.6, cex.lab=1.6, ylim=c(0, y_uplim), xlab='Time', ylab='Hazard') par(new=TRUE) plot(seq(0, 5.4, 0.01)[seq(0, 5.4, 0.01)<=t_up], ht0b[seq(0, 5.4, 0.01)<=t_up], type="l", col="green", lty=1, lwd=3, cex.axis=1.6, cex.lab=1.6, ylim=c(0, y_uplim), xlab='Time', ylab='Hazard') par(new=TRUE) plot(seq(0, 5.4, 0.01)[seq(0, 5.4, 0.01)<=t_up], mpl_h0Ti4[seq(0, 5.4, 0.01)<=t_up], type="l", col="red", lty=4, lwd=3, cex.axis=1.6, cex.lab=1.6, ylim=c(0, y_uplim), xlab='Time', ylab='Hazard') par(new=TRUE) plot(seq(0, 5.4, 0.01)[seq(0, 5.4, 0.01)<=t_up], mpl_h0Ti3[seq(0, 5.4, 0.01)<=t_up], type="l", col="blue", lty=4, lwd=3, cex.axis=1.6, cex.lab=1.6, ylim=c(0, y_uplim), xlab='Time', ylab='Hazard') dev.off()