Note that for two of the models, the replication of Table 5.3 will be based on the smaller dataset (which I name “imi”) than the one shown in Stata and in the book. This is because the exact discrete Cox model using the larger dataset (“imi2”) is too computionally intense for R to estimate within a reasonable time. In addition, the Weibull model cannot handle start-stop Surv objects. Nevertheless, the results using the smaller dataset are essentially the same so I still reproduce them below.

First, we prepare both the large and small datasets for analysis by renaming certain variables to avoid conflicting with R’s base functions, creating the Surv object, and creating mean centered variables for analysis of the larger dataset.

#-------Dataset

# This is the dataset used in the Stata version, but because it takes so long to run, one can get the same results with the smaller dataset. This is important when running the exact discrete Cox model, because that model is computationally intense. 

imi2 <- read.dta("~/Dropbox/event_history/dta/omifinal2.dta")


imi <- read.dta("~/Dropbox/event_history/dta/omismall.dta")


#-------Data prepation

# Rename dependent variable for both datasets to avoid confusion with base R functions

names(imi2)[names(imi2) == 'break'] <- 'breakdown'

# Create survival object for dv in the larger dataset

dv2 <- Surv(imi2$'_t0', imi2$'_t', imi2$event)

# Create survival object for dv in the smller dataset

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

# We want to only run the models on the data that has no missing observations for the following logit equation. 

#m1 <- glm(fail ~ ctg + ali + idem + tdem + pbal + breakdown, data = imi2)

imi2 <- imi2 %>% 
  drop_na(fail, ctg, ali, idem, tdem, pbal, breakdown)


# Lastly, we mean center the pbal, idem, and tdem variables.

imi2$idemmean <- imi2$idem - mean(imi2$idem)

imi2$tdemmean <- imi2$tdem - mean(imi2$tdem)

imi2$pbalmean <- imi2$pbal - mean(imi2$pbal)

Table 5.3

We first run the Cox exact discrete model (which is essentially the same as the conditoinal logit except it does not drop observations) and the Weibull model using the smaller dataset.

# Cox exact discrete

ed_mod <- coxph(dv ~ pbal + ctg + ali + idem + tdem + breakdown, data = imi, ties = "exact")

# Weibull A.F.T. 

weib_mod <- survreg(dv ~ pbal + ctg + ali + idem + tdem + breakdown, data = imi, dist = "weibull")

Again, we have to convert the Weibull A.F.T. to a P.H. parameterization.

# Conversion from Weibull A.F.T.  to P.H. 

intercept <- rep(1, 656)

imi$intercept <- intercept

weib_mod1 <- survreg(dv ~ 0 + intercept +  pbal + ctg + ali + idem + tdem + breakdown, data = imi, dist = "weibull")

weib_mod2 <- survreg(dv ~ 0 + pbal + intercept +  ctg + ali + idem + tdem + breakdown, data = imi, dist = "weibull")

weib_mod1_ph <- ConvertWeibull(weib_mod1)

weib_mod2_ph <- ConvertWeibull(weib_mod2)

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

# Delete rows that are redundant/unnecessary 
weib_ph <- weib_ph[-c(1,2,9,12,13,14,15,16),]

weib_ph <- weib_ph[c("intercept", "pbal", "ctg", "ali", "idem", "tdem",
                     "breakdown", "gamma"),]

weib_ph
##              Estimate          SE
## intercept -1.32153217 0.151345243
## pbal      -0.50188655 0.152589537
## ctg       -0.28344841 0.101406267
## ali        0.25609560 0.095594264
## idem       0.01247340 0.006120223
## tdem       0.02285714 0.006710263
## breakdown -0.44311104 0.197807373
## gamma      0.65543389 0.022078980

Let’s format these models into a presentable table. I show how to do this for the exact discrete model and repeat the same steps for the Weibull (not shown).

# Create texreg object for exact discrete/conditional logit model. 

tex_reg_converter <- function(original_model){
  
  
  coefficient.names <- c("Relative Capabilities", "Territorial Contiguity", "Intervenor Allied to Target", "Intervenor Democracy", "Target Democracy",
                       "Breakdown of Authority")
  coefficients <- original_model$coef 
  se <- sqrt(diag(vcov(original_model)))
  n <- original_model$n
  lik <- logLik(original_model)
  gof <- c(lik, n)
  gof.names <- c("Log Likelihood", "N" )
  decimal.places <- c(TRUE, FALSE)
  
  createTexreg(
    coef.names = coefficient.names, coef = coefficients, se = se,
    gof.names = gof.names, gof = gof, gof.decimal = decimal.places)
  
}

ed_mod <- tex_reg_converter(ed_mod)

Now we can run the logit model with the lowess term using the large dataset.

log_mod <- glm(event ~ pbalmean + ctg + ali + idemmean + tdemmean + breakdown + lowesst2, data = imi2, family = "binomial")

We format the results of the logit and put them together with the Cox and Weibull results.

tex_reg_converter <- function(original_model){

coefficient.names <- c("Constant","Relative Capabilities", "Territorial Contiguity", "Intervenor Allied to Target", "Intervenor Democracy", "Target Democracy", "Breakdown of Authority", "Duration Dependency")
coefficients <- original_model$coef 
se <- sqrt(diag(vcov(original_model)))
n <- 9374
lik <- logLik(original_model)
gof <- c(lik, n)
gof.names <- c("Log Likelihood", "N" )
decimal.places <- c(TRUE, FALSE)

createTexreg(
    coef.names = coefficient.names, coef = coefficients, se = se,
    gof.names = gof.names, gof = gof, gof.decimal = decimal.places)
  
}

logit_mod <- tex_reg_converter(log_mod)
# Put all the tables together

htmlreg(list(ed_mod, logit_mod, weib_mod_ph), stars = c(), caption = "Models of Militarized Interventions", caption.above = TRUE, custom.model.names = c("Conditional Logit", "Logit", "Weibull"), center = FALSE, single.row = TRUE)
Models of Militarized Interventions
Conditional Logit Logit Weibull
Relative Capabilities -0.48 (0.16) -0.43 (0.16) -0.50 (0.15)
Territorial Contiguity -0.26 (0.11) -0.25 (0.11) -0.28 (0.10)
Intervenor Allied to Target 0.24 (0.10) 0.22 (0.10) 0.26 (0.10)
Intervenor Democracy 0.01 (0.01) 0.01 (0.01) 0.01 (0.01)
Target Democracy 0.02 (0.01) 0.02 (0.01) 0.02 (0.01)
Breakdown of Authority -0.46 (0.21) -0.43 (0.20) -0.44 (0.20)
Constant -4.07 (0.13) -1.32 (0.15)
Duration Dependency 16.20 (0.95)
Log Likelihood -1591.49 -1779.31 -1931.50
N 520 9374 520
Shape Parameter 0.66