# Load Carpenter (2002) data
data("CarpenterFdaData")
# Load survival package
library(survival)
# Run basic model
# From Keele (2010) replication data
M1 <- coxph(Surv(acttime, censor) ~ prevgenx + lethal + deathrt1 +
acutediz + hosp01 + pspline(hospdisc, df = 4) +
pspline(hhosleng, df = 4) + mandiz01 + femdiz01 +
peddiz01 + orphdum + natreg + vandavg3 + wpnoavg3 +
pspline(condavg3, df = 4) + pspline(orderent, df = 4) +
pspline(stafcder, df = 4), data = CarpenterFdaData)
# Simulate Relative Hazards for orderent
Sim1 <- coxsimSpline(M1, bspline = "pspline(stafcder, df = 4)",
bdata = CarpenterFdaData$stafcder,
qi = "Hazard Ratio",
Xj = seq(1100, 1700, by = 10),
Xl = seq(1099, 1699, by = 10), spin = TRUE, nsim = 100)
# Plot relative hazard
simGG(Sim1, alpha = 0.5)
## Not run:
# # Simulate Hazard Rate for orderent
# Sim2 <- coxsimSpline(M1, bspline = "pspline(orderent, df = 4)",
# bdata = CarpenterFdaData$orderent,
# qi = "Hazard Rate",
# Xj = seq(1, 30, by = 2), ci = 0.9, nsim = 10)
#
# # Create a time grid plot
# # Find all points in time where baseline hazard was found
# unique(Sim2$sims$Time)
#
# # Round time values so they can be exactly matched with FacetTime
# Sim2$sims$Time <- round(Sim2$sims$Time, digits = 2)
#
# # Create plot
# simGG(Sim2, FacetTime = c(6.21, 25.68, 100.64, 202.36),
# type = 'ribbons', alpha = 0.5)
#
# # Simulated Fitted Values of stafcder
# Sim3 <- coxsimSpline(M1, bspline = "pspline(stafcder, df = 4)",
# bdata = CarpenterFdaData$stafcder,
# qi = "Hazard Ratio",
# Xj = seq(1100, 1700, by = 10),
# Xl = seq(1099, 1699, by = 10), ci = 0.90)
#
# # Plot simulated Hazard Ratios
# simGG(Sim3, xlab = "\nFDA Drug Review Staff", type = 'lines', alpha = 0.2)
# simGG(Sim3, xlab = "\nFDA Drug Review Staff", alpha = 0.2,
# SmoothSpline = TRUE, type = 'points')
# ## End(Not run)
Run the code above in your browser using DataLab