Let’s load the packages we will need and the UN dataset.

#---- Packages 
library(foreign)
library(survival)
library(survminer)
library(flexsurv)
library(texreg)
library(ggplot2)
library(SurvRegCensCov)
#-------Dataset
un <- read.dta("~/Dropbox/event_history/dta/UNFINAL.dta")

We create the survival object for our dependent variable, using the syntax Surv(time, event).

#--------Data prepation

dv <- Surv(un$'_t', un$'_d')

Table 3.1

First we run the exponential and Weibull A.F.T models

# Exponential Model

exp_mod <- survreg(dv ~ civil + interst, data = un, dist = "exp")

# Weibull A.F.T

weib_mod <- survreg(dv ~ civil + interst, data = un, dist = "weib")

The Weibull Prop. Hazards model is a little more involved to run, because the default in R is the accelerated failure time (A.F.T.) parameterization. To convert the Weibull A.F.T. to Prop. Hazards, we will use the ConvertWeibull function. Because this function does not calculate the intercept or standard errors, however, we will first have to create an intercept term manually by entering a vector of 1s in the dataset. We will then enter the intercept term as a covariate in the model.

intercept <- rep(1,58)

un$intercept <- intercept

We run the Weibull A.F.T models with this intercept term. The ConvertWeibull function will drop the first term so we have to two identical models with the terms ordered differently to get the conversions for all the variables.

weib_mod1 <- survreg(dv ~ 0 + civil + intercept + interst, data = un)

weib_mod2 <- survreg(dv ~ 0 + intercept + civil + interst, data = un)

weib_mod1_ph <- ConvertWeibull(weib_mod1)

weib_mod2_ph <- ConvertWeibull(weib_mod2)

We can extract the coefficients and standard errors from these two models, and delete the redundant rows.

#Extract the coefficients and standard errors

weib_ph <- rbind(weib_mod1_ph$vars, weib_mod2_ph$vars)

# Delete rows that are redundant/unnecessary 
weib_ph <- weib_ph[-c(1,2,5,8),]

weib_ph <- weib_ph[c("intercept", "civil", "interst", "gamma"),]

weib_ph
##             Estimate         SE
## intercept -3.4599094 0.49528585
## civil      0.8879243 0.38320176
## interst   -1.4014415 0.51178181
## gamma      0.8068950 0.09988463
Weibull Model of U.N. Peacekeeping Missions
Exponential Model Weibull A.F.T. Weibull Prop. Hazards
Constant 4.35 (0.21) 4.29 (0.27) -3.46 (0.50)
Civil War -1.17 (0.36) -1.10 (0.45) 0.89 (0.38)
Interstate Conflict 1.64 (0.50) 1.74 (0.62) -1.40 (0.51)
Shape Parameter 1.00 1.24 0.81
Log Likelihood -202.85 -201.15 -201.15
N 54 54 54

Figure 3.1

# Using the Weibull A.F.T. model we ran earlier, we first calculate lambda for each conflict type, using the fact that lambda = exp(-beta'x)
lambda_civil = unname(exp(-(weib_mod$coef[1] + weib_mod$coef[2])))

lambda_interst = unname(exp(-(weib_mod$coef[1] + weib_mod$coef[3])))

lambda_icw <- unname(exp(-(weib_mod$coef[1])))
# We can generate the hazard rate for each covariate profile, knowing that h(t) = 
# lambda * p * (lambda * t)^(p-1)

p <- 1/weib_mod$scale
t <- un$'_t'

hazard_civil <- lambda_civil * p * (lambda_civil * t)^(p-1)

hazard_interstate <- lambda_interst * p * (lambda_interst * t)^(p-1)

hazard_icw <- lambda_icw * p * (lambda_icw * t)^(p-1)
#Plot hazard rates

ggplot(data = un, aes(x = t)) + 
  #geom_point(aes(y = hazard_civil)) +
  geom_line(aes(y = hazard_civil, colour = "Civil War")) +
  #geom_point(aes(y = hazard_interstate)) +
  geom_line(aes(y = hazard_interstate, colour = "Interstate Conflict")) +
  #geom_point(aes(y = hazard_icw)) +
  geom_line(aes(y = hazard_icw, colour = "Internationalized Civil War")) +
  theme_bw() +
  ylab(label = "Hazard Rates") +  xlab("Duration of U.N. Peacekeeping Missions") +
  labs(colour = "Type of Conflict") +
  ggtitle("Weibull Hazard Rates by Type of Conflict")