# NOT RUN {
library(missDeaths)
ratetable = survexp.us
sim = md.simparams() +
md.sex("sex", 0.5) +
md.binom("Z1", 0.5) +
md.mvnorm(c("Z2", "Z3"), c(100, 0), matrix(c(225, 3, 2, 1), 2, 2)) +
md.sample("Z4", c(1, 2, 3, 4), c(0.25, 0.25, 0.25, 0.25)) +
md.uniform("birth", as.Date("1925-1-1"), as.Date("1950-1-1")) +
md.uniform("start", as.Date("2000-1-1"), as.Date("2005-1-1")) +
md.death("death", ratetable, "sex", "birth", "start") +
md.eval("age", "as.numeric(start - birth)/365.2425", 80, FALSE) +
md.exp("event", "start", c("age", "sex", "Z1", "Z2"),
c(0.1, 2, 1, 0.01), 0.05/365.2425)
data = md.simulate(sim, 1000)
#construct a complete right censored data set
complete = md.censor(data, "start", as.Date("2010-1-1"), c("event", "death"))
#construct an observed right censored data set with non-reported deaths
observed = complete
observed$time[which(complete$status == 2)] = observed$maxtime[which(complete$status == 2)]
observed$status[which(complete$status == 2)] = 0
#estimate coefficients of the proportional hazards model
D = md.D(age=observed$age, sex=observed$sex, year=observed$start)
md.survcox(observed, Surv(time, status) ~ age + sex + Z1 + Z2,
observed$maxtime, D, ratetable, iterations=4, R=50)
#estimate net- and event-free survival
np = md.survnp(observed$time, observed$status, observed$maxtime, D, ratetable)
w = list(list(time=np$time, est=np$surv.net, var=(np$std.err.net)^2))
timepoints(w, times=c(3,9)*365.2425)
plot(np$time/365.2425, np$surv.net)
plot(np$time/365.2425, np$surv.efs)
# }
Run the code above in your browser using DataCamp Workspace