#################
# 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)