public
Created

Relative Non-proportional Hazard Graph in R

  • Download Gist
relativeNPLGraph.R
R
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92
#################
# 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)

More Info

For more information about this code please see this blog post

Please sign in to comment on this gist.

Something went wrong with that request. Please try again.