# NOT RUN {
##-- Copula types
copula3 <- 'frank'
##-- 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
##-- control inputs for the coxph_mpl_dc function
control <- coxph_mpl_dc.control(ordSp = 4,
binCount = 100,
tau = 0.8, copula = copula3,
pent = 'penalty_mspl', smpart = 'REML',
penc = 'penalty_mspl', smparc = 'REML',
cat.smpar = 'No' )
##-- Fitting cox ph hazard model for T using MPL and an correct copula
#with REML smoothing parameters
coxMPLests5 <- coxph_mpl_dc(surv, cova, control, )
mpl_beta_phi_zp5 <- coxMPLests5$mpl_beta_phi_zp
mpl_h0t5 <- coxMPLests5$mpl_h0t
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)
##-- 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=100,
tau=0.8, copula=copula3,
pent='penalty_mspl', smpart=0, penc='penalty_mspl', smparc=0,
cat.smpar = 'No')
mpl_beta_phi_zp3 <- coxMPLests3$mpl_beta_phi_zp
mpl_h0t3 <- coxMPLests3$mpl_h0t
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
Ti<-seq(0, 5.4, 0.01)[seq(0, 5.4, 0.01)<=t_up]
h0Ti<-ht0b[seq(0, 5.4, 0.01)<=t_up]
h0Ti5<-mpl_h0Ti5[seq(0, 5.4, 0.01)<=t_up]
h0Ti3<-mpl_h0Ti3[seq(0, 5.4, 0.01)<=t_up]
plot( x = Ti, y = h0Ti5,
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')
lines(x = Ti, y = h0Ti,
col="green",
lty=1, lwd=3, cex.axis=1.6, cex.lab=1.6, ylim=c(0, y_uplim)
)
lines(x = Ti, y = h0Ti3,
col="blue",
lty=4, lwd=3, cex.axis=1.6, cex.lab=1.6, ylim=c(0, y_uplim)
)
# }
# NOT RUN {
# }
Run the code above in your browser using DataLab