Middle95HR.R

#################
# Use R to recreate part of QMV the relative non-proportional hazard estimation graph (fig 2) from Licht (2011) 
# Updated to keep only the middle 95% of simulations
# Christopher Gandrud (http://christophergandrud.blogspot.com/)
# 30 December 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),]

# Drop bottom 2.5% of observations
TVSimPerc <- ddply(TVSim, .(Time), transform, Lower = HRqmv < quantile(HRqmv, c(0.025)))

# Drop top 2.5% of observations
TVSimPerc <- ddply(TVSimPerc, .(Time), transform, Upper = HRqmv > quantile(HRqmv, c(0.975)))

# Remove variables outside of the middle 95%
TVSimPerc <- subset(TVSimPerc, Lower == FALSE & Upper == FALSE)

#### Graph Simulated Combined Hazard Ratios ####
ggplot(TVSimPerc, 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, middle 95%\n") +
  theme_bw(base_size = 15)