11_ppcs_regressions.Rmd

---
title: "PPCS Regressions"
author: "Brenda Fried"
date: "8/1/2019"
output: html_document
---

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

```{r}
if (!require("pacman")) install.packages("pacman")
pacman::p_load(tidyverse, pROC, ROCR, scales, here)

load(here("clean_data", "summarized_ppcs.RData"))
```

```{r}
# no control
logit <- glm(force ~ civilian_race, family = "binomial", data= merged_ppcs)
exp(coef(logit))
summary(logit)

# add civilian demographics
logit2 <- glm(force ~ civilian_race + civilian_gender + civilian_income + civilian_employed + population_size + I(civilian_age^2), family='binomial', data = merged_ppcs)
summary(logit2)
exp(coef(logit2))

# add encounter characteristics 
logit3 <- glm(force ~ civilian_race + civilian_gender + civilian_income + civilian_employed + population_size + I(civilian_age^2) + time_of_encounter + type_of_incident + off_hispanic + off_white + off_black + off_other + off_split, family='binomial', data = merged_ppcs)
summary(logit3)
exp(coef(logit3))

# add civilian behavior
logit4 <- glm(force ~ civilian_race + civilian_gender + civilian_income + civilian_employed + population_size + I(civilian_age^2) + time_of_encounter + type_of_incident + off_hispanic + off_white + off_black + off_other + off_split + civilian_behavior, family='binomial', data = merged_ppcs)
summary(logit4)
exp(coef(logit4))

# add year
logit5 <- glm(force ~ civilian_race + civilian_gender + civilian_income + civilian_employed + population_size + I(civilian_age^2) + time_of_encounter + type_of_incident + off_hispanic + off_white + off_black + off_other + off_split + civilian_behavior + year, family='binomial', data = merged_ppcs)
summary(logit5)
exp(coef(logit5))

logit_without_race <- glm(force ~ civilian_gender + civilian_income + civilian_employed + population_size + I(civilian_age^2) + time_of_encounter + type_of_incident + off_hispanic + off_white + off_black + off_other + off_split + civilian_behavior + year, family='binomial', data = merged_ppcs)
```

Filter data because of missing variables
Seen from the summary of the models
```{r}
merged_ppcs <- merged_ppcs %>% filter(type_of_incident != "3", !(year %in% c(1996, 1999, 2008)))
```

# Add Predictions
```{r}
# add predictions setting the threshold to 0.03
ppcs_perf_data <- data.frame(actual = merged_ppcs$force, 
                             prob = predict(logit5, merged_ppcs, type = "response")) %>%
  mutate(pred = ifelse(prob > 0.03, 1, 0))
```

# ROC and AUC
```{r}
# create ROC Curve
prob <- as.numeric(round(ppcs_perf_data$prob))
actual <- as.numeric(as.character(ppcs_perf_data$actual))

roc_obj <- roc(prob,actual)
auc(roc_obj)

pred_roc <- prediction(ppcs_perf_data$prob, ppcs_perf_data$actual)
eval <- performance(pred_roc,  measure='tpr', x.measure='fpr')

roc_eval <- data.frame(fpr=unlist(eval@x.values), tpr=unlist(eval@y.values))
ggplot(data=roc_eval, aes(x=fpr, y=tpr)) + 
  theme_bw()+
  geom_line(color = "#451878") +
  geom_abline(linetype=2, color = "#451878") +
  scale_x_continuous(labels=percent, lim=c(0,1)) +
  scale_y_continuous(labels=percent, lim=c(0,1)) +
  xlab('Probability of a False Alarm') +
  ylab('Probability of Detecting Use of Force') +
  labs(title = "ROC Curve of the In-sample Data") +
  theme(axis.title.x = element_text(face = "bold", size=20, color = "#0d693e"),
        axis.title.y = element_text(face = "bold", size=20, color = "#0d693e"),
        axis.text.x = element_text(face = "bold", size=15, hjust = 1),
        axis.text.y = element_text(face = "bold",size=18),
        plot.title = element_text(hjust = 0.5,size=25, color = "#374687",face = "bold"),
        legend.title = element_text(color = "#0d693e", face = "bold", size = 18 ),
        legend.text = element_text(size = 15 ))+
  ggsave("figures/ppcs_roc_curve.png",  width = 12, height = 10, dpi = 150, units = "in", device='png')

table(actual = ppcs_perf_data$actual, predicted = ppcs_perf_data$pred)

```

# Model Performance measures
```{r}
# accuracy: fraction of correct classifications
ppcs_perf_data %>%
  summarize(acc = mean(pred == actual))

# precision: fraction of positive predictions that are actually true
ppcs_perf_data %>%
  filter(pred == 1) %>%
  summarize(prec = mean(actual == 1))

# recall: fraction of true examples that we predicted to be positive
# aka true positive rate, sensitivity
ppcs_perf_data %>%
  filter(actual == 1) %>%
  summarize(recall = mean(pred == 1))

# false positive rate: fraction of false examples that we predicted to be positive
ppcs_perf_data %>%
  filter(actual == 0) %>%
  summarize(fpr = mean(pred == 1))

```


## Running the model without Race as a predictor
```{r}
ppcs_perf_data2 <- data.frame(actual = merged_ppcs$force, 
                             prob = predict(logit5, merged_ppcs, type = "response")) %>%
  mutate(pred = ifelse(prob > 0.03, 1, 0))
```

```{r}
# create ROC Curve
prob <- as.numeric(round(ppcs_perf_data2$prob))
actual <- as.numeric(as.character(ppcs_perf_data2$actual))

roc_obj <- roc(prob,actual)
auc(roc_obj)

pred_roc <- prediction(ppcs_perf_data2$prob, ppcs_perf_data2$actual)
eval <- performance(pred_roc,  measure='tpr', x.measure='fpr')

roc_eval <- data.frame(fpr=unlist(eval@x.values), tpr=unlist(eval@y.values))
ggplot(data=roc_eval, aes(x=fpr, y=tpr)) + 
  theme_bw()+
  geom_line(color = "#451878") +
  geom_abline(linetype=2, color = "#451878") +
  scale_x_continuous(labels=percent, lim=c(0,1)) +
  scale_y_continuous(labels=percent, lim=c(0,1)) +
  xlab('Probability of a False Alarm') +
  ylab('Probability of Detecting Use of Force') +
  labs(title = "ROC Curve of the In-sample Data") +
  theme(axis.title.x = element_text(face = "bold", size=20, color = "#0d693e"),
        axis.title.y = element_text(face = "bold", size=20, color = "#0d693e"),
        axis.text.x = element_text(face = "bold", size=15, hjust = 1),
        axis.text.y = element_text(face = "bold",size=18),
        plot.title = element_text(hjust = 0.5,size=25, color = "#374687",face = "bold"),
        legend.title = element_text(color = "#0d693e", face = "bold", size = 18 ),
        legend.text = element_text(size = 15 ))+
  ggsave("figures/ppcs_roc_curve_no_race.png",  width = 12, height = 10, dpi = 150, units = "in", device='png')
```


## Levels of Force Regressions
# Regression for at least grabbed
added columns according to fryer
```{r}
merged_ppcs <- merged_ppcs %>% mutate(
    at_least_grabbed_fry = ifelse(
      grab_push == 1 |
      handcuffed == 1|
      point_gun == 1 |
      hit_kick == 1  |
      pepper_stun == 1, 1, 0),
    
    at_least_handcuffed_fry = case_when(
     (handcuffed == 1|
      point_gun == 1 |
      hit_kick == 1  |
      pepper_stun == 1) ~ 1,
     (grab_push == 1) ~NA_real_,
     TRUE ~ 0),
    
    at_least_point_gun_fry = case_when(
      (point_gun == 1 |
      hit_kick == 1  |
      pepper_stun == 1) ~ 1,
      (grab_push == 1 | handcuffed == 1) ~ NA_real_,
    TRUE ~ 0),

    at_least_kick_pepper_stun_fry = case_when(
      (hit_kick == 1  |
      pepper_stun == 1) ~ 1,
      (grab_push == 1 | handcuffed == 1 | point_gun == 1) ~ NA_real_,
      TRUE ~ 0),
    
    at_least_grabbed_fry = as.factor(at_least_grabbed_fry),
    at_least_handcuffed_fry = as.factor(at_least_handcuffed_fry),
    at_least_point_gun_fry = as.factor(at_least_point_gun_fry),
    at_least_kick_pepper_stun_fry = as.factor(at_least_kick_pepper_stun_fry)
)

```

# At least Grabbed
5 models for added controls
```{r}
at_least_grab_model1 <- glm(at_least_grabbed_fry ~ civilian_race, family = "binomial", data= merged_ppcs)
exp(coef(at_least_grab_model1))
summary(at_least_grab_model1)
```

```{r}
at_least_grab_model2 <- glm(at_least_grabbed_fry ~ civilian_race + civilian_gender + civilian_employed + civilian_income + population_size + I(civilian_age^2), family = "binomial", data = merged_ppcs)

exp(coef(at_least_grab_model2))
summary(at_least_grab_model2)
```

```{r}
at_least_grab_model3 <- glm(at_least_grabbed_fry ~ civilian_race + civilian_gender + civilian_employed + civilian_income + population_size + I(civilian_age^2) + time_of_encounter + type_of_incident + off_white + off_black + off_hispanic + off_other + off_hispanic + off_split, family = "binomial", data = merged_ppcs)

exp(coef(at_least_grab_model3))
summary(at_least_grab_model3)

summary(merged_ppcs)
```

```{r}
at_least_grab_model4 <- glm(at_least_grabbed_fry  ~ civilian_race + civilian_gender + civilian_employed + civilian_income + population_size + I(civilian_age^2) + time_of_encounter + type_of_incident + off_white + off_black + off_hispanic + off_other + off_hispanic + off_split + civilian_behavior, family = "binomial", data = merged_ppcs)

exp(coef(at_least_grab_model4))
summary(at_least_grab_model4)
```

```{r}

at_least_grab_model5 <- glm(at_least_grabbed ~ civilian_race + civilian_gender + civilian_employed + civilian_income + population_size + I(civilian_age^2) + time_of_encounter + type_of_incident + off_white + off_black + off_hispanic + off_other + off_hispanic + off_split + civilian_behavior + year, family = "binomial", data = merged_ppcs)

exp(coef(at_least_grab_model5))
summary(at_least_grab_model5)
```


#At least handcuffed
```{r}
at_least_handcuffed_model1 <- glm(at_least_handcuffed_fry ~ civilian_race, family = "binomial", data = merged_ppcs)

summary(at_least_handcuffed_model1)
exp(coef(at_least_handcuffed_model1))
```

```{r}
at_least_handcuffed_model2 <- glm(at_least_handcuffed_fry ~ civilian_race+ civilian_gender + civilian_employed + civilian_income + population_size + I(civilian_age^2), family = "binomial", data = merged_ppcs)

summary(at_least_handcuffed_model2)
exp(coef(at_least_handcuffed_model2))
```

```{r}
at_least_handcuffed_model3 <- glm(at_least_handcuffed_fry ~ civilian_race + civilian_gender + civilian_employed + civilian_income + population_size + I(civilian_age^2) + time_of_encounter + type_of_incident + off_white + off_black + off_other + off_hispanic + off_split, family = "binomial", data = merged_ppcs)

exp(coef(at_least_handcuffed_model3))
summary(at_least_handcuffed_model3)
```

```{r}
at_least_handcuffed_model4 <- glm(at_least_handcuffed_fry ~ civilian_race + civilian_gender + civilian_employed + civilian_income + population_size + I(civilian_age^2) + time_of_encounter + type_of_incident + off_white + off_black + off_other + off_hispanic + off_split + civilian_behavior, family = "binomial", data = merged_ppcs)

exp(coef(at_least_handcuffed_model4))
summary(at_least_handcuffed_model4)
```
```{r}
at_least_handcuffed_model5 <- glm(at_least_handcuffed_fry ~ civilian_race + civilian_gender + civilian_employed + civilian_income + population_size + I(civilian_age^2) + time_of_encounter + type_of_incident + off_white + off_black + off_other + off_hispanic + off_split + civilian_behavior + year, family = "binomial", data = merged_ppcs)

exp(coef(at_least_handcuffed_model5))
summary(at_least_handcuffed_model5)
```

# At least point Gun Regressions
```{r}
at_least_point_gun_model1 <- glm(at_least_point_gun_fry ~ civilian_race, family = "binomial", data = merged_ppcs)

summary(at_least_point_gun_model1)
exp(coef(at_least_point_gun_model1))
```

```{r}
at_least_point_gun_model2 <- glm(at_least_point_gun_fry ~ civilian_race + civilian_gender + civilian_employed + civilian_income + population_size + I(civilian_age^2), family = "binomial", data = merged_ppcs)

summary(at_least_point_gun_model2)
exp(coef(at_least_point_gun_model2))
```

```{r}
at_least_point_gun_model3 <- glm(at_least_point_gun ~ civilian_race + civilian_gender + civilian_employed + civilian_income + population_size + I(civilian_age^2) + time_of_encounter + type_of_incident + off_white + off_black + off_other + off_hispanic + off_split, family = "binomial", data = merged_ppcs)

exp(coef(at_least_point_gun_model3))
summary(at_least_point_gun_model3)
```

```{r}
at_least_point_gun_model4 <- glm(at_least_point_gun ~ civilian_race + civilian_gender + civilian_employed + civilian_income + population_size + I(civilian_age^2) + time_of_encounter + type_of_incident + off_white + off_black + off_other + off_hispanic + off_split + civilian_behavior, family = "binomial", data = merged_ppcs)

exp(coef(at_least_point_gun_model4))
summary(at_least_point_gun_model4)
```

```{r}
at_least_point_gun_model5 <- glm(at_least_point_gun_fry ~ civilian_race + civilian_gender + civilian_employed + civilian_income + population_size + I(civilian_age^2) + time_of_encounter + type_of_incident + off_white + off_black + off_other + off_hispanic + off_split + civilian_behavior + year, family = "binomial", data = merged_ppcs)

exp(coef(at_least_point_gun_model5))
summary(at_least_point_gun_model5)

```

# At least kick and pepper spray
```{r}
at_least_kick_pepper_stun_model1 <- glm(at_least_kick_pepper_stun_fry ~ civilian_race, family = "binomial", data = merged_ppcs)

summary(at_least_kick_pepper_stun_model1)
exp(coef(at_least_kick_pepper_stun_model1))
```

```{r}
at_least_kick_pepper_stun_model2 <- glm(at_least_kick_pepper_stun_fry ~ civilian_race + civilian_gender + civilian_employed + civilian_income + population_size + I(civilian_age^2), family = "binomial", data = merged_ppcs)

summary(at_least_kick_pepper_stun_model2)
exp(coef(at_least_kick_pepper_stun_model2))
```

```{r}
at_least_kick_pepper_stun_model3 <- glm(at_least_kick_pepper_stun_fry ~ civilian_race + civilian_gender + civilian_employed + civilian_income + population_size + I(civilian_age^2) + time_of_encounter + type_of_incident + off_white + off_black + off_other + off_hispanic + off_split, family = "binomial", data = merged_ppcs)

exp(coef(at_least_kick_pepper_stun_model3))
summary(at_least_kick_pepper_stun_model3)
```

```{r}
at_least_kick_pepper_stun_model4 <- glm(at_least_kick_pepper_stun_fry ~ civilian_race + civilian_gender + civilian_employed + civilian_income + population_size + I(civilian_age^2) + time_of_encounter + type_of_incident + off_white + off_black + off_other + off_hispanic + off_split + civilian_behavior, family = "binomial", data = merged_ppcs)

exp(coef(at_least_kick_pepper_stun_model4))
summary(at_least_kick_pepper_stun_model4)
```

```{r}
at_least_kick_pepper_stun_model5 <- glm(at_least_kick_pepper_stun_fry ~ civilian_race + civilian_gender + civilian_employed + civilian_income + population_size + I(civilian_age^2) + time_of_encounter + type_of_incident + off_white + off_black + off_other + off_hispanic + off_split + civilian_behavior + year, family = "binomial", data = merged_ppcs)

exp(coef(at_least_kick_pepper_stun_model5))
summary(at_least_kick_pepper_stun_model5)
```