Created

Embed URL

HTTPS clone URL

SSH clone URL

You can clone with HTTPS or SSH.

Download Gist

Relative Non-proportional Hazard Graph in R

View relativeNPLGraph.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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Something went wrong with that request. Please try again.