################# | |
# Use R to recreate part of QMV the relative non-proportional hazard estimation graph (fig 2) from Licht (2011) | |
# Christopher Gandrud (http://christophergandrud.blogspot.com/) | |
# 4 October 2012 | |
################# | |
#### Set Up #### | |
# Load required libraries | |
library(survival) | |
library(MSBVAR) | |
library(reshape) | |
library(reshape2) | |
library(ggplot2) | |
# Load Golub & Steunenberg (2007) Data | |
## Originally downloaded from http://hdl.handle.net/1902.1/15633 | |
GS <- read.table("GolubEUPdata.tab", header = TRUE) | |
# Create log time variable | |
GS$LT <- log(GS$end) | |
# Create natural log time interactions | |
attach(GS) | |
GS$Lqmv <- LT * qmv | |
GS$Lqmvpostsea <- LT * qmvpostsea | |
GS$Lcoop <- LT * coop | |
GS$Lcodec <- LT * codec | |
GS$Lthatcher <- LT * thatcher | |
GS$Lbacklog <- LT * backlog | |
detach(GS) | |
#### Run Cox PH Model #### | |
# Note, this model does not exactly match Licht (2011), but is pretty close | |
M1 <- coxph(Surv(begin, end, event) ~ | |
qmv + qmvpostsea + qmvpostteu + | |
coop + codec + eu9 + eu10 + eu12 + | |
eu15 + thatcher + agenda + backlog + | |
Lqmv + Lqmvpostsea + Lcoop + Lcodec + | |
Lthatcher + Lbacklog, | |
data = GS, | |
ties = "efron") | |
#### Simulate Relative Hazards #### | |
# Create coefficient matrix | |
Coef <- matrix(M1$coefficients) | |
# Create variance-covariance matrix | |
VC <- vcov(M1) | |
# Simulate Using MSBVAR | |
Drawn <- rmultnorm(n = 1000, mu = Coef, vmat = VC) | |
#### Create Combined Coefficient #### | |
# Keep qmv and Lqmv | |
Drawn <- data.frame(Drawn[, c(1, 13)]) | |
# Create Merge Variable (Simulation Number) | |
Drawn$ID <- 1:1000 | |
# Create range of years over which the combined coefficient will be created | |
Range <- log(seq(from = 80, to = 2000, by = 15)) | |
# Create Combined time interactionsCoefficient | |
TVSim <- outer(Drawn[,2], Range) | |
TVSim <- data.frame(melt(TVSim)) | |
TVSim <- rename(TVSim, c(X1 = "ID", X2 = "Time", value = "TVR")) | |
# Merge in the non-time interacted coefficient and combine | |
TVSim <- merge(Drawn, TVSim, by = "ID") | |
TVSim$CCqmv <- TVSim$qmv + TVSim$TVR | |
#### Create Combined Relative Hazard #### | |
TVSim$HRqmv <- exp(TVSim$CCqmv) | |
# Order Variables by time | |
TVSim <- TVSim[order(TVSim$Time),] | |
#### Graph Simulated Combined Hazard Ratios #### | |
ggplot(TVSim, aes(Time, HRqmv)) + | |
geom_point(shape = 21, alpha = I(0.01), colour = "#FA9FB5", size = 5) + | |
geom_smooth() + | |
geom_hline(aes(yintercept = 1), linetype = "dotted") + | |
scale_y_continuous(breaks = c(-1, 0, 1, 3, 6)) + | |
scale_x_continuous(breaks = c(0, 129), labels = c(80, 2000)) + | |
xlab("Time in Days") + ylab("Simulated QMV Relative Hazard\n") + | |
ggtitle("Simulated Relative Hazard for QMV from times 80-2000\n | |
Based on Licht (2011) Fig. 2\n") + | |
theme_bw(base_size = 15) |
This comment has been minimized.
Show comment
Hide comment
This comment has been minimized.
Show comment Hide comment
More InfoFor more information about this code please see this blog post |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
More Info
For more information about this code please see this blog post