Jeff Compare2

JeffNPLHazard2.R

###############
# NPL Relative Hazard 2
# Christopher Gandrud
# 6 October 2012
###############

# Load required libraries
library(foreign)
library(survival)
library(MSBVAR)
library(reshape)
library(reshape2)
library(ggplot2)

# Load Data
Data1 <- read.dta("PS_Data_Current.dta")

#### Clean up base data set ####
# Keep only variables used in the model
vars <- c("begin", "end", "Partisan_Spell", "BEGINbc", "executive_dom_fall",
          "cumulative_crisis", "polity2", "Growth", "periphery_LAm", "imfcode")

# Keep complete cases
CompleteData <- Data1[complete.cases(Data1[vars]),]
CompleteData <- CompleteData[vars]


# Create  interactions
CompleteData$ED_TVC <- CompleteData$executive_dom_fall * CompleteData$end

# Run Cox PH Model
M1 <- coxph(Surv(begin, end, Partisan_Spell) ~ 
                   ED_TVC + BEGINbc*executive_dom_fall + 
                   cumulative_crisis + polity2 + Growth + 
                   periphery_LAm  +  
                   frailty(imfcode), 
                   data = CompleteData, 
                   control = coxph.control(
                     eps = 1e-09, iter.max = 100, outer.max = 100))

summary(M1)

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

# Keep qmv and Lqmv
Drawn <- data.frame(Drawn[, c(1, 3, 8)])

# Create Merge Variable (Simulation Number)
Drawn$ID <- 1:1000

# Create range of years over which the combined coefficient will be
# created
Range <- seq(from = 1821, to = 2008, by = 2)

# Multiply the time-varying coeff by time
TVSim <- outer(Drawn[, 1], Range)
TVSim <- data.frame(melt(TVSim))

TVSim <- rename(TVSim, c(X1 = "ID", X2 = "Time", value = "TVR"))

#### Create Combined Coefficient
# Combine base variables
Drawn$Base <- Drawn[,2] + Drawn[,3]
BaseData <- Drawn[, c("ID", "Base")]

# Add to time interaction to the base variables
TVSim <- merge(BaseData, TVSim, by = "ID")
TVSim$CC <- TVSim$Base + TVSim$TVR

# Create Combined Hazard Ratio
TVSim$HRCC <- exp(TVSim$CC)

# Order Variables
TVSim <- TVSim[order(TVSim$Time), ]
TVSim <- TVSim[order(TVSim$ID), ]

# Remove values outside of the 2.5% and 97.5% quantiles
# Find 2.5% and 97.5% quantiles for HRCC
TVSimPerc <- ddply(TVSim, .(Time), transform, Lower = HRCC < quantile(HRCC, c(0.025)))

TVSimPerc <- ddply(TVSimPerc, .(Time), transform, Upper = HRCC > quantile(HRCC, c(0.975)))

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

#### Graph Simulated Combined Hazard Ratios ####
# Plot
ggplot(TVSimPerc, aes(x = Time, y = HRCC)) +
  geom_point(alpha = I(0.06) ,colour = "gray60") + 
  geom_smooth() + 
  geom_hline(aes(yintercept = 1), linetype = "dotted") + 
  scale_y_continuous(limits = c(0, 4)) + 
  scale_x_continuous(breaks = c(0, 94), labels = c(1821, 2008)) + 
  xlab("") + ylab("Simulated Relative Hazard\n") + 
  theme_bw(base_size = 15)