success_of_coups.Rmd

---
title: "Machine Learning Models to Predict Coup Success"
date: December 16, 2019
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


## Getting started

*Introduction*: The goal of this script is to find the strongest predictors of coup success, coded in the dataset as `realized_coup`, and then to build a the best possible predictive model of coup success. The dataset used is publicly available from the Cline Center at the University of Illinois, and it includes over 50 coup features for more than 1,000 coups since the 1940s.

```{r load}

library(tidyverse)
require(stargazer)
library(stats)
library(dplyr)
library(janitor)
library(boot)
library(knitr)
library(e1071)
library(rpart)
library(rpart.plot)
library(ggplot2)
library(randomForest)
library(gbm)
library(Hmisc)

# Create a dataframe: coup_df

coup_csv <- read_csv("coupdata.csv")
colnames(coup_csv) %<>% stringr::str_replace_all("\\s","_") %>% tolower

coup_df <- data.frame(coup_csv)
colnames(coup_df)
coup_df <- coup_df %>%
  filter(!is.na(realized_coup))

coup_df <- coup_df %>% clean_names()

# Manually edit colnames of dataframe

colnames(coup_df)[which(names(coup_df) == "ambigious_coup")] <- "ambiguous_coup"
colnames(coup_df)[which(names(coup_df) == "were_others_then_the_incumbent_injured")] <- "were_others_than_the_incumbent_injured"
```

## Construct potentially relevant variables

New binary category variables include decade (seventies, sixties, etc.) and season (spring, summer and so on).

```{r variables, warning=FALSE}

# Eliminate nulls from the dataset

sum(is.na(coup_df))

coup_df <- coup_df %>%
  mutate(were_ordinary_citizens_involved = (is.na(were_ordinary_citizens_involved) <- 0)) %>%
  mutate(were_ordinary_citizens_involved = as.numeric(were_ordinary_citizens_involved))

sum(is.na(coup_df$were_ordinary_citizens_involved))

# Comment: We assume that except where ordinary citizens were known to be involved, all others event should be coded as having negligible or no ordinary citizen involvement.

# Code coup type as numeric instead of character ("conspiracy", "attempt," etc.)

coup_df <- coup_df %>%
  mutate(coup_type_numeric = type_of_coup)

coup_df$coup_type_numeric <- gsub("conspiracy", "1", coup_df$coup_type_numeric)
coup_df$coup_type_numeric <- gsub("attempt", "2", coup_df$coup_type_numeric)
coup_df$coup_type_numeric <- gsub("coup", "3", coup_df$coup_type_numeric)
coup_df$coup_type_numeric <- as.numeric(coup_df$coup_type_numeric)

# Create season variables

head(coup_df$month_of_event)
typeof(coup_df$month_of_event)

coup_df <- coup_df %>%
  mutate(winter = case_when(
    month_of_event == (12 | 1:2) ~ 1,
    TRUE ~ 0))

coup_df <- coup_df %>%
  mutate(spring = case_when(
    month_of_event == (3:5) ~ 1,
    TRUE ~ 0))

coup_df <- coup_df %>%
  mutate(summer = case_when(
    month_of_event == (6:8) ~ 1,
    TRUE ~ 0))

coup_df <- coup_df %>%
  mutate(autumn = case_when(
    month_of_event == (9:11) ~ 1,
    TRUE ~ 0))

sum(coup_df$summer) == sum(coup_df$spring)

# Create decade variables

typeof(coup_df$year)
head(coup_df$year)

coup_df <- coup_df %>%
  mutate(forties = case_when(
    (year > 1939) & (year < 1950) ~ 1,
    TRUE  ~ 0)) %>%
  mutate(fifties = case_when(
    (year > 1949) & (year < 1960) ~ 1,
    TRUE ~ 0)) %>%
  mutate(sixties = case_when(
    (year > 1959) & (year < 1970) ~ 1,
    TRUE  ~ 0)) %>%
  mutate(seventies = case_when(
    (year > 1969) & (year < 1980) ~ 1,
    TRUE  ~ 0)) %>%
  mutate(eighties = case_when(
    (year > 1979) & (year < 1990) ~ 1,
    TRUE  ~ 0)) %>%
  mutate(nineties = case_when(
    (year > 1989) & (year < 2000) ~ 1,
    TRUE  ~ 0)) %>%
  mutate(aughties = case_when(
    (year > 1999) & (year < 2010) ~ 1,
    TRUE  ~ 0))

sum(coup_df$forties)
sum(coup_df$fifties)
sum(coup_df$sixties)
sum(coup_df$seventies)
sum(coup_df$eighties)
sum(coup_df$nineties)
sum(coup_df$aughties)

```


## Export the data as csv

For use in Tableau or Python, export the post-data wrangling clean dataset.

```{r export, warning=FALSE}

write.csv(coup_df,"coup_data_tableau_python.csv", row.names = FALSE)

```

## Drop irrelevant variables

Use select(-x) to drop variables that are irrrelevant to coup success.

```{r drop, warning=FALSE}

coup_df <- coup_df %>% select(-x54, -x53, -country, -coup_id, -cow_code, -day_of_event)
coup_df <- coup_df %>% select(-year, -month_of_event, -type_of_coup, -unrealized, 
                              -coup_conspiracies, -attempted_coup, -coup_type_numeric)


```

## Split the data into train and test

```{r split, warning=FALSE}

set.seed(1)
train <- sample_frac(coup_df, 0.7)
test <- dplyr::setdiff(coup_df, train)

sum(is.na(test))


```

## Create a simple logistic regression model for coup violence with a promising predictor: `was_a_foreign_gov_involved`

```{r foreign, warning=FALSE}

mean(coup_df$realized_coup)
filter(coup_df, was_a_foreign_gov_involved == 1) %>% summarise(mean(realized_coup), na.rm = TRUE)

# Over 60% of foreign-backed coups succeed, while only about 40% of attempted coups overall succeed

logit_form_1 <- formula(realized_coup ~ was_a_foreign_gov_involved)

foreign_gov_fit <- glm(logit_form_1, data = train, family = binomial(link = 'logit'))

summary(foreign_gov_fit)

# Comment: A coup's status as being foreign-orchestrated is a strong predictor of success, based on the logistic regression

train <- train %>% 
  mutate(foreign_gov_logit = foreign_gov_fit$fitted.values, 
         foreign_gov_success_pred = as.numeric(foreign_gov_logit > 0.5))

error_rate_train <- train %>% 
  mutate(log_error = realized_coup != foreign_gov_success_pred) %>% 
  summarise(log_error_rate = sum(log_error)/n())
error_rate_train

train %>% 
  group_by(realized_coup, foreign_gov_success_pred) %>% 
  summarise(n = n()) %>% 
  spread(realized_coup, n)

# Comment: The model had a 37.5% failure rate in the training data. What about the test data?

test <- test %>%
  mutate(foreign_gov_pred = predict(foreign_gov_fit,., type = "response"), 
         foreign_gov_success_pred = as.numeric(foreign_gov_pred > 0.5))

error_rate_test <- test %>% 
  mutate(log_error = realized_coup != foreign_gov_success_pred) %>% 
  summarise(log_error_rate = sum(log_error)/n())
error_rate_test

# Comment: The model had a 44.1% failure rate in the test data. Thus, despite a real correlation between foreign involvement
# and coup success, this cannot be the sole basis of a predictive model.

test %>% 
  group_by(realized_coup, foreign_gov_success_pred) %>% 
  summarise(n = n()) %>% 
  spread(realized_coup, n)

# Remove predictions from test, train

train <- train %>% select(-foreign_gov_success_pred, -foreign_gov_logit)
test <- test %>% select(-foreign_gov_success_pred, -foreign_gov_pred)


```

## Create a full regression model of coup success for all potential predictors

```{r full_model, warning=FALSE}


logit_form_all <- formula(realized_coup ~ forties + fifties + sixties + seventies + eighties + nineties + aughties +
                            spring + summer + autumn + winter + dissident_actions + foreign_coup + counter_coup + auto_coup + 
                            ambiguous_coup + were_military_actors_involved_in_the_coup +
                            were_other_security_officials_involved_in_the_coup + were_rebel_soldiers_involved_in_the_coup + 
                            were_non_military_government_officials_involved_in_the_coup + 
                            were_non_government_political_actors_involved_in_the_coup + were_political_radicals_involved_in_the_coup +
                            were_business_leaders_involved_in_the_coup + were_organized_labor_activists_involved_in_the_coup +
                            were_students_or_academics_involved_in_the_coup + were_religious_leaders_involved_in_the_coup + 
                            were_ethnic_group_leaders_involved + were_political_radicals_involved_in_the_coup + 
                            were_business_leaders_involved_in_the_coup + were_ordinary_citizens_involved + were_mercenaries_involved +
                            were_non_government_foreign_actors_involved + was_a_foreign_gov_involved + were_small_arms_used + 
                            were_explosive_devices_used + were_military_grade_weapons_used)


lm_coup_all <- glm(logit_form_all, data = train, family = binomial(link = 'logit'))

summary(lm_coup_all)

stargazer(lm_coup_all, type = "text")


# Comments: The full model does contain a number of significant predictors. At the highest level of significance (p < 0.001),
# predictors include dissident involvement (negative effect), civilian officials (positive), and  
# ordinary citizen involvement (positive). At a lower level of significance, foreign involvement  (positive), auto-coup (positive), 
# rebel soldiers (positive), and military officers (positive) all contribute to a coup's success.

train <- train %>% 
  mutate(coup_logit_all = lm_coup_all$fitted.values, 
         realized_coup_pred = as.numeric(coup_logit_all > 0.5))

error_rate_train <- train %>% 
  mutate(log_error = realized_coup != realized_coup_pred) %>% 
  summarise(log_error_rate = sum(log_error)/n())
error_rate_train

train %>% 
  group_by(realized_coup, realized_coup_pred) %>% 
  summarise(n = n()) %>% 
  spread(realized_coup, n)

# Comment: The model has a 26.1% failure rate in the training data. What about the test data?

test <- test %>%
  mutate(coup_logit_all = predict(lm_coup_all,., type = "response"), 
         realized_coup_pred = as.numeric(coup_logit_all > 0.5))

error_rate_test <- test %>% 
  mutate(log_error = realized_coup != realized_coup_pred) %>% 
  summarise(log_error_rate = sum(log_error)/n())
error_rate_test

# Comment: The model had a 33.3% failure rate in the test data. This is hardly an improvement on foreign involvement alone.

test %>% 
  group_by(realized_coup, realized_coup_pred) %>% 
  summarise(n = n()) %>% 
  spread(realized_coup, n)

# Remove predictions from test data

train <- train %>% select(-realized_coup_pred, -coup_logit_all)
test <- test %>% select(-realized_coup_pred, -coup_logit_all)


```

## Forward Stepwise Model Selection Method: Choosing the best combination of predictors

The following algorithm selects the logistic regression model with the lowest error rate by testing all possible combinations of predictors.

```{r stepwise, warning = FALSE}

# Function to calculate error rate for any logistic regression model

cv_fun <- function(f) {
  glmfit <- glm(f, data = train, family = binomial(link = 'logit'))
  stepwise_pred <- predict(glmfit, newdata = test, type = "response")
  stepwise_binary <- as.numeric(stepwise_pred > 0.5)
  log_error <- case_when((test$realized_coup == stepwise_binary) ~ 0,
                         TRUE ~ 1)
  mean(log_error)
}

# Initial settings
X <- list("military_coup", "rebel_coup", "palace_coup", 
          "popular_revolt", "dissident_actions", "foreign_coup", "internationally_mediated_transitions", "forced_resignation",
          "counter_coup", "auto_coup","ambiguous_coup", 
          "were_military_actors_involved_in_the_coup",
          "were_other_security_officials_involved_in_the_coup", "were_rebel_soldiers_involved_in_the_coup",
          "were_non_military_government_officials_involved_in_the_coup", "were_others_than_the_incumbent_killed",
          "were_non_government_political_actors_involved_in_the_coup", "were_political_radicals_involved_in_the_coup",
          "were_business_leaders_involved_in_the_coup", "were_organized_labor_activists_involved_in_the_coup",
          "were_students_or_academics_involved_in_the_coup", "were_religious_leaders_involved_in_the_coup",
          "were_ethnic_group_leaders_involved", "were_ordinary_citizens_involved", "were_mercenaries_involved",
          "were_non_government_foreign_actors_involved", "was_a_foreign_gov_involved", "winter",
          "spring", "summer", "autumn", "forties", "fifties", "sixties", "seventies", "eighties", "nineties", "aughties")
length(X)
f_best = vector("list", length(X) + 1)
f_best_cv = vector("list", length(X) + 1)
i <- 1

while(length(X) > 1) {
  if (i > 1) {
    # Best model from previous iteration
    f_old <- f_best[[i-1]]
    
    # Unused controls
    X <- as.character(f_old)[3] %>% 
      strsplit(split = " + ", fixed = TRUE) %>% 
      unlist() %>% 
      setdiff(X, .)
    
    # All models adding one unused control to f_old
    f_new <- paste(as.character(f_old)[2], "~", as.character(f_old)[3], "+", X) %>% 
      lapply(as.formula, env = .GlobalEnv)
    
    # Cross-Validation for all new models
    f_new_cv <- lapply(f_new, cv_fun)
    
    # Save best model and Error Rate
    f_best[[i]] <- f_new[[which.min(f_new_cv)]]
    f_best_cv[[i]] <- f_new_cv[[which.min(f_new_cv)]]
  } else {
    # Constant-only model and Error Rate
    f_best[[i]] <- formula(realized_coup ~ 1)
    f_best_cv[[i]] <- cv_fun(formula(realized_coup ~ 1))
  }
  i <- i + 1
  print(length(X))
}

# Which model has the lowest error rate?

f_best_forward <- f_best[[which.min(f_best_cv)]]
f_best_forward
f_best_forward_error <- f_best_cv[[which.min(f_best_cv)]]
f_best_forward_error


```


## Examine the best model from the selection algorithm

It is important to see if the model can be tweaked to lower false negatives, since coups are rare and significant enough that we should err on the side of predicting too many, not too few. This tweaking can be done by changing the threshold at which a value counts as a positive prediction from 0.5 to something lower, like 0.25.

To do so, we must look closely at the model's confusion matrices.

```{r selected_model, warning = FALSE}

logit_form_stepwise <- formula(realized_coup ~ 1 + palace_coup + popular_revolt + was_a_foreign_gov_involved + dissident_actions +
                                 auto_coup + internationally_mediated_transitions + military_coup + 
                                 were_non_government_foreign_actors_involved)

lm_coup_stepwise <- glm(logit_form_stepwise, data = train, family = binomial(link = 'logit'))

summary(lm_coup_stepwise)

train <- train %>% 
  mutate(step_coup_pred = lm_coup_stepwise$fitted.values, 
         step_coup_pred = as.numeric(step_coup_pred > 0.5))

error_rate_train <- train %>% 
  mutate(log_error = realized_coup != step_coup_pred) %>% 
  summarise(log_error_rate = sum(log_error)/n())
error_rate_train

train %>% 
  group_by(realized_coup, step_coup_pred) %>% 
  summarise(n = n()) %>% 
  spread(realized_coup, n)

sum(train$realized_coup)

# Comment: Despite an overall error rate that sounds moderately high (26.2%), the model produces useful predictions, predicting coup success in over half of coups (156/305).

test <- test %>%
  mutate(step_coup_pred = predict(lm_coup_stepwise,., type = "response"), 
         step_coup_pred = as.numeric(step_coup_pred > 0.5))

error_rate_test <- test %>% 
  mutate(log_error = realized_coup != step_coup_pred) %>% 
  summarise(log_error_rate = sum(log_error)/n())
error_rate_test

test %>% 
  group_by(realized_coup, step_coup_pred) %>% 
  summarise(n = n()) %>% 
  spread(realized_coup, n)

# Interestingly, while the overall error rate is very slightly higher in the test data, the rate of false negatives is lower, and the model predicts success in ~2/3 of successful coups.

# Remove predictions from test data

train <- train %>% select(-step_coup_pred)
test <- test %>% select(-step_coup_pred)


```

## Naive Bayes Classifier

```{r naive_bayes, warning=FALSE}
library(naivebayes)

X <- train %>% select(forties, fifties, sixties, seventies, eighties, nineties, aughties,
                            spring, summer, autumn, winter, dissident_actions, foreign_coup, counter_coup, auto_coup, 
                            ambiguous_coup, were_military_actors_involved_in_the_coup,
                            were_other_security_officials_involved_in_the_coup, were_rebel_soldiers_involved_in_the_coup, 
                            were_non_military_government_officials_involved_in_the_coup, 
                            were_non_government_political_actors_involved_in_the_coup, were_political_radicals_involved_in_the_coup,
                            were_business_leaders_involved_in_the_coup, were_organized_labor_activists_involved_in_the_coup,
                            were_students_or_academics_involved_in_the_coup, were_religious_leaders_involved_in_the_coup, 
                            were_ethnic_group_leaders_involved, were_political_radicals_involved_in_the_coup, 
                            were_business_leaders_involved_in_the_coup, were_ordinary_citizens_involved, were_mercenaries_involved,
                            were_non_government_foreign_actors_involved, was_a_foreign_gov_involved, were_small_arms_used, 
                            were_explosive_devices_used, were_military_grade_weapons_used)

X <- data.matrix(X, rownames.force = NA)

Y <- train %>% select(realized_coup) %>% mutate(realized_coup = eval(realized_coup > 0.5))
Y <- as.vector(t(Y))

Bern_Model <- bernoulli_naive_bayes(X, Y)

# Model summary

Bern_Model

# Test the model

A <- test %>% select(forties, fifties, sixties, seventies, eighties, nineties, aughties,
                            spring, summer, autumn, winter, dissident_actions, foreign_coup, counter_coup, auto_coup, 
                            ambiguous_coup, were_military_actors_involved_in_the_coup,
                            were_other_security_officials_involved_in_the_coup, were_rebel_soldiers_involved_in_the_coup, 
                            were_non_military_government_officials_involved_in_the_coup, 
                            were_non_government_political_actors_involved_in_the_coup, were_political_radicals_involved_in_the_coup,
                            were_business_leaders_involved_in_the_coup, were_organized_labor_activists_involved_in_the_coup,
                            were_students_or_academics_involved_in_the_coup, were_religious_leaders_involved_in_the_coup, 
                            were_ethnic_group_leaders_involved, were_political_radicals_involved_in_the_coup, 
                            were_business_leaders_involved_in_the_coup, were_ordinary_citizens_involved, were_mercenaries_involved,
                            were_non_government_foreign_actors_involved, was_a_foreign_gov_involved, were_small_arms_used, 
                            were_explosive_devices_used, were_military_grade_weapons_used)

A <- data.matrix(A, rownames.force = NA)

B <- test %>% select(realized_coup) %>% mutate(realized_coup = eval(realized_coup > 0.5))
B <- as.vector(t(B))

nb_pred <- predict(Bern_Model, A)
table(nb_pred,B)

nb_error_rate <- sum(nb_pred != B)/(sum(nb_pred == B) + sum(nb_pred != B))
nb_error_rate

# A model trained using the Python implementation of Bernoulli Naive Bayes from sklearn yields better predictions

```


## Tree-based Classification: Decision Tree Models

```{r tree, warning=FALSE}

# Use the full formula with all predictors employed in the full logistic regression: `logit_form_all`

# Train the decision tree

tree_realized <- rpart(logit_form_all, data = train)

tree_realized

# Visualize the tree

plot(tree_realized, uniform = TRUE)
text(tree_realized, pretty=0)

#  Prune the tree

index <- which.min(tree_realized$cptable[ , "xerror"])
tree_min <- tree_realized$cptable[index, "CP"]
prune_tree <- prune(tree_realized, cp = tree_min)

prune_tree

# Comment: No pruning is needed

# Visualize the pruned tree

plot(prune_tree, uniform = TRUE)
text(prune_tree, pretty=0)

prp(prune_tree, extra = 1, box.palette = "auto")

# Predict using the test data and transform into predictions

tree_predict <- predict(tree_realized, newdata = test)
tree_predict <- as.numeric(tree_predict > 0.5)

# Means squared error (MSE)

mean((tree_predict - test$realized_coup)^2)

# Comment: MSE is different from the error rate in previous models, so let's use a different metric: error rate

# Compare predicted success coups with actual successful coups, i.e., the error rate on successful coups

tree_predict <- predict(tree_realized, newdata = test)
tree_predict <- as.numeric(tree_predict > 0.5)

error_rate_test <- test %>% 
  mutate(log_error = realized_coup != tree_predict) %>% 
  summarise(log_error_rate = sum(log_error)/n())
error_rate_test

test %>% 
  mutate(tree_predict = tree_predict) %>%
  group_by(realized_coup, tree_predict) %>% 
  summarise(n = n()) %>% 
  spread(realized_coup, n)

# Comment: Of 106 successful coups, this model predicted around half. The overall error rate is 37%.


```

## Using Random Forest to predict coup success

Random Forest is a model building method that builds a decision tree by randomly leaving out whole variables to prevent overfiting.

```{r random, warning=FALSE}

rf_fit <- randomForest(logit_form_all, data = train, mtry = 6, importance = TRUE)

summary(rf_fit)

# MSE

rf_test <- predict(rf_fit, newdata = test)

rf_mse <- mean((rf_test - test$realized_coup)^2)

rf_mse

# Visualize error rates

plot(rf_fit, uniform = TRUE, main = "Errors By Tree Size, Random Forest Models, `coup_success`")

# Prediction Error Rate Test: 0.5 prediction threshold

rf_test_binary_1 <- as.numeric(rf_test > 0.5)

error_rate_test <- test %>% 
  mutate(log_error = realized_coup != rf_test_binary_1) %>% 
  summarise(log_error_rate = sum(log_error)/n())
error_rate_test

test %>% 
  mutate(rf_test_binary_1 = rf_test_binary_1) %>%
  group_by(realized_coup, rf_test_binary_1) %>% 
  summarise(n = n()) %>% 
  spread(realized_coup, n)

# Comment: These results are good. The overall error rate is around 32%, and 2/3 of successful coups are predicted. But what happens if the prediction threshold is lowered to 0.25?

rf_test_binary_2 <- as.numeric(rf_test > 0.25)

error_rate_test <- test %>% 
  mutate(log_error = realized_coup != rf_test_binary_2) %>% 
  summarise(log_error_rate = sum(log_error)/n())
error_rate_test

test %>% 
  mutate(rf_test_binary_2 = rf_test_binary_2) %>%
  group_by(realized_coup, rf_test_binary_2) %>% 
  summarise(n = n()) %>% 
  spread(realized_coup, n)

# Comment: At this lower prediction level, the error rate rises slightly from 37% to 39.6%, but now 86% of violent coups are predicted, with a tolerable rise in false positives.

# Determine node purity and MSE by predictor

# varImpPlot(rf_fit)
impToPlot <- importance(rf_fit, scale=FALSE)

# write.csv(impToPlot,"varImpPlot_success.csv", row.names = TRUE)


```


## Bonus: Find model predictions for two recent coups, Egypt (2013) and Turkey (2016)

I have coded the data related to these two coups myself, based on contemporary reports (outside the dataset). Note that the Egyptian coup succeeded, whereas Turkey's failed.

```{r bonus, warning=FALSE}

# Load the datasets and deselect 

egypt_csv <- read_csv("egypt_coup.csv")
egypt_df <- data.frame(egypt_csv)

turkey_csv <- read_csv("turkey_coup.csv")
turkey_df <- data.frame(turkey_csv)

length(turkey_df) == length(test)

# Predict 2013 Egyptian coup outcome with algorithm-selected logistic regression

egypt_pred <- predict(lm_coup_stepwise, newdata = egypt_df, type = "response")

egypt_pred

# Predict 2013 Egyptian coup outcome with the Naive Bayes model

E <- data.matrix(egypt_df, rownames.force = NA)

nb_egypt <- predict(Bern_Model, E)

nb_egypt

# Predict 2013 Egyptian coup outcome with the Decision Tree model

tree_egypt <- predict(tree_realized, newdata = egypt_df)

tree_egypt

# Predict 2013 Egyptian coup outcome with the Random Forest model

rf_egypt <- predict(rf_fit, newdata = egypt_df)

rf_egypt

# Predict 2016 Turkish coup outcome with algorithm-selected logistic regression

turkey_pred <- predict(lm_coup_stepwise, newdata = turkey_df, type = "response")

turkey_pred

# Predict 2016 Turkish coup outcome with the Naive Bayes model

T <- data.matrix(turkey_df, rownames.force = NA)

nb_turkey <- predict(Bern_Model, T)

nb_turkey

# Predict 2016 Turkish coup outcome with the Decision Tree model

tree_turkey <- predict(tree_realized, newdata = turkey_df)

tree_turkey

# Predict 2016 Turkish coup outcome with Random Forest

rf_turkey <- predict(rf_fit, newdata = turkey_df)

rf_turkey

```

Clearly, these models have predictive power. Using a prediction threshold of 0.5, all three models correctly predicted the outcome of  the Egyptian and Turkish coups. In both cases, the stepwise selected logistic regression was most predictive.

# *Conclusion*: Machine learning models are effective in predicting coup success, though not perfectly. A logistic regression built through a forward stepwise selection algorithm produced the best results, but tree-based classication methods, whether a Decision Tree or Random Forest model, performs nearly as well.