Created
October 4, 2012 04:28
-
-
Save christophergandrud/3831482 to your computer and use it in GitHub Desktop.
Relative Non-proportional Hazard Graph in R
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
################# | |
# 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) |
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