--- subtitle: "homicide and assault analysis" author: "Tyler Lane" date: "25/02/2018" output: word_document editor_options: chunk_output_type: console --- ```{r setup, echo=FALSE} knitr::opts_chunk$set(echo = TRUE) library(tidyverse) library(ggpubr) library(janitor) library(lubridate) library(nlme) library(see) library(zoo) library(forecast) library(scales) library(broom.mixed) # load data df <- read_rds("crime.rds") protested <- read_rds("final protested list.rds") ``` # functions ```{r} # function to transform data for analyses its <- function(x, y) { # create index # homicide data d <- df %>% filter(city == x) %>% select(city, qtr, hom_rate, assault_rate) %>% # add basic ITS terms mutate( level = ifelse(qtr >= y, 1, 0)) %>% group_by(level) %>% mutate(time = ifelse(level == 1, row_number() %>% as.numeric(), rev(row_number() %>% as.numeric() - 1) * -1), trend = ifelse(level == 1, time, 0)) # create index to make homicide and assault proportional i_d <- d %>% filter(qtr %in% c(as.yearqtr("2013 Q3"), as.yearqtr("2013 Q4"), as.yearqtr("2014 Q1"), as.yearqtr("2014 Q2"))) index <- mean(i_d$hom_rate) / mean(i_d$assault_rate) # index rates, finish reformatting d %>% mutate(assault_rate = assault_rate * index) %>% # stacking rates pivot_longer(-c(city, qtr, level, trend, time), names_to = "group", values_to = "rate") %>% mutate( sin1 = sin(2 * pi * 1 * time / 4), sin2 = sin(2 * pi * 2 * time / 4), cos1 = cos(2 * pi * 1 * time / 4), cos2 = cos(2 * pi * 2 * time / 4), group = recode(group, "hom_rate" = "Homicide rate", "assault_rate" = "Assault rate (indexed)") %>% fct_relevel("Homicide rate", "Assault rate (indexed)"), level_assault = ifelse(group == "Assault rate (indexed)", level, 0), trend_assault = ifelse(group == "Assault rate (indexed)", trend, 0), time_assault = ifelse(group == "Assault rate (indexed)", time, 0)) %>% arrange(group) %>% ungroup() %>% return() } # function to plot ITS analyses gg_cities <- function(d, model, event) { d %>% # add fitted datapoints mutate(fit = ifelse(qtr != event, exp(predict(model, d)), NA), fit_ds = ifelse(qtr != event, exp(predict(model, d %>% mutate(sin1 = 0, sin2 = 0, cos1 = 0, cos2 = 0))), NA)) %>% # plot results ggplot(aes(qtr, rate, colour = as.character(group))) + geom_vline(xintercept = as.yearqtr(event), size = 3, alpha = .5, colour = "#95a5a6") + geom_point2(size = 1.3) + scale_colour_manual(values = c("#c83200", "#000080")) + # drawing assault points geom_line(data = . %>% filter(qtr < event, group == "Assault rate (indexed)"), aes(qtr, fit_ds), alpha = .5, size = 1.2, colour = "#c83200") + geom_line(data = . %>% filter(group == "Assault rate (indexed)"), aes(qtr, fit_ds), alpha = .5, size = 1.2, colour = "#c83200") + geom_line(data = . %>% filter(group == "Assault rate (indexed)"), aes(qtr, fit, colour = group), alpha = .5) + # drawing homicide points geom_line(data = . %>% filter(group == "Homicide rate"), aes(qtr, fit_ds), size = 1.2, colour = "#000080") + geom_line(data = . %>% filter(group == "Homicide rate"), aes(qtr, fit, colour = group), alpha = .5) + labs(subtitle = d$city, x = NULL, y = NULL) + theme_lucid() + scale_x_yearqtr(format = "%Y") + coord_cartesian( ylim = c(0, NA)) + theme(legend.title = element_blank(), axis.text = element_text(size = 8), text = element_text(size = 12)) } # autocorrelation check autocor <- function(x) { par(mfrow = c(1, 2)) ggarrange( ggAcf(residuals(x), lag.max = 4) + theme_bw() + labs(title = "ACF", y = NULL, x = NULL), ggPacf(residuals(x), lag.max = 4) + theme_bw() + labs(title = "PACF", y = NULL, x = NULL), ncol = 2) } # function to check models look <- function(x) { x %>% tidy(conf.int = TRUE, conf.level = 0.99) %>% clean_names() %>% transmute(term, `est (99% ci)` = paste0(round(estimate * 100), " (", round(conf_low * 100), ",", round(conf_high * 100), ")"), p = round(p_value, 3), sig = case_when(p_value <= .001 ~ "***", p_value <= .01 ~ "**", p_value <= .05 ~ "*", p_value >= .05 ~ "")) %>% return() } ``` # Albuquerque ```{r} # view city events protested %>% filter(city == "Albuquerque") ``` ## John Edward O'Keefe ```{r} # extract data alb_jeo <- its("Albuquerque", "2015 Q1") # analysis alb_jeo_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = alb_jeo %>% filter(qtr != "2015 Q1")) look(alb_jeo_m) # non-significant seasonal terms excluded alb_jeo_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", data = alb_jeo %>% filter(qtr != "2015 Q1")) look(alb_jeo_m) # autocorrelation check autocor(alb_jeo_m) # adding ARMA terms alb_jeo_p3 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", correlation = corARMA(p = 3, form = ~time | group), data = alb_jeo %>% filter(qtr != "2015 Q1")) alb_jeo_q3 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", correlation = corARMA(q = 3, form = ~time | group), data = alb_jeo %>% filter(qtr != "2015 Q1")) alb_jeo_p3q3 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", correlation = corARMA(p = 3, q = 3, form = ~time | group), data = alb_jeo %>% filter(qtr != "2015 Q1")) # compare model AICs AIC(alb_jeo_m, alb_jeo_p3, alb_jeo_q3, alb_jeo_p3q3) # best fit model outputs look(alb_jeo_p3q3) # assign to standard object for plotting alb_jeo_FINAL <- alb_jeo_p3q3 ``` ## Edgar Camacho-Alvarado ```{r} # extract data alb_eca <- its("Albuquerque", "2016 Q1") # analysis alb_eca_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = alb_eca %>% filter(qtr != "2016 Q1")) look(alb_eca_m) # non-significant seasonal terms excluded alb_eca_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", data = alb_eca %>% filter(qtr != "2016 Q1")) look(alb_eca_m) # autocorrelation check autocor(alb_eca_m) # adding ARMA terms alb_eca_p4 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", correlation = corARMA(p = 4, form = ~time | group), data = alb_eca %>% filter(qtr != "2016 Q1")) alb_eca_q4 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", correlation = corARMA(q = 4, form = ~time | group), data = alb_eca %>% filter(qtr != "2016 Q1")) alb_eca_p4q4 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", correlation = corARMA(p = 4, q = 4, form = ~time | group), data = alb_eca %>% filter(qtr != "2016 Q1")) # compare model AICs AIC(alb_eca_m, alb_eca_p4, alb_eca_q4, alb_eca_p4q4) # best fit model outputs look(alb_eca_q4) # assign to standard object for plotting alb_eca_FINAL <- alb_eca_q4 ``` ## compare events ```{r} AIC(alb_jeo_FINAL, alb_eca_FINAL) alb_jeo_FINAL ``` ## plot ```{r} alb_plot <- gg_cities( d = alb_jeo, model = alb_jeo_FINAL, event = as.yearqtr("2015 Q1")) alb_plot ``` # Anaheim-Santa Ana edge city ```{r} # view city events protested %>% filter(city %in% c("Anaheim", "Santa Ana")) ``` ## Gustavo Najera - Anaheim ```{r} # extract data asa <- its("Anaheim-Santa Ana", "2016 Q1") # analysis asa_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = asa %>% filter(qtr != "2016 Q1")) look(asa_m) # non-significant seasonal terms excluded asa_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group, method = "ML", data = asa %>% filter(qtr != "2016 Q1")) look(asa_m) # autocorrelation check autocor(asa_m) # adding ARMA terms asa_q2 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group, method = "ML", correlation = corARMA(q = 2, form = ~time | group), data = asa %>% filter(qtr != "2016 Q1")) # compare model AICs AIC(asa_m, asa_q2) # best fit model outputs look(asa_q2) # assign to standard object for plotting asa_FINAL <- asa_q2 ``` ## plot ```{r} asa_plot <- gg_cities( d = asa, model = asa_FINAL, event = as.yearqtr("2016 Q1")) asa_plot ``` # Atlanta ```{r} # view city events protested %>% filter(city == "Atlanta") ``` ## Anthony Hill - shooting ```{r} # extract data atl_ah1 <- its("Atlanta", "2015 Q1") %>% # last quarter missing; exclude filter(qtr != "2019 Q4") # analysis atl_ah1_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = atl_ah1 %>% filter(qtr != "2015 Q1")) look(atl_ah1_m) # non-significant seasonal terms excluded atl_ah1_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos2, method = "ML", data = atl_ah1 %>% filter(qtr != "2015 Q1")) look(atl_ah1_m) # autocorrelation check autocor(atl_ah1_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(atl_ah1_m) # assign to standard object for plotting atl_ah1_FINAL <- atl_ah1_m ``` ## Anthony Hill - officer indicted ```{r} # extract data atl_ah2 <- its("Atlanta", "2016 Q1") %>% # last quarter missing; exclude filter(qtr != "2019 Q4") # analysis atl_ah2_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = atl_ah2 %>% filter(qtr != "2016 Q1")) look(atl_ah2_m) # non-significant seasonal terms excluded atl_ah2_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos2, method = "ML", data = atl_ah2 %>% filter(qtr != "2016 Q1")) look(atl_ah2_m) # autocorrelation check autocor(atl_ah2_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(atl_ah2_m) # assign to standard object for plotting atl_ah2_FINAL <- atl_ah2_m ``` ## compare events ```{r} AIC(atl_ah1_FINAL, atl_ah2_FINAL) atl_ah2_FINAL ``` ## plot ```{r} atl_plot <- gg_cities( d = atl_ah2, model = atl_ah2_FINAL, event = as.yearqtr("2016 Q1")) atl_plot ``` # Aurora ```{r} # view city events protested %>% filter(city == "Aurora") ``` ## Naeschylus Vinzant ```{r} # extract data aur <- its("Aurora", "2015 Q1") %>% # add dummy variable to control for movie theatre shooting mutate(movie_shooting = ifelse(qtr == "2012 Q3" & group == "Homicide rate", 1, 0)) # analysis aur_m <- gls( log(rate) ~ level + trend + time + movie_shooting + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = aur %>% filter(qtr != "2015 Q1")) look(aur_m) # non-significant seasonal terms excluded aur_m <- gls( log(rate) ~ level + trend + time + movie_shooting + level_assault + trend_assault + time_assault + group, method = "ML", data = aur %>% filter(qtr != "2015 Q1")) look(aur_m) # autocorrelation check autocor(aur_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(aur_m) # assign to standard object for plotting aur_FINAL <- aur_m ``` ## plot ```{r} aur_plot <- gg_cities( d = aur %>% mutate(movie_shooting = 0), model = aur_FINAL, event = as.yearqtr("2015 Q1")) aur_plot ``` # Austin ```{r} # view city events protested %>% filter(city == "Austin") ``` ## David Joseph ```{r} # extract data aus <- its("Austin", "2016 Q1") # analysis aus_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = aus %>% filter(qtr != "2016 Q1")) look(aus_m) # non-significant seasonal terms excluded aus_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group, method = "ML", data = aus %>% filter(qtr != "2016 Q1")) look(aus_m) # autocorrelation check autocor(aus_m) # adding ARMA terms aus_p2 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group, method = "ML", correlation = corARMA(p = 2, form = ~time | group), data = aus %>% filter(qtr != "2016 Q1")) aus_q2 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group, method = "ML", correlation = corARMA(q = 2, form = ~time | group), data = aus %>% filter(qtr != "2016 Q1")) aus_p2q2 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group, method = "ML", correlation = corARMA(p = 2, q = 2, form = ~time | group), data = aus %>% filter(qtr != "2016 Q1")) # compare model AICs AIC(aus_m, aus_p2, aus_q2, aus_p2q2) # best fit model outputs look(aus_p2q2) # assign to standard object for plotting aus_FINAL <- aus_p2q2 ``` ## plot ```{r} aus_plot <- gg_cities( d = aus, model = aus_FINAL, event = as.yearqtr("2016 Q1")) aus_plot ``` # Bakersfield ```{r} # view city events protested %>% filter(city == "Bakersfield") ``` ## Francisco Serna ```{r} # extract data bak <- its("Bakersfield", "2016 Q4") # analysis bak_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = bak %>% filter(qtr != "2016 Q4")) look(bak_m) # non-significant seasonal terms excluded bak_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1, method = "ML", data = bak %>% filter(qtr != "2016 Q4")) look(bak_m) # autocorrelation check autocor(bak_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(bak_m) # assign to standard object for plotting bak_FINAL <- bak_m ``` ## plot ```{r} bak_plot <- gg_cities( d = bak, model = bak_FINAL, event = as.yearqtr("2016 Q4")) bak_plot ``` # Baltimore ```{r} # view city events protested %>% filter(city == "Baltimore") ``` ## Freddie Gray ```{r} # extract data bal <- its("Baltimore", "2015 Q2") # analysis bal_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = bal %>% filter(qtr != "2015 Q2")) look(bal_m) # non-significant seasonal terms excluded bal_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", data = bal %>% filter(qtr != "2015 Q2")) look(bal_m) # autocorrelation check autocor(bal_m) # adding ARMA terms bal_p3 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", correlation = corARMA(p = 3, form = ~time | group), data = bal %>% filter(qtr != "2015 Q2")) # compare model AICs AIC(bal_m, bal_p3) # best fit model outputs look(bal_m) # assign to standard object for plotting bal_FINAL <- bal_m ``` ## plot ```{r} bal_plot <- gg_cities( d = bal, model = bal_FINAL, event = as.yearqtr("2015 Q2")) bal_plot ``` # Baton Rouge ```{r} # view city events protested %>% filter(city == "Baton Rouge") ``` ## Alton Sterling ```{r} # extract data br <- its("Baton Rouge", "2016 Q3") # analysis br_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = br %>% filter(qtr != "2016 Q3")) look(br_m) # non-significant seasonal terms excluded br_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group, method = "ML", data = br %>% filter(qtr != "2016 Q3")) look(br_m) # autocorrelation check autocor(br_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(br_m) # assign to standard object for plotting br_FINAL <- br_m ``` ## plot ```{r} br_plot <- gg_cities( d = br, model = br_FINAL, event = as.yearqtr("2016 Q3")) br_plot ``` # Charleston ```{r} # view city events protested %>% filter(str_detect(city, "Charleston")) ``` ## Walter Scott ```{r} # extract data cha_sc <- its("Charleston-North Charleston", "2015 Q2") # analysis cha_sc_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = cha_sc %>% filter(qtr != "2015 Q2")) look(cha_sc_m) # non-significant seasonal terms excluded cha_sc_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + cos1, method = "ML", data = cha_sc %>% filter(qtr != "2015 Q2")) look(cha_sc_m) # autocorrelation check autocor(cha_sc_m) # compare model AICs cha_sc_p1 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + cos1, method = "ML", correlation = corARMA(p = 1, form = ~time | group), data = cha_sc %>% filter(qtr != "2015 Q2")) cha_sc_q1 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + cos1, method = "ML", correlation = corARMA(q = 1, form = ~time | group), data = cha_sc %>% filter(qtr != "2015 Q2")) cha_sc_p1q1 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + cos1, method = "ML", correlation = corARMA(p = 1, q = 1, form = ~time | group), data = cha_sc %>% filter(qtr != "2015 Q2")) # compare model AICs AIC(cha_sc_m, cha_sc_p1, cha_sc_q1, cha_sc_p1q1) # best fit model outputs look(cha_sc_p1q1) # assign to standard object for plotting cha_sc_FINAL <- cha_sc_p1q1 ``` ## plot ```{r} cha_sc_plot <- gg_cities( d = cha_sc, model = cha_sc_FINAL, event = as.yearqtr("2015 Q2")) cha_sc_plot ``` # Charlotte ```{r} # view city events protested %>% filter(city == "Charlotte") ``` ## Keith Lamont Scott ```{r} # extract data cha_nc <- its("Charlotte", "2016 Q3") # analysis cha_nc_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = cha_nc %>% filter(qtr != "2016 Q3")) look(cha_nc_m) # non-significant seasonal terms excluded cha_nc_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1, method = "ML", data = cha_nc %>% filter(qtr != "2016 Q3")) look(cha_nc_m) # autocorrelation check autocor(cha_nc_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(cha_nc_m) # assign to standard object for plotting cha_nc_FINAL <- cha_nc_m ``` ## plot ```{r} cha_nc_plot <- gg_cities( d = cha_nc, model = cha_nc_FINAL, event = as.yearqtr("2016 Q3")) cha_nc_plot ``` # Chicago ```{r} # view city events protested %>% filter(str_detect(city, "Chicago")) ``` ## Pierre Loury ```{r} # extract data chi_pl <- its("Chicago", "2016 Q2") # analysis chi_pl_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = chi_pl %>% filter(qtr != "2016 Q2")) look(chi_pl_m) # non-significant seasonal terms excluded chi_pl_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", data = chi_pl %>% filter(qtr != "2016 Q2")) look(chi_pl_m) # autocorrelation check autocor(chi_pl_m) # adding ARMA terms chi_pl_p2 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", correlation = corARMA(p = 2, form = ~time | group), data = chi_pl %>% filter(qtr != "2016 Q2")) chi_pl_q2 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", correlation = corARMA(q = 2, form = ~time | group), data = chi_pl %>% filter(qtr != "2016 Q2")) chi_pl_p2q2 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", correlation = corARMA(p = 2, q = 2, form = ~time | group), data = chi_pl %>% filter(qtr != "2016 Q2")) # compare model AICs AIC(chi_pl_m, chi_pl_p2, chi_pl_q2, chi_pl_p2q2) # best fit model outputs look(chi_pl_p2q2) # assign to standard object for plotting chi_pl_FINAL <- chi_pl_p2q2 ``` ## Paul O'Neal ```{r} # extract data chi_po <- its("Chicago", "2016 Q3") # analysis chi_po_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = chi_po %>% filter(qtr != "2016 Q3")) look(chi_po_m) # non-significant seasonal terms excluded chi_po_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", data = chi_po %>% filter(qtr != "2016 Q3")) look(chi_po_m) # autocorrelation check autocor(chi_po_m) # adding ARMA terms chi_po_p1 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", correlation = corARMA(p = 1, form = ~time | group), data = chi_po %>% filter(qtr != "2016 Q3")) chi_po_q1 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", correlation = corARMA(q = 1, form = ~time | group), data = chi_po %>% filter(qtr != "2016 Q3")) chi_po_p1q1 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", correlation = corARMA(p = 1, q = 1, form = ~time | group), data = chi_po %>% filter(qtr != "2016 Q3")) # compare model AICs AIC(chi_po_m, chi_po_p1, chi_po_q1, chi_po_p1q1) # best fit model outputs look(chi_po_p1) # assign to standard object for plotting chi_po_FINAL <- chi_po_p1 ``` ## Laquan McDonald - release of footage ```{r} # extract data chi_lm <- its("Chicago", "2015 Q4") # analysis chi_lm_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = chi_lm %>% filter(qtr != "2015 Q4")) look(chi_lm_m) # non-significant seasonal terms excluded chi_lm_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1 + cos2, method = "ML", data = chi_lm %>% filter(qtr != "2015 Q4")) look(chi_lm_m) # autocorrelation check autocor(chi_lm_m) # adding ARMA terms chi_lm_p1 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1 + cos2, method = "ML", correlation = corARMA(p = 1, form = ~time | group), data = chi_lm %>% filter(qtr != "2015 Q4")) chi_lm_q1 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1 + cos2, method = "ML", correlation = corARMA(q = 1, form = ~time | group), data = chi_lm %>% filter(qtr != "2015 Q4")) chi_lm_p1q1 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1 + cos2, method = "ML", correlation = corARMA(p = 1, q = 1, form = ~time | group), data = chi_lm %>% filter(qtr != "2015 Q4")) # compare model AICs AIC(chi_lm_m, chi_lm_p1, chi_lm_q1, chi_lm_p1q1) # best fit model outputs look(chi_lm_q1) # assign to standard object for plotting chi_lm_FINAL <- chi_lm_q1 ``` ## compare events ```{r} AIC(chi_pl_FINAL, chi_po_FINAL, chi_lm_FINAL) chi_lm_FINAL ``` ## plot ```{r} chi_plot <- gg_cities( d = chi_lm, model = chi_lm_FINAL, event = as.yearqtr("2015 Q4")) chi_plot ``` # Cincinnati ```{r} # view city events protested %>% filter(str_detect(city, "Cincinnati")) ``` ## QuanDavier Hicks ```{r} # extract data cin_qdh <- its("Cincinnati", "2015 Q2") # analysis cin_qdh_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = cin_qdh %>% filter(qtr != "2015 Q2")) look(cin_qdh_m) # non-significant seasonal terms excluded cin_qdh_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", data = cin_qdh %>% filter(qtr != "2015 Q2")) look(cin_qdh_m) # autocorrelation check autocor(cin_qdh_m) # adding ARMA terms cin_qdh_p4 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, correlation = corARMA(p = 4, form = ~time | group), method = "ML", data = cin_qdh %>% filter(qtr != "2015 Q2")) cin_qdh_q4 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, correlation = corARMA(q = 4, form = ~time | group), method = "ML", data = cin_qdh %>% filter(qtr != "2015 Q2")) cin_qdh_p4q4 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, correlation = corARMA(p = 4, q = 4, form = ~time | group), method = "ML", data = cin_qdh %>% filter(qtr != "2015 Q2")) # compare model AICs AIC(cin_qdh_m, cin_qdh_p4, cin_qdh_q4, cin_qdh_p4q4) # best fit model outputs look(cin_qdh_p4q4) # assign to standard object for plotting cin_qdh_FINAL <- cin_qdh_p4q4 ``` ## Paul Gaston ```{r} # extract data cin_pg <- its("Cincinnati", "2016 Q1") # analysis cin_pg_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = cin_pg %>% filter(qtr != "2016 Q1")) look(cin_pg_m) # non-significant seasonal terms excluded cin_pg_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", data = cin_pg %>% filter(qtr != "2016 Q1")) look(cin_pg_m) # autocorrelation check autocor(cin_pg_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(cin_pg_m) # assign to standard object for plotting cin_pg_FINAL <- cin_pg_m ``` ## Samuel DuBose ```{r} # extract data cin_sd <- its("Cincinnati", "2015 Q3") # analysis cin_sd_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = cin_sd %>% filter(qtr != "2015 Q3")) look(cin_sd_m) # non-significant seasonal terms excluded cin_sd_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", data = cin_sd %>% filter(qtr != "2015 Q3")) look(cin_sd_m) # autocorrelation check autocor(cin_sd_m) # adding ARMA terms cin_sd_p3 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, correlation = corARMA(p = 3, form = ~time | group), method = "ML", data = cin_sd %>% filter(qtr != "2015 Q3")) cin_sd_q3 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, correlation = corARMA(q = 3, form = ~time | group), method = "ML", data = cin_sd %>% filter(qtr != "2015 Q3")) cin_sd_p3q3 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, correlation = corARMA(p = 3, q = 3, form = ~time | group), method = "ML", data = cin_sd %>% filter(qtr != "2015 Q3")) # compare model AICs AIC(cin_sd_m, cin_sd_p3, cin_sd_q3, cin_sd_p3q3) # best fit model outputs look(cin_sd_p3q3) # assign to standard object for plotting cin_sd_FINAL <- cin_sd_p3q3 ``` ## compare events ```{r} AIC(cin_qdh_FINAL, cin_pg_FINAL, cin_sd_FINAL) cin_sd_FINAL ``` ## plot ```{r} cin_plot <- gg_cities( d = cin_sd, model = cin_sd_FINAL, event = as.yearqtr("2015 Q3")) cin_plot ``` # Cleveland ```{r} # view city events protested %>% filter(str_detect(city, "Cleveland")) ``` ## Brandon Jones ```{r} # extract data cle_bj <- its("Cleveland", "2015 Q1") # analysis cle_bj_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = cle_bj %>% filter(qtr != "2015 Q1")) look(cle_bj_m) # non-significant seasonal terms excluded cle_bj_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1, method = "ML", data = cle_bj %>% filter(qtr != "2015 Q1")) look(cle_bj_m) # autocorrelation check autocor(cle_bj_m) # adding ARMA terms cle_bj_p3 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1, method = "ML", correlation = corARMA(p = 3, form = ~time | group), data = cle_bj %>% filter(qtr != "2015 Q1")) cle_bj_q3 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1, method = "ML", correlation = corARMA(q = 3, form = ~time | group), data = cle_bj %>% filter(qtr != "2015 Q1")) cle_bj_p3q3 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1, method = "ML", correlation = corARMA(p = 3, q = 3, form = ~time | group), data = cle_bj %>% filter(qtr != "2015 Q1")) # compare model AICs AIC(cle_bj_m, cle_bj_p3, cle_bj_q3, cle_bj_p3q3) # best fit model outputs look(cle_bj_p3q3) # assign to standard object for plotting cle_bj_FINAL <- cle_bj_p3q3 ``` ## Tamir Rice - shooting ```{r} # extract data cle_tr1 <- its("Cleveland", "2014 Q4") # analysis cle_tr1_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = cle_tr1 %>% filter(qtr != "2014 Q4")) look(cle_tr1_m) # non-significant seasonal terms excluded cle_tr1_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + cos1, method = "ML", data = cle_tr1 %>% filter(qtr != "2014 Q4")) look(cle_tr1_m) # autocorrelation check autocor(cle_tr1_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(cle_tr1_m) # assign to standard object for plotting cle_tr1_FINAL <- cle_tr1_m ``` ## Tamir Rice - officers not indicted ```{r} # extract data cle_tr2 <- its("Cleveland", "2015 Q4") # analysis cle_tr2_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = cle_tr2 %>% filter(qtr != "2015 Q4")) look(cle_tr2_m) # non-significant seasonal terms excluded cle_tr2_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + cos1, method = "ML", data = cle_tr2 %>% filter(qtr != "2015 Q4")) look(cle_tr2_m) # autocorrelation check autocor(cle_tr2_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(cle_tr2_m) # assign to standard object for plotting cle_tr2_FINAL <- cle_tr2_m ``` ## compare events ```{r} AIC(cle_bj_FINAL, cle_tr1_FINAL, cle_tr2_FINAL) cle_tr2_FINAL ``` ## plot ```{r} cle_plot <- gg_cities( d = cle_tr2, model = cle_tr2_FINAL, event = as.yearqtr("2015 Q4")) cle_plot ``` # Columbus ```{r} # view city events protested %>% filter(city == "Columbus") ``` ## Henry Green - following death, demanding investigation ```{r} # extract data col_hg <- its("Columbus", "2016 Q2") # analysis col_hg_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = col_hg %>% filter(qtr != "2016 Q2")) look(col_hg_m) # non-significant seasonal terms excluded col_hg_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", data = col_hg %>% filter(qtr != "2016 Q2")) look(col_hg_m) # autocorrelation check autocor(col_hg_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(col_hg_m) # assign to standard object for plotting col_hg_FINAL <- col_hg_m ``` ## Henry Green/Tyre King - following decision not to indict officer responsible for Green, and following death of King in same quarter ```{r} # extract data col_hgtk <- its("Columbus", "2016 Q3") # analysis col_hgtk_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = col_hgtk %>% filter(qtr != "2016 Q3")) look(col_hgtk_m) # non-significant seasonal terms excluded col_hgtk_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", data = col_hgtk %>% filter(qtr != "2016 Q3")) look(col_hgtk_m) # autocorrelation check autocor(col_hgtk_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(col_hgtk_m) # assign to standard object for plotting col_hgtk_FINAL <- col_hgtk_m ``` ## compare events ```{r} AIC(col_hg_FINAL, col_hgtk_FINAL) col_hgtk_FINAL ``` ## plot ```{r} col_plot <- gg_cities( d = col_hgtk, model = col_hgtk_FINAL, event = as.yearqtr("2016 Q3")) col_plot ``` # Dallas-Forth Worth-Arlington ```{r} # view city events protested %>% filter(city %in% c("Arlington", "Fort Worth", "Grapevine")) ``` ## Christian Taylor (Arlington) - shooting ```{r} # extract data dfwa_ct1 <- its("Dallas-Fort Worth-Arlington", "2015 Q3") # analysis dfwa_ct1_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = dfwa_ct1 %>% filter(qtr != "2015 Q3")) look(dfwa_ct1_m) # non-significant seasonal terms excluded dfwa_ct1_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", data = dfwa_ct1 %>% filter(qtr != "2015 Q3")) look(dfwa_ct1_m) # autocorrelation check autocor(dfwa_ct1_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(dfwa_ct1_m) # assign to standard object for plotting dfwa_ct1_FINAL <- dfwa_ct1_m ``` ## Christian Taylor (Arlington) - officer not indicted ```{r} # extract data dfwa_ct2 <- its("Dallas-Fort Worth-Arlington", "2016 Q2") # analysis dfwa_ct2_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = dfwa_ct2 %>% filter(qtr != "2016 Q2")) look(dfwa_ct2_m) # non-significant seasonal terms excluded dfwa_ct2_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", data = dfwa_ct2 %>% filter(qtr != "2016 Q2")) look(dfwa_ct2_m) # autocorrelation check autocor(dfwa_ct2_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(dfwa_ct2_m) # assign to standard object for plotting dfwa_ct2_FINAL <- dfwa_ct2_m ``` ## Kelvin Goldston (Fort Worth) - shooting/Rubén García Villalpando (Grapevine) - call for release of footage ```{r} # extract data dfwa_kgrgv <- its("Dallas-Fort Worth-Arlington", "2015 Q2") # analysis dfwa_kgrgv_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = dfwa_kgrgv %>% filter(qtr != "2015 Q2")) look(dfwa_kgrgv_m) # non-significant seasonal terms excluded dfwa_kgrgv_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", data = dfwa_kgrgv %>% filter(qtr != "2015 Q2")) look(dfwa_kgrgv_m) # autocorrelation check autocor(dfwa_kgrgv_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(dfwa_kgrgv_m) # assign to standard object for plotting dfwa_kgrgv_FINAL <- dfwa_kgrgv_m ``` ## Rubén García Villalpando (Grapevine) - death ```{r} # extract data dfwa_rgv <- its("Dallas-Fort Worth-Arlington", "2015 Q1") # analysis dfwa_rgv_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = dfwa_rgv %>% filter(qtr != "2015 Q1")) look(dfwa_rgv_m) # non-significant seasonal terms excluded dfwa_rgv_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", data = dfwa_rgv %>% filter(qtr != "2015 Q1")) look(dfwa_rgv_m) # autocorrelation check autocor(dfwa_rgv_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(dfwa_rgv_m) # assign to standard object for plotting dfwa_rgv_FINAL <- dfwa_rgv_m ``` ## compare events ```{r} AIC(dfwa_ct1_FINAL, dfwa_ct2_FINAL, dfwa_kgrgv_FINAL, dfwa_rgv_FINAL) dfwa_ct1_FINAL ``` ## plot ```{r} dfwa_plot <- gg_cities( d = dfwa_ct1, model = dfwa_ct1_FINAL, event = as.yearqtr("2015 Q3")) dfwa_plot ``` # Denver ```{r} # view city events protested %>% filter(city == "Denver") ``` ## Jessie Hernandez ```{r} # extract data den_jh <- its("Detroit", "2015 Q1") # analysis den_jh_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = den_jh %>% filter(qtr != "2015 Q1")) look(den_jh_m) # non-significant seasonal terms excluded den_jh_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", data = den_jh %>% filter(qtr != "2015 Q1")) look(den_jh_m) # autocorrelation check autocor(den_jh_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(den_jh_m) # assign to standard object for plotting den_jh_FINAL <- den_jh_m ``` ## Paul Castaway ```{r} # extract data den_pc <- its("Detroit", "2015 Q3") # analysis den_pc_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = den_pc %>% filter(qtr != "2015 Q3")) look(den_pc_m) # non-significant seasonal terms excluded den_pc_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", data = den_pc %>% filter(qtr != "2015 Q3")) look(den_pc_m) # autocorrelation check autocor(den_pc_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(den_pc_m) # assign to standard object for plotting den_pc_FINAL <- den_pc_m ``` ## Ramone Lonergan ```{r} # extract data den_rl <- its("Detroit", "2016 Q1") # analysis den_rl_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = den_rl %>% filter(qtr != "2016 Q1")) look(den_rl_m) # non-significant seasonal terms excluded den_rl_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", data = den_rl %>% filter(qtr != "2016 Q1")) look(den_rl_m) # autocorrelation check autocor(den_rl_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(den_rl_m) # assign to standard object for plotting den_rl_FINAL <- den_rl_m ``` ## compare events ```{r} AIC(den_jh_FINAL, den_pc_FINAL, den_rl_FINAL) den_rl_FINAL ``` ## plot ```{r} den_plot <- gg_cities( d = den_rl, model = den_rl_FINAL, event = as.yearqtr("2016 Q1")) den_plot ``` # Detroit ```{r} # view city events protested %>% filter(city == "Detroit") ``` ## Terrance Kellom ```{r} # extract data det <- its("Detroit", "2015 Q2") # analysis det_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = det %>% filter(qtr != "2015 Q2")) look(det_m) # non-significant seasonal terms excluded det_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", data = det %>% filter(qtr != "2015 Q2")) look(det_m) # autocorrelation check autocor(det_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(det_m) # assign to standard object for plotting det_FINAL <- det_m ``` ## plot ```{r} det_plot <- gg_cities( d = det, model = det_FINAL, event = as.yearqtr("2015 Q2")) det_plot ``` # Durham ```{r} # view city events protested %>% filter(city == "Durham") ``` ## Frank Nathaniel Clark ```{r} # extract data dur <- its("Durham", "2016 Q4") # analysis dur_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = dur %>% filter(qtr != "2016 Q4")) look(dur_m) # non-significant seasonal terms excluded dur_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group, method = "ML", data = dur %>% filter(qtr != "2016 Q4")) look(dur_m) # autocorrelation check autocor(dur_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(dur_m) # assign to standard object for plotting dur_FINAL <- dur_m ``` ## plot ```{r} dur_plot <- gg_cities( d = dur, model = dur_FINAL, event = as.yearqtr("2016 Q4")) dur_plot ``` # Fresno ```{r} # view city events protested %>% filter(city == "Fresno") ``` ## Freddy Centeno ```{r} # extract data fre_fc <- its("Fresno", "2015 Q3") # analysis fre_fc_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = fre_fc %>% filter(qtr != "2015 Q3")) look(fre_fc_m) # non-significant seasonal terms excluded fre_fc_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + cos1, method = "ML", data = fre_fc %>% filter(qtr != "2015 Q3")) look(fre_fc_m) # autocorrelation check autocor(fre_fc_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(fre_fc_m) # assign to standard object for plotting fre_fc_FINAL <- fre_fc_m ``` ## Dylan Noble ```{r} # extract data fre_dn <- its("Fresno", "2016 Q2") # analysis fre_dn_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = fre_dn %>% filter(qtr != "2016 Q2")) look(fre_dn_m) # non-significant seasonal terms excluded fre_dn_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group, method = "ML", data = fre_dn %>% filter(qtr != "2016 Q2")) look(fre_dn_m) # autocorrelation check autocor(fre_dn_m) # adding ARMA terms fre_dn_p2 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group, method = "ML", correlation = corARMA(p = 2, form = ~time | group), data = fre_dn %>% filter(qtr != "2016 Q2")) fre_dn_q2 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group, method = "ML", correlation = corARMA(q = 2, form = ~time | group), data = fre_dn %>% filter(qtr != "2016 Q2")) fre_dn_p2q2 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group, method = "ML", correlation = corARMA(p = 2, q = 2, form = ~time | group), data = fre_dn %>% filter(qtr != "2016 Q2")) # compare model AICs AIC(fre_dn_m, fre_dn_p2, fre_dn_q2, fre_dn_p2q2) # best fit model outputs look(fre_dn_q2) # assign to standard object for plotting fre_dn_FINAL <- fre_dn_q2 ``` ## compare events ```{r} AIC(fre_fc_FINAL, fre_dn_FINAL) fre_dn_FINAL ``` ## plot ```{r} fre_plot <- gg_cities( d = fre_dn, model = fre_dn_FINAL, event = as.yearqtr("2016 Q2")) fre_plot ``` # Houston ```{r} # view city events protested %>% filter(city %in% c("Hampstead", "Houston")) ``` ## Sandra Bland (Hampstead) ```{r} # extract data hou_sb <- its("Houston", "2015 Q3") # analysis hou_sb_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = hou_sb %>% filter(qtr != "2015 Q3")) look(hou_sb_m) # non-significant seasonal terms excluded hou_sb_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group, method = "ML", data = hou_sb %>% filter(qtr != "2015 Q3")) look(hou_sb_m) # autocorrelation check autocor(hou_sb_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(hou_sb_m) # assign to standard object for plotting hou_sb_FINAL <- hou_sb_m ``` ## Ashtian Barnes - death, little evidence of protest ```{r} # extract data hou_ab1 <- its("Houston", "2016 Q2") # analysis hou_ab1_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = hou_ab1 %>% filter(qtr != "2016 Q2")) look(hou_ab1_m) # non-significant seasonal terms excluded hou_ab1_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + cos1, method = "ML", data = hou_ab1 %>% filter(qtr != "2016 Q2")) look(hou_ab1_m) # autocorrelation check autocor(hou_ab1_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(hou_ab1_m) # assign to standard object for plotting hou_ab1_FINAL <- hou_ab1_m ``` ## Ashtian Barnes - protest followed acquittal of officers ```{r} # extract data hou_ab2 <- its("Houston", "2016 Q3") # analysis hou_ab2_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = hou_ab2 %>% filter(qtr != "2016 Q3")) look(hou_ab2_m) # non-significant seasonal terms excluded hou_ab2_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group, method = "ML", data = hou_ab2 %>% filter(qtr != "2016 Q3")) look(hou_ab2_m) # autocorrelation check autocor(hou_ab2_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(hou_ab2_m) # assign to standard object for plotting hou_ab2_FINAL <- hou_ab2_m ``` ## compare events ```{r} AIC(hou_sb_FINAL, hou_ab1_FINAL, hou_ab2_FINAL) hou_sb_FINAL ``` ## plot ```{r} hou_plot <- gg_cities( d = hou_sb, model = hou_sb_FINAL, event = as.yearqtr("2015 Q3")) hou_plot ``` # Indianapolis ```{r} # view city events protested %>% filter(city == "Indianapolis") ``` ## Mack Long - death ```{r} # extract data ind_ml1 <- its("Indianapolis", "2015 Q2") # analysis ind_ml1_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = ind_ml1 %>% filter(qtr != "2015 Q2")) look(ind_ml1_m) # non-significant seasonal terms excluded ind_ml1_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", data = ind_ml1 %>% filter(qtr != "2015 Q2")) look(ind_ml1_m) # autocorrelation check autocor(ind_ml1_m) # adding ARMA terms ind_ml1_q2 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", correlation = corARMA(q = 2, form = ~time | group), data = ind_ml1 %>% filter(qtr != "2015 Q2")) # compare model AICs AIC(ind_ml1_m, ind_ml1_q2) # best fit model outputs look(ind_ml1_q2) # assign to standard object for plotting ind_ml1_FINAL <- ind_ml1_q2 ``` ## Mack Long - demand for footage release/Andre Green ```{r} # extract data ind_ml2 <- its("Indianapolis", "2015 Q3") # analysis ind_ml2_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = ind_ml2 %>% filter(qtr != "2015 Q3")) look(ind_ml2_m) # non-significant seasonal terms excluded ind_ml2_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", data = ind_ml2 %>% filter(qtr != "2015 Q3")) look(ind_ml2_m) # autocorrelation check autocor(ind_ml2_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs # assign to standard object for plotting ind_ml2_FINAL <- ind_ml2_m ``` ## compare events ```{r} AIC(ind_ml1_FINAL, ind_ml2_FINAL) ind_ml1_FINAL ``` ## plot ```{r} ind_plot <- gg_cities( d = ind_ml1, model = ind_ml1_FINAL, event = as.yearqtr("2015 Q2")) ind_plot ``` # Long Beach ```{r} # view city events protested %>% filter(city == "Long Beach") ``` ## Feras Morad ```{r} # extract data lb_fm <- its("Long Beach", "2015 Q2") # analysis lb_fm_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = lb_fm %>% filter(qtr != "2015 Q2")) look(lb_fm_m) # non-significant seasonal terms excluded lb_fm_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + cos1, method = "ML", data = lb_fm %>% filter(qtr != "2015 Q2")) look(lb_fm_m) # autocorrelation check autocor(lb_fm_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(lb_fm_m) # assign to standard object for plotting lb_fm_FINAL <- lb_fm_m ``` ## Barry Prak ```{r} # extract data lb_bp <- its("Long Beach", "2016 Q2") # analysis lb_bp_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = lb_bp %>% filter(qtr != "2016 Q2")) look(lb_bp_m) # non-significant seasonal terms excluded lb_bp_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + cos1, method = "ML", data = lb_bp %>% filter(qtr != "2016 Q2")) look(lb_bp_m) # autocorrelation check autocor(lb_bp_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(lb_bp_m) # assign to standard object for plotting lb_bp_FINAL <- lb_bp_m ``` ## compare events ```{r} AIC(lb_fm_FINAL, lb_bp_FINAL) lb_fm_FINAL ``` ## plot ```{r} lb_plot <- gg_cities( d = lb_fm, model = lb_fm_FINAL, event = as.yearqtr("2015 Q2")) lb_plot ``` # Los Angeles ```{r} # view city events protested %>% filter(city == "Los Angeles") ``` ## Charly Leundeu Keunang ```{r} # extract data la_clk <- its("Los Angeles", "2015 Q1") # analysis la_clk_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = la_clk %>% filter(qtr != "2015 Q1")) look(la_clk_m) # non-significant seasonal terms excluded la_clk_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", data = la_clk %>% filter(qtr != "2015 Q1")) look(la_clk_m) # autocorrelation check autocor(la_clk_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(la_clk_m) # assign to standard object for plotting la_clk_FINAL <- la_clk_m ``` ## Redel Jones - death ```{r} # extract data la_rj1 <- its("Los Angeles", "2015 Q3") # analysis la_rj1_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = la_rj1 %>% filter(qtr != "2015 Q3")) look(la_rj1_m) # non-significant seasonal terms excluded la_rj1_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", data = la_rj1 %>% filter(qtr != "2015 Q3")) look(la_rj1_m) # autocorrelation check autocor(la_rj1_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(la_rj1_m) # assign to standard object for plotting la_rj1_FINAL <- la_rj1_m ``` ## Redel Jones - involved officer not indicted ```{r} # extract data la_rj2 <- its("Los Angeles", "2016 Q2") # analysis la_rj2_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = la_rj2 %>% filter(qtr != "2016 Q2")) look(la_rj2_m) # non-significant seasonal terms excluded la_rj2_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", data = la_rj2 %>% filter(qtr != "2016 Q2")) look(la_rj2_m) # autocorrelation check autocor(la_rj2_m) # adding ARMA terms la_rj2_p1 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", correlation = corARMA(p = 1, form = ~time | group), data = la_rj2 %>% filter(qtr != "2016 Q2")) la_rj2_q1 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", correlation = corARMA(q = 1, form = ~time | group), data = la_rj2 %>% filter(qtr != "2016 Q2")) la_rj2_p1q1 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", correlation = corARMA(p = 1, q = 1, form = ~time | group), data = la_rj2 %>% filter(qtr != "2016 Q2")) # compare model AICs AIC(la_rj2_m, la_rj2_p1, la_rj2_q1, la_rj2_p1q1) # best fit model outputs look(la_rj2_p1q1) # assign to standard object for plotting la_rj2_FINAL <- la_rj2_p1q1 ``` ## Kenney Watkins ```{r} # extract data la_kw <- its("Los Angeles", "2016 Q3") # analysis la_kw_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = la_kw %>% filter(qtr != "2016 Q3")) look(la_kw_m) # non-significant seasonal terms excluded la_kw_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", data = la_kw %>% filter(qtr != "2016 Q3")) look(la_kw_m) # autocorrelation check autocor(la_kw_m) # adding ARMA terms la_kw_p1 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", correlation = corARMA(p = 1, form = ~time | group), data = la_kw %>% filter(qtr != "2016 Q3")) la_kw_q1 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", correlation = corARMA(q = 1, form = ~time | group), data = la_kw %>% filter(qtr != "2016 Q3")) la_kw_p1q1 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", correlation = corARMA(p = 1, q = 1, form = ~time | group), data = la_kw %>% filter(qtr != "2016 Q3")) # compare model AICs AIC(la_kw_m, la_kw_p1, la_kw_q1, la_kw_p1q1) # best fit model outputs look(la_kw_p1) # assign to standard object for plotting la_kw_FINAL <- la_kw_p1 ``` ## Carnell Snell ```{r} # extract data la_cs <- its("Los Angeles", "2016 Q4") # analysis la_cs_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = la_cs %>% filter(qtr != "2016 Q4")) look(la_cs_m) # non-significant seasonal terms excluded la_cs_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", data = la_cs %>% filter(qtr != "2016 Q4")) look(la_cs_m) # autocorrelation check autocor(la_cs_m) # adding ARMA terms la_cs_p1 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", correlation = corARMA(p = 1, form = ~time | group), data = la_cs %>% filter(qtr != "2016 Q4")) la_cs_q1 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", correlation = corARMA(q = 1, form = ~time | group), data = la_cs %>% filter(qtr != "2016 Q4")) la_cs_p1q1 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", correlation = corARMA(p = 1, q = 1, form = ~time | group), data = la_cs %>% filter(qtr != "2016 Q4")) # compare model AICs AIC(la_cs_m, la_cs_p1, la_cs_q1, la_cs_p1q1) # best fit model outputs look(la_cs_p1) # assign to standard object for plotting la_cs_FINAL <- la_cs_p1 ``` ## Ezell Ford ```{r} # extract data la_ef <- its("Los Angeles", "2014 Q3") # analysis la_ef_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = la_ef %>% filter(qtr != "2014 Q3")) look(la_ef_m) # non-significant seasonal terms excluded la_ef_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", data = la_ef %>% filter(qtr != "2014 Q3")) look(la_ef_m) # autocorrelation check autocor(la_ef_m) # adding ARMA terms la_ef_p2 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", correlation = corARMA(p = 2, form = ~time | group), data = la_ef %>% filter(qtr != "2014 Q3")) la_ef_q2 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", correlation = corARMA(q = 2, form = ~time | group), data = la_ef %>% filter(qtr != "2014 Q3")) la_ef_p2q2 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", correlation = corARMA(p = 2, q = 2, form = ~time | group), data = la_ef %>% filter(qtr != "2014 Q3")) # compare model AICs AIC(la_ef_m, la_ef_p2, la_ef_q2, la_ef_p2q2) # best fit model outputs look(la_ef_p2q2) # assign to standard object for plotting la_ef_FINAL <- la_ef_p2q2 ``` ## compare events ```{r} AIC(la_clk_FINAL, la_rj1_FINAL, la_rj2_FINAL, la_kw_FINAL, la_cs_FINAL, la_ef_FINAL) la_ef_FINAL ``` ## plot ```{r} la_plot <- gg_cities( d = la_ef, model = la_ef_FINAL, event = as.yearqtr("2014 Q3")) la_plot ``` # Louisville ```{r} # view city events protested %>% filter(city == "Louisville") ``` ## Deng Manyoun ```{r} # extract data lou <- its("Louisville", "2015 Q2") # analysis lou_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = lou %>% filter(qtr != "2015 Q2")) look(lou_m) # non-significant seasonal terms excluded lou_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + cos1, method = "ML", data = lou %>% filter(qtr != "2015 Q2")) look(lou_m) # autocorrelation check autocor(lou_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(lou_m) # assign to standard object for plotting lou_FINAL <- lou_m ``` ## plot ```{r} lou_plot <- gg_cities( d = lou, model = lou_FINAL, event = as.yearqtr("2015 Q1")) lou_plot ``` # Memphis ```{r} # view city events protested %>% filter(city == "Memphis") ``` ## Darrius Stewart - death ```{r} # extract data mem_ds1 <- its("Memphis", "2015 Q3") # analysis mem_ds1_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = mem_ds1 %>% filter(qtr != "2015 Q3")) look(mem_ds1_m) # non-significant seasonal terms excluded mem_ds1_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1, method = "ML", data = mem_ds1 %>% filter(qtr != "2015 Q3")) look(mem_ds1_m) # autocorrelation check autocor(mem_ds1_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(mem_ds1_m) # assign to standard object for plotting mem_ds1_FINAL <- mem_ds1_m ``` ## Darrius Stewart - officer not indicted ```{r} # extract data mem_ds2 <- its("Memphis", "2015 Q4") # analysis mem_ds2_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = mem_ds2 %>% filter(qtr != "2015 Q4")) look(mem_ds2_m) # non-significant seasonal terms excluded mem_ds2_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + cos1, method = "ML", data = mem_ds2 %>% filter(qtr != "2015 Q4")) look(mem_ds2_m) # autocorrelation check autocor(mem_ds2_m) # adding ARMA terms mem_ds2_p3 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + cos1, method = "ML", correlation = corARMA(p = 3, form = ~time | group), data = mem_ds2 %>% filter(qtr != "2015 Q4")) mem_ds2_q3 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + cos1, method = "ML", correlation = corARMA(q = 3, form = ~time | group), data = mem_ds2 %>% filter(qtr != "2015 Q4")) mem_ds2_p3q3 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + cos1, method = "ML", correlation = corARMA(p = 3, q = 3, form = ~time | group), data = mem_ds2 %>% filter(qtr != "2015 Q4")) # compare model AICs AIC(mem_ds2_m, mem_ds2_p3, mem_ds2_q3, mem_ds2_p3q3) # best fit model outputs look(mem_ds2_q3) # assign to standard object for plotting mem_ds2_FINAL <- mem_ds2_q3 ``` ## compare events ```{r} AIC(mem_ds1_FINAL, mem_ds2_FINAL) mem_ds2_FINAL ``` ## plot ```{r} mem_plot <- gg_cities( d = mem_ds2, model = mem_ds2_FINAL, event = as.yearqtr("2015 Q4")) mem_plot ``` # Milwaukee ```{r} # view city events protested %>% filter(city == "Milwaukee") ``` ## Sylville Smith ```{r} # extract data mil <- its("Milwaukee", "2016 Q3") # analysis mil_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = mil %>% filter(qtr != "2016 Q3")) look(mil_m) # non-significant seasonal terms excluded mil_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", data = mil %>% filter(qtr != "2016 Q3")) look(mil_m) # autocorrelation check autocor(mil_m) # adding ARMA terms mil_p2 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", correlation = corARMA(p = 2, form = ~time | group), data = mil %>% filter(qtr != "2016 Q3")) mil_q4 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", correlation = corARMA(q = 4, form = ~time | group), data = mil %>% filter(qtr != "2016 Q3")) # mil_p2q4 <- gls( # log(rate) ~ level + trend + time + # level_assault + trend_assault + time_assault + group + # sin1 + cos1, # method = "ML", # correlation = corARMA(p = 2, q = 4, form = ~time | group), # data = mil %>% # filter(qtr != "2016 Q3")) # compare model AICs AIC(mil_m, mil_p2, mil_q4) # best fit model outputs look(mil_q4) # assign to standard object for plotting mil_FINAL <- mil_q4 ``` ## plot ```{r} mil_plot <- gg_cities( d = mil, model = mil_FINAL, event = as.yearqtr("2016 Q3")) mil_plot ``` # Minneapolis-St Paul ```{r} # view city events protested %>% filter(city %in% c("Falcon Heights", "Minneapolis", "St. Paul")) ``` ## Philando Castile (Falcon Heights) ```{r} # extract data msp_pc <- its("Minneapolis-St Paul", "2016 Q3") %>% filter(year(qtr) >= 2013) # analysis msp_pc_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = msp_pc %>% filter(qtr != "2016 Q3")) look(msp_pc_m) # non-significant seasonal terms excluded msp_pc_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1, method = "ML", data = msp_pc %>% filter(qtr != "2016 Q3")) look(msp_pc_m) # autocorrelation check autocor(msp_pc_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(msp_pc_m) # assign to standard object for plotting msp_pc_FINAL <- msp_pc_m ``` ## Jamar Clark (Minneapolis) ```{r} # extract data msp_jc <- its("Minneapolis-St Paul", "2015 Q4") %>% filter(year(qtr) >= 2013) # analysis msp_jc_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = msp_jc %>% filter(qtr != "2015 Q4")) look(msp_jc_m) # non-significant seasonal terms excluded msp_jc_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + cos1, method = "ML", data = msp_jc %>% filter(qtr != "2015 Q4")) look(msp_jc_m) # autocorrelation check autocor(msp_jc_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(msp_jc_m) # assign to standard object for plotting msp_jc_FINAL <- msp_jc_m ``` ## Philip Quinn (St Paul) - death ```{r} # extract data msp_pq1 <- its("Minneapolis-St Paul", "2015 Q3") %>% filter(year(qtr) >= 2013) # analysis msp_pq1_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = msp_pq1 %>% filter(qtr != "2015 Q3")) look(msp_pq1_m) # non-significant seasonal terms excluded msp_pq1_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1, method = "ML", data = msp_pq1 %>% filter(qtr != "2015 Q3")) look(msp_pq1_m) # autocorrelation check autocor(msp_pq1_m) # adding ARMA terms msp_pq1_p2 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1, method = "ML", correlation = corARMA(p = 2, form = ~time | group), data = msp_pq1 %>% filter(qtr != "2015 Q3")) msp_pq1_q2 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1, method = "ML", correlation = corARMA(q = 2, form = ~time | group), data = msp_pq1 %>% filter(qtr != "2015 Q3")) msp_pq1_p2q2 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1, method = "ML", correlation = corARMA(p = 2, q = 2, form = ~time | group), data = msp_pq1 %>% filter(qtr != "2015 Q3")) # compare model AICs AIC(msp_pq1_m, msp_pq1_p2, msp_pq1_q2, msp_pq1_p2q2) # best fit model outputs look(msp_pq1_m) # assign to standard object for plotting msp_pq1_FINAL <- msp_pq1_m ``` ## Philip Quinn (St Paul) ```{r} # extract data msp_pq2 <- its("Minneapolis-St Paul", "2016 Q1") %>% filter(year(qtr) >= 2013) # analysis msp_pq2_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = msp_pq2 %>% filter(qtr != "2016 Q1")) look(msp_pq2_m) # non-significant seasonal terms excluded msp_pq2_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1, method = "ML", data = msp_pq2 %>% filter(qtr != "2016 Q1")) look(msp_pq2_m) # autocorrelation check autocor(msp_pq2_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(msp_pq2_m) # assign to standard object for plotting msp_pq2_FINAL <- msp_pq2_m ``` ## compare events ```{r} AIC(msp_pc_FINAL, msp_jc_FINAL, msp_pq1_FINAL, msp_pq2_FINAL) msp_jc_FINAL ``` ## plot ```{r} msp_plot <- gg_cities( d = msp_jc, model = msp_jc_FINAL, event = as.yearqtr("2015 Q4")) msp_plot ``` # New Orleans ```{r} # view city events protested %>% filter(city == "New Orleans") ``` ## Eric Harris ```{r} # extract data no <- its("New Orleans", "2016 Q1") # analysis no_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = no %>% filter(qtr != "2016 Q1")) look(no_m) # non-significant seasonal terms excluded no_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group, method = "ML", data = no %>% filter(qtr != "2016 Q1")) look(no_m) # autocorrelation check autocor(no_m) # adding ARMA terms no_p2 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group, method = "ML", correlation = corARMA(p = 2, form = ~time | group), data = no %>% filter(qtr != "2016 Q1")) no_q2 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group, method = "ML", correlation = corARMA(q = 2, form = ~time | group), data = no %>% filter(qtr != "2016 Q1")) no_p2q2 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group, method = "ML", correlation = corARMA(p = 2, q = 2, form = ~time | group), data = no %>% filter(qtr != "2016 Q1")) # compare model AICs AIC(no_m, no_p2, no_q2, no_p2q2) # best fit model outputs look(no_q2) # assign to standard object for plotting no_FINAL <- no_q2 ``` ## plot ```{r} no_plot <- gg_cities( d = no, model = no_FINAL, event = as.yearqtr("2016 Q1")) no_plot ``` # New York ```{r} # view city events protested %>% filter(city == "New York") ``` ## Eric Garner - death ```{r} # extract data ny_eg1 <- its("New York", "2014 Q3") # analysis ny_eg1_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = ny_eg1 %>% filter(qtr != "2014 Q3")) look(ny_eg1_m) # non-significant seasonal terms excluded ny_eg1_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1, method = "ML", data = ny_eg1 %>% filter(qtr != "2014 Q3")) look(ny_eg1_m) # autocorrelation check autocor(ny_eg1_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(ny_eg1_m) # assign to standard object for plotting ny_eg1_FINAL <- ny_eg1_m ``` ## Eric Garner - officer acquitted ```{r} # extract data ny_eg2 <- its("New York", "2014 Q4") # analysis ny_eg2_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = ny_eg2 %>% filter(qtr != "2014 Q4")) look(ny_eg2_m) # non-significant seasonal terms excluded ny_eg2_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + cos1, method = "ML", data = ny_eg2 %>% filter(qtr != "2014 Q4")) look(ny_eg2_m) # autocorrelation check autocor(ny_eg2_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(ny_eg2_m) # assign to standard object for plotting ny_eg2_FINAL <- ny_eg2_m ``` ## Deborah Danner ```{r} # extract data ny_dd <- its("New York", "2016 Q4") # analysis ny_dd_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = ny_dd %>% filter(qtr != "2016 Q4")) look(ny_dd_m) # non-significant seasonal terms excluded ny_dd_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + cos1, method = "ML", data = ny_dd %>% filter(qtr != "2016 Q4")) look(ny_dd_m) # autocorrelation check autocor(ny_dd_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(ny_dd_m) # assign to standard object for plotting ny_dd_FINAL <- ny_dd_m ``` ## compare events ```{r} AIC(ny_eg1_FINAL, ny_eg2_FINAL, ny_dd_FINAL) ny_eg2_FINAL ``` ## plot ```{r} ny_plot <- gg_cities( d = ny_eg2, model = ny_eg2_FINAL, event = as.yearqtr("2014 Q4")) ny_plot ``` # Oakland ```{r} # view city events protested %>% filter(city == "Oakland") ``` ## Demouria Hogg ```{r} # extract data oak_dh <- its("Oakland", "2015 Q2") # analysis oak_dh_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = oak_dh %>% filter(qtr != "2016 Q2")) look(oak_dh_m) # non-significant seasonal terms excluded oak_dh_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group, method = "ML", data = oak_dh %>% filter(qtr != "2016 Q2")) look(oak_dh_m) # autocorrelation check autocor(oak_dh_m) # adding ARMA terms oak_dh_p2 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group, method = "ML", correlation = corARMA(p = 2, form = ~time | group), data = oak_dh %>% filter(qtr != "2016 Q2")) oak_dh_q2 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group, method = "ML", correlation = corARMA(q = 2, form = ~time | group), data = oak_dh %>% filter(qtr != "2016 Q2")) oak_dh_p2q2 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group, method = "ML", correlation = corARMA(p = 2, q = 2, form = ~time | group), data = oak_dh %>% filter(qtr != "2016 Q2")) # compare model AICs AIC(oak_dh_m, oak_dh_p2, oak_dh_q2, oak_dh_p2q2) # best fit model outputs look(oak_dh_p2q2) # assign to standard object for plotting oak_dh_FINAL <- oak_dh_p2q2 ``` ## Nathaniel Wilks ```{r} # extract data oak_nw <- its("Oakland", "2015 Q3") # analysis oak_nw_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = oak_nw %>% filter(qtr != "2015 Q3")) look(oak_nw_m) # non-significant seasonal terms excluded oak_nw_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group, method = "ML", data = oak_nw %>% filter(qtr != "2015 Q3")) look(oak_nw_m) # autocorrelation check autocor(oak_nw_m) # adding ARMA terms oak_nw_p4 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group, method = "ML", correlation = corARMA(p = 4, form = ~time | group), data = oak_nw %>% filter(qtr != "2015 Q3")) oak_nw_q4 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group, method = "ML", correlation = corARMA(q = 4, form = ~time | group), data = oak_nw %>% filter(qtr != "2015 Q3")) oak_nw_p4q4 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group, method = "ML", correlation = corARMA(p = 4, q = 4, form = ~time | group), data = oak_nw %>% filter(qtr != "2015 Q3")) # compare model AICs AIC(oak_nw_m, oak_nw_p4, oak_nw_q4, oak_nw_p4q4) # best fit model outputs look(oak_nw_p4q4) # assign to standard object for plotting oak_nw_FINAL <- oak_nw_p4q4 ``` ## compare events ```{r} AIC(oak_dh_FINAL, oak_nw_FINAL) oak_dh_FINAL ``` ## plot ```{r} oak_plot <- gg_cities( d = oak_dh, model = oak_dh_FINAL, event = as.yearqtr("2015 Q2")) oak_plot ``` # Philadelphia ```{r} # view city events protested %>% filter(city == "Philadelphia") ``` ## Christopher Sowell - protests after a few days, slipped into next quarter ```{r} # extract data phi_cs1 <- its("Philadelphia", "2016 Q3") # analysis phi_cs1_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = phi_cs1 %>% filter(qtr != "2016 Q3")) look(phi_cs1_m) # non-significant seasonal terms excluded phi_cs1_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", data = phi_cs1 %>% filter(qtr != "2016 Q3")) look(phi_cs1_m) # autocorrelation check autocor(phi_cs1_m) # adding ARMA terms phi_cs1_p1 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, correlation = corARMA(p = 1, form = ~time | group), method = "ML", data = phi_cs1 %>% filter(qtr != "2016 Q3")) phi_cs1_q1 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, correlation = corARMA(q = 1, form = ~time | group), method = "ML", data = phi_cs1 %>% filter(qtr != "2016 Q3")) phi_cs1_p1q1 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, correlation = corARMA(p = 1, q = 1, form = ~time | group), method = "ML", data = phi_cs1 %>% filter(qtr != "2016 Q3")) # compare model AICs AIC(phi_cs1_m, phi_cs1_p1, phi_cs1_q1, phi_cs1_p1q1) # best fit model outputs look(phi_cs1_q1) # assign to standard object for plotting phi_cs1_FINAL <- phi_cs1_q1 ``` ## Christopher Sowell - protests after a few days, slipped into next quarter ```{r} # extract data phi_cs2 <- its("Philadelphia", "2016 Q4") # analysis phi_cs2_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = phi_cs2 %>% filter(qtr != "2016 Q4")) look(phi_cs2_m) # non-significant seasonal terms excluded phi_cs2_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + cos1, method = "ML", data = phi_cs2 %>% filter(qtr != "2016 Q4")) look(phi_cs2_m) # autocorrelation check autocor(phi_cs2_m) # adding ARMA terms phi_cs2_p1 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + cos1, correlation = corARMA(p = 1, form = ~time | group), method = "ML", data = phi_cs2 %>% filter(qtr != "2016 Q4")) phi_cs2_q1 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + cos1, correlation = corARMA(q = 1, form = ~time | group), method = "ML", data = phi_cs2 %>% filter(qtr != "2016 Q4")) phi_cs2_p1q1 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + cos1, correlation = corARMA(p = 1, q = 1, form = ~time | group), method = "ML", data = phi_cs2 %>% filter(qtr != "2016 Q4")) # compare model AICs AIC(phi_cs2_m, phi_cs2_p1, phi_cs2_q1, phi_cs2_p1q1) # best fit model outputs look(phi_cs2_q1) # assign to standard object for plotting phi_cs2_FINAL <- phi_cs2_q1 ``` ## Brandon Tate-Brown - death ```{r} # extract data phi_btb1 <- its("Philadelphia", "2014 Q4") # analysis phi_btb1_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = phi_btb1 %>% filter(qtr != "2014 Q4")) look(phi_btb1_m) # non-significant seasonal terms excluded phi_btb1_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + cos1, method = "ML", data = phi_btb1 %>% filter(qtr != "2014 Q4")) look(phi_btb1_m) # autocorrelation check autocor(phi_btb1_m) # adding ARMA terms phi_btb1_p1 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + cos1, method = "ML", correlation = corARMA(p = 1, form = ~time | group), data = phi_btb1 %>% filter(qtr != "2016 Q3")) phi_btb1_q2 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + cos1, method = "ML", correlation = corARMA(p = 2, form = ~time | group), data = phi_btb1 %>% filter(qtr != "2016 Q3")) # phi_btb1_p1q2 <- gls( # log(rate) ~ level + trend + time + # level_assault + trend_assault + time_assault + group + # cos1, # method = "ML", # correlation = corARMA(p = 1, q = 2, form = ~time | group), # data = phi_btb1 %>% # filter(qtr != "2016 Q3")) # compare model AICs AIC(phi_btb1_m, phi_btb1_p1, phi_btb1_q2) # best fit model outputs look(phi_btb1_q2) # assign to standard object for plotting phi_btb1_FINAL <- phi_btb1_q2 ``` ## Brandon Tate-Brown - charges against officer dropped ```{r} # extract data phi_btb2 <- its("Philadelphia", "2015 Q1") # analysis phi_btb2_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = phi_btb2 %>% filter(qtr != "2015 Q1")) look(phi_btb2_m) # non-significant seasonal terms excluded phi_btb2_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1, method = "ML", data = phi_btb2 %>% filter(qtr != "2015 Q1")) look(phi_btb2_m) # autocorrelation check autocor(phi_btb2_m) # adding ARMA terms phi_btb2_p1 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1, method = "ML", correlation = corARMA(p = 1, form = ~time | group), data = phi_btb2 %>% filter(qtr != "2015 Q1")) phi_btb2_q1 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1, method = "ML", correlation = corARMA(q = 1, form = ~time | group), data = phi_btb2 %>% filter(qtr != "2015 Q1")) phi_btb2_p1q1 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1, method = "ML", correlation = corARMA(p = 1, q = 1, form = ~time | group), data = phi_btb2 %>% filter(qtr != "2015 Q1")) # compare model AICs AIC(phi_btb2_m, phi_btb2_p1, phi_btb2_q1, phi_btb2_p1q1) # best fit model outputs look(phi_btb2_q1) # assign to standard object for plotting phi_btb2_FINAL <- phi_btb2_q1 ``` ## compare events ```{r} AIC(phi_cs1_FINAL, phi_cs2_FINAL, phi_btb1_FINAL, phi_btb2_FINAL) phi_cs1_FINAL ``` ## plot ```{r} phi_plot <- gg_cities( d = phi_cs1, model = phi_cs1_FINAL, event = as.yearqtr("2016 Q3")) phi_plot ``` # Greater Phoenix ```{r} # view city events protested %>% filter(city %in% c("Mesa", "Tempe")) ``` ## Daniel Shaver/Kayden Clarke (both Mesa) ```{r} # extract data pho_dskc <- its("Greater Phoenix", "2016 Q1") # analysis pho_dskc_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = pho_dskc %>% filter(qtr != "2016 Q1")) look(pho_dskc_m) # non-significant seasonal terms excluded pho_dskc_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos2, method = "ML", data = pho_dskc %>% filter(qtr != "2016 Q1")) look(pho_dskc_m) # autocorrelation check autocor(pho_dskc_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(pho_dskc_m) # assign to standard object for plotting pho_dskc_FINAL <- pho_dskc_m ``` ## Dalvin Hollins ```{r} # extract data pho_dh <- its("Greater Phoenix", "2016 Q3") # analysis pho_dh_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = pho_dh %>% filter(qtr != "2016 Q3")) look(pho_dh_m) # non-significant seasonal terms excluded pho_dh_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos2, method = "ML", data = pho_dh %>% filter(qtr != "2016 Q3")) look(pho_dh_m) # autocorrelation check autocor(pho_dh_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(pho_dh_m) # assign to standard object for plotting pho_dh_FINAL <- pho_dh_m ``` ## compare events ```{r} AIC(pho_dskc_FINAL, pho_dh_FINAL) pho_dskc_FINAL ``` ## plot ```{r} pho_plot <- gg_cities( d = pho_dskc, model = pho_dskc_FINAL, event = as.yearqtr("2016 Q1")) pho_plot ``` # Raleigh ```{r} # view city events protested %>% filter(city == "Raleigh") ``` ## Akiel Denkins ```{r} # extract data ral <- its("Raleigh", "2016 Q1") # analysis ral_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = ral %>% filter(qtr != "2016 Q1")) look(ral_m) # non-significant seasonal terms excluded ral_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1, method = "ML", data = ral %>% filter(qtr != "2016 Q1")) look(ral_m) # autocorrelation check autocor(ral_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(ral_m) # assign to standard object for plotting ral_FINAL <- ral_m ``` ## plot ```{r} ral_plot <- gg_cities( d = ral, model = ral_FINAL, event = as.yearqtr("2016 Q1")) ral_plot ``` # Sacramento ```{r} # view city events protested %>% filter(city == "Sacramento") ``` ## Dazion Flenaugh ```{r} # extract data sac_df <- its("Sacramento", "2016 Q2") # analysis sac_df_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = sac_df %>% filter(qtr != "2016 Q2")) look(sac_df_m) # non-significant seasonal terms excluded sac_df_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + cos1 + cos2, method = "ML", data = sac_df %>% filter(qtr != "2016 Q2")) look(sac_df_m) # autocorrelation check autocor(sac_df_m) # adding ARMA terms sac_df_p1 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + cos1 + cos2, method = "ML", correlation = corARMA(p = 1, form = ~time | group), data = sac_df %>% filter(qtr != "2016 Q2")) sac_df_q1 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + cos1 + cos2, method = "ML", correlation = corARMA(q = 1, form = ~time | group), data = sac_df %>% filter(qtr != "2016 Q2")) sac_df_p1q1 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + cos1 + cos2, method = "ML", correlation = corARMA(p = 1, q = 1, form = ~time | group), data = sac_df %>% filter(qtr != "2016 Q2")) # compare model AICs AIC(sac_df_m, sac_df_p1, sac_df_q1, sac_df_p1q1) # best fit model outputs look(sac_df_p1q1) # assign to standard object for plotting sac_df_FINAL <- sac_df_p1q1 ``` ## Joseph Mann ```{r} # extract data sac_jm <- its("Sacramento", "2016 Q3") # analysis sac_jm_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = sac_jm %>% filter(qtr != "2016 Q3")) look(sac_jm_m) # non-significant seasonal terms excluded sac_jm_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + cos2, method = "ML", data = sac_jm %>% filter(qtr != "2016 Q3")) look(sac_jm_m) # autocorrelation check autocor(sac_jm_m) # adding ARMA terms sac_jm_p1 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + cos2, method = "ML", correlation = corARMA(p = 1, form = ~time | group), data = sac_jm %>% filter(qtr != "2016 Q3")) sac_jm_q1 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + cos2, method = "ML", correlation = corARMA(q = 1, form = ~time | group), data = sac_jm %>% filter(qtr != "2016 Q3")) sac_jm_p1q1 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + cos2, method = "ML", correlation = corARMA(p = 1, q = 1, form = ~time | group), data = sac_jm %>% filter(qtr != "2016 Q3")) # compare model AICs AIC(sac_jm_m, sac_jm_p1, sac_jm_q1, sac_jm_p1q1) # best fit model outputs look(sac_jm_p1q1) # assign to standard object for plotting sac_jm_FINAL <- sac_jm_p1q1 ``` ## compare events ```{r} AIC(sac_df_FINAL, sac_jm_FINAL) sac_df_FINAL ``` ## plot ```{r} sac_plot <- gg_cities( d = sac_df, model = sac_df_FINAL, event = as.yearqtr("2016 Q2")) sac_plot ``` # San Antonio ```{r} # view city events protested %>% filter(city == "San Antonio") ``` ## Antronie Scott ```{r} # extract data sa <- its("San Antonio", "2016 Q1") # analysis sa_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = sa %>% filter(qtr != "2016 Q1")) look(sa_m) # non-significant seasonal terms excluded sa_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1, method = "ML", data = sa %>% filter(qtr != "2016 Q1")) look(sa_m) # autocorrelation check autocor(sa_m) # adding ARMA terms sa_p4 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1, correlation = corARMA(p = 4, form = ~time | group), method = "ML", data = sa %>% filter(qtr != "2016 Q1")) sa_q4 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1, correlation = corARMA(q = 4, form = ~time | group), method = "ML", data = sa %>% filter(qtr != "2016 Q1")) sa_p4q4 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1, correlation = corARMA(p = 4, q = 4, form = ~time | group), method = "ML", data = sa %>% filter(qtr != "2016 Q1")) # compare model AICs AIC(sa_m, sa_p4, sa_q4, sa_p4q4) # best fit model outputs look(sa_p4q4) # assign to standard object for plotting sa_FINAL <- sa_p4q4 ``` ## plot ```{r} sa_plot <- gg_cities( d = sa, model = sa_FINAL, event = as.yearqtr("2016 Q1")) sa_plot ``` # San Diego ```{r} # view city events protested %>% filter(city %in% c("El Cajon", "San Diego")) ``` ## Alfred Olango (El Cajon) ```{r} # extract data sd_ao <- its("San Diego", "2016 Q3") # analysis sd_ao_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = sd_ao %>% filter(qtr != "2016 Q3")) look(sd_ao_m) # non-significant seasonal terms excluded sd_ao_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + cos2, method = "ML", data = sd_ao %>% filter(qtr != "2016 Q3")) look(sd_ao_m) # autocorrelation check autocor(sd_ao_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(sd_ao_m) # assign to standard object for plotting sd_ao_FINAL <- sd_ao_m ``` ## Lamontez Jones ```{r} # extract data sd_lj <- its("San Diego", "2015 Q4") # analysis sd_lj_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = sd_lj %>% filter(qtr != "2015 Q4")) look(sd_lj_m) # non-significant seasonal terms excluded sd_lj_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + cos2, method = "ML", data = sd_lj %>% filter(qtr != "2015 Q4")) look(sd_lj_m) # autocorrelation check autocor(sd_lj_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(sd_lj_m) # assign to standard object for plotting sd_lj_FINAL <- sd_lj_m ``` ## compare events ```{r} AIC(sd_ao_FINAL, sd_lj_FINAL) sd_lj_FINAL ``` ## plot ```{r} sd_plot <- gg_cities( d = sd_lj, model = sd_lj_FINAL, event = as.yearqtr("2015 Q4")) sd_plot ``` # San Francisco ```{r} # view city events protested %>% filter(city == "San Francisco") ``` ## Amilcar Perez-Lopez/Alice Brown ```{r} # extract data sf_aplab <- its("San Francisco", "2015 Q1") # analysis sf_aplab_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = sf_aplab %>% filter(qtr != "2015 Q1")) look(sf_aplab_m) # non-significant seasonal terms excluded sf_aplab_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group, method = "ML", data = sf_aplab %>% filter(qtr != "2015 Q1")) look(sf_aplab_m) # autocorrelation check autocor(sf_aplab_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(sf_aplab_m) # assign to standard object for plotting sf_aplab_FINAL <- sf_aplab_m ``` ## Mario Woods - death ```{r} # extract data sf_mw1 <- its("San Francisco", "2015 Q4") # analysis sf_mw1_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = sf_mw1 %>% filter(qtr != "2015 Q4")) look(sf_mw1_m) # non-significant seasonal terms excluded sf_mw1_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group, method = "ML", data = sf_mw1 %>% filter(qtr != "2015 Q4")) look(sf_mw1_m) # autocorrelation check autocor(sf_mw1_m) # adding ARMA terms sf_mw1_p3 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group, method = "ML", correlation = corARMA(p = 3, form = ~time | group), data = sf_mw1 %>% filter(qtr != "2015 Q4")) sf_mw1_q3 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group, method = "ML", correlation = corARMA(q = 3, form = ~time | group), data = sf_mw1 %>% filter(qtr != "2015 Q4")) sf_mw1_p3q3 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group, method = "ML", correlation = corARMA(p = 3, q = 3, form = ~time | group), data = sf_mw1 %>% filter(qtr != "2015 Q4")) # compare model AICs AIC(sf_mw1_m, sf_mw1_p3, sf_mw1_q3, sf_mw1_p3q3) # best fit model outputs look(sf_mw1_q3) # assign to standard object for plotting sf_mw1_FINAL <- sf_mw1_q3 ``` ## Mario Woods - protest timed to coincide with super bowl ```{r} # extract data sf_mw2 <- its("San Francisco", "2016 Q1") # analysis sf_mw2_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = sf_mw2 %>% filter(qtr != "2016 Q1")) look(sf_mw2_m) # non-significant seasonal terms excluded sf_mw2_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group, method = "ML", data = sf_mw2 %>% filter(qtr != "2016 Q1")) look(sf_mw2_m) # autocorrelation check autocor(sf_mw2_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(sf_mw2_m) # assign to standard object for plotting sf_mw2_FINAL <- sf_mw2_m ``` ## Jessica Nelson-Williams ```{r} # extract data sf_jnw <- its("San Francisco", "2016 Q2") # analysis sf_jnw_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = sf_jnw %>% filter(qtr != "2016 Q2")) look(sf_jnw_m) # non-significant seasonal terms excluded sf_jnw_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group, method = "ML", data = sf_jnw %>% filter(qtr != "2016 Q2")) look(sf_jnw_m) # autocorrelation check autocor(sf_jnw_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(sf_jnw_m) # assign to standard object for plotting sf_jnw_FINAL <- sf_jnw_m ``` ## compare events ```{r} AIC(sf_aplab_FINAL, sf_mw1_FINAL, sf_mw2_FINAL, sf_jnw_FINAL) sf_mw1_FINAL ``` ## plot ```{r} sf_plot <- gg_cities( d = sf_mw1, model = sf_mw1_FINAL, event = as.yearqtr("2015 Q4")) sf_plot ``` # Seattle ```{r} # view city events protested %>% filter(city == "Seattle") ``` ## Che Taylor ```{r} # extract data sea <- its("Seattle", "2016 Q1") # analysis sea_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = sea %>% filter(qtr != "2016 Q1")) look(sea_m) # non-significant seasonal terms excluded sea_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1, method = "ML", data = sea %>% filter(qtr != "2016 Q1")) look(sea_m) # autocorrelation check autocor(sea_m) # adding ARMA terms sea_p4 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1, method = "ML", correlation = corARMA(p = 4, form = ~time | group), data = sea %>% filter(qtr != "2016 Q1")) sea_q4 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1, method = "ML", correlation = corARMA(q = 4, form = ~time | group), data = sea %>% filter(qtr != "2016 Q1")) sea_p4q4 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1, method = "ML", correlation = corARMA(p = 4, q = 4, form = ~time | group), data = sea %>% filter(qtr != "2016 Q1")) # compare model AICs AIC(sea_m, sea_p4, sea_q4, sea_p4q4) # best fit model outputs look(sea_p4q4) # assign to standard object for plotting sea_FINAL <- sea_p4q4 ``` ## plot ```{r} sea_plot <- gg_cities( d = sea, model = sea_FINAL, event = as.yearqtr("2016 Q1")) sea_plot ``` # St Louis ```{r} # view city events protested %>% filter(city %in% c("Ferguson", "St. Louis")) ``` ## Michael Brown (Ferguson) - shooting ```{r} # extract data stl_mb1 <- its("St Louis", "2014 Q3") # analysis stl_mb1_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = stl_mb1 %>% filter(qtr != "2014 Q3")) look(stl_mb1_m) # non-significant seasonal terms excluded stl_mb1_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", data = stl_mb1 %>% filter(qtr != "2014 Q3")) look(stl_mb1_m) # autocorrelation check autocor(stl_mb1_m) # adding ARMA terms stl_mb1_p2 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", correlation = corARMA(p = 2, form = ~time | group), data = stl_mb1 %>% filter(qtr != "2014 Q3")) stl_mb1_q2 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", correlation = corARMA(q = 2, form = ~time | group), data = stl_mb1 %>% filter(qtr != "2014 Q3")) stl_mb1_p2q2 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", correlation = corARMA(p = 2, q = 2, form = ~time | group), data = stl_mb1 %>% filter(qtr != "2014 Q3")) # compare model AICs AIC(stl_mb1_m, stl_mb1_p2, stl_mb1_q2, stl_mb1_p2q2) # best fit model outputs look(stl_mb1_q2) # assign to standard object for plotting stl_mb1_FINAL <- stl_mb1_q2 ``` ## Michael Brown (Ferguson) - officers not indicted ```{r} # extract data stl_mb2 <- its("St Louis", "2014 Q4") # analysis stl_mb2_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = stl_mb2 %>% filter(qtr != "2014 Q4")) look(stl_mb2_m) # non-significant seasonal terms excluded stl_mb2_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", data = stl_mb2 %>% filter(qtr != "2014 Q4")) look(stl_mb2_m) # autocorrelation check autocor(stl_mb2_m) # adding ARMA terms stl_mb2_p2 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", correlation = corARMA(p = 2, form = ~time | group), data = stl_mb2 %>% filter(qtr != "2014 Q4")) stl_mb2_q2 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", correlation = corARMA(q = 2, form = ~time | group), data = stl_mb2 %>% filter(qtr != "2014 Q4")) stl_mb2_p2q2 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", correlation = corARMA(p = 2, q = 2, form = ~time | group), data = stl_mb2 %>% filter(qtr != "2014 Q4")) # compare model AICs AIC(stl_mb2_m, stl_mb2_p2, stl_mb2_q2, stl_mb2_p2q2) # best fit model outputs look(stl_mb2_q2) # assign to standard object for plotting stl_mb2_FINAL <- stl_mb2_q2 ``` ## Isaac Holmes ```{r} # extract data stl_ih <- its("St Louis", "2015 Q1") # analysis stl_ih_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = stl_ih %>% filter(qtr != "2015 Q1")) look(stl_ih_m) # non-significant seasonal terms excluded stl_ih_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", data = stl_ih %>% filter(qtr != "2015 Q1")) look(stl_ih_m) # autocorrelation check autocor(stl_ih_m) # adding ARMA terms stl_ih_p2 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", correlation = corARMA(p = 2, form = ~time | group), data = stl_ih %>% filter(qtr != "2015 Q1")) stl_ih_q2 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", correlation = corARMA(q = 2, form = ~time | group), data = stl_ih %>% filter(qtr != "2015 Q1")) stl_ih_p2q2 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", correlation = corARMA(p = 2, q = 2, form = ~time | group), data = stl_ih %>% filter(qtr != "2015 Q1")) # compare model AICs AIC(stl_ih_m, stl_ih_p2, stl_ih_q2, stl_ih_p2q2) # best fit model outputs look(stl_ih_p2q2) # assign to standard object for plotting stl_ih_FINAL <- stl_ih_p2q2 ``` ## Crayton West ```{r} # extract data stl_cw <- its("St Louis", "2016 Q1") # analysis stl_cw_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = stl_cw %>% filter(qtr != "2016 Q3")) look(stl_cw_m) # non-significant seasonal terms excluded stl_cw_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", data = stl_cw %>% filter(qtr != "2016 Q3")) look(stl_cw_m) # autocorrelation check autocor(stl_cw_m) # adding ARMA terms stl_cw_p2 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", correlation = corARMA(p = 2, form = ~time | group), data = stl_cw %>% filter(qtr != "2016 Q3")) stl_cw_q2 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", correlation = corARMA(q = 2, form = ~time | group), data = stl_cw %>% filter(qtr != "2016 Q3")) stl_cw_p2q2 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + cos1, method = "ML", correlation = corARMA(p = 2, q = 2, form = ~time | group), data = stl_cw %>% filter(qtr != "2016 Q3")) # compare model AICs AIC(stl_cw_m, stl_cw_p2, stl_cw_q2, stl_cw_p2q2) # best fit model outputs look(stl_cw_p2q2) # assign to standard object for plotting stl_cw_FINAL <- stl_cw_p2q2 ``` ## compare events ```{r} AIC(stl_mb1_FINAL, stl_mb2_FINAL, stl_ih_FINAL, stl_cw_FINAL) stl_mb2_FINAL ``` ## plot ```{r} stl_plot <- gg_cities( d = stl_mb2, model = stl_mb2_FINAL, event = as.yearqtr("2014 Q4")) stl_plot ``` # Stockton ```{r} # view city events protested %>% filter(city == "Stockton") ``` ## Matautu Nuu ```{r} # extract data sto <- its("Stockton", "2015 Q1") # analysis sto_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = sto %>% filter(qtr != "2015 Q1")) look(sto_m) # non-significant seasonal terms excluded sto_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1, method = "ML", data = sto %>% filter(qtr != "2015 Q1")) look(sto_m) # autocorrelation check autocor(sto_m) # adding ARMA terms sto_p2 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1, method = "ML", correlation = corARMA(p = 2, form = ~time | group), data = sto %>% filter(qtr != "2015 Q1")) sto_q2 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1, method = "ML", correlation = corARMA(q = 2, form = ~time | group), data = sto %>% filter(qtr != "2015 Q1")) sto_p2q2 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1, method = "ML", correlation = corARMA(p = 2, q = 2, form = ~time | group), data = sto %>% filter(qtr != "2015 Q1")) # compare model AICs AIC(sto_m, sto_p2, sto_q2, sto_p2q2) # best fit model outputs look(sto_p2q2) # assign to standard object for plotting sto_FINAL <- sto_p2q2 ``` ## plot ```{r} sto_plot <- gg_cities( d = sto, model = sto_FINAL, event = as.yearqtr("2015 Q1")) sto_plot ``` # Tampa ```{r} # view city events protested %>% filter(city == "Tampa") ``` ## Levonia Riggins ```{r} # extract data tam <- its("Tampa", "2016 Q3") # analysis tam_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = tam %>% filter(qtr != "2016 Q3")) look(tam_m) # non-significant seasonal terms excluded tam_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group, method = "ML", data = tam %>% filter(qtr != "2016 Q3")) look(tam_m) # autocorrelation check autocor(tam_m) # adding ARMA terms tam_p1 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group, method = "ML", correlation = corARMA(p = 1, form = ~ time | group), data = tam %>% filter(qtr != "2016 Q3")) tam_q1 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group, method = "ML", correlation = corARMA(q = 1, form = ~ time | group), data = tam %>% filter(qtr != "2016 Q3")) tam_p1q1 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group, method = "ML", correlation = corARMA(p = 1, q = 1, form = ~ time | group), data = tam %>% filter(qtr != "2016 Q3")) # compare model AICs AIC(tam_m, tam_p1, tam_q1, tam_p1q1) # best fit model outputs look(tam_p1q1) # assign to standard object for plotting tam_FINAL <- tam_p1q1 ``` ## plot ```{r} tam_plot <- gg_cities( d = tam, model = tam_FINAL, event = as.yearqtr("2016 Q3")) tam_plot ``` # Tulsa ```{r} # view city events protested %>% filter(city == "Tulsa") ``` ## Eric Harris ```{r} # extract data tul_eh <- its("Tulsa", "2015 Q2") # analysis tul_eh_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = tul_eh %>% filter(qtr != "2015 Q2")) look(tul_eh_m) # non-significant seasonal terms excluded tul_eh_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + cos1, method = "ML", data = tul_eh %>% filter(qtr != "2015 Q2")) look(tul_eh_m) # autocorrelation check autocor(tul_eh_m) # adding ARMA terms tul_eh_p4 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + cos1, method = "ML", correlation = corARMA(p = 4, form = ~time | group), data = tul_eh %>% filter(qtr != "2015 Q2")) tul_eh_q4 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + cos1, method = "ML", correlation = corARMA(q = 4, form = ~time | group), data = tul_eh %>% filter(qtr != "2015 Q2")) tul_eh_p4q4 <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + cos1, method = "ML", correlation = corARMA(p = 4, q = 4, form = ~time | group), data = tul_eh %>% filter(qtr != "2015 Q2")) # compare model AICs AIC(tul_eh_m, tul_eh_p4, tul_eh_q4, tul_eh_p4q4) # best fit model outputs look(tul_eh_p4q4) # assign to standard object for plotting tul_eh_FINAL <- tul_eh_p4 ``` ## Terence Crutcher ```{r} # extract data tul_tc <- its("Tulsa", "2016 Q3") # analysis tul_tc_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = tul_tc %>% filter(qtr != "2016 Q3")) look(tul_tc_m) # non-significant seasonal terms excluded tul_tc_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1, method = "ML", data = tul_tc %>% filter(qtr != "2016 Q3")) look(tul_tc_m) # autocorrelation check autocor(tul_tc_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(tul_tc_m) # assign to standard object for plotting tul_tc_FINAL <- tul_tc_m ``` ## compare events ```{r} AIC(tul_eh_FINAL, tul_tc_FINAL) tul_eh_FINAL ``` ## plot ```{r} tul_plot <- gg_cities( d = tul_eh, model = tul_eh_FINAL, event = as.yearqtr("2015 Q2")) tul_plot ``` # Washington DC ```{r} # view city events protested %>% filter(city == "Washington") ``` ## Terrence Sterling - death ```{r} # extract data was_ts1 <- its("Washington DC", "2016 Q3") # analysis was_ts1_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = was_ts1 %>% filter(qtr != "2016 Q3")) look(was_ts1_m) # non-significant seasonal terms excluded was_ts1_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1, method = "ML", data = was_ts1 %>% filter(qtr != "2016 Q3")) look(was_ts1_m) # autocorrelation check autocor(was_ts1_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(was_ts1_m) # assign to standard object for plotting was_ts1_FINAL <- was_ts1_m ``` ## Terrence Sterling - initial protests ```{r} # extract data was_ts2 <- its("Washington DC", "2016 Q4") # analysis was_ts2_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = was_ts2 %>% filter(qtr != "2016 Q4")) look(was_ts2_m) # non-significant seasonal terms excluded was_ts2_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + cos1, method = "ML", data = was_ts2 %>% filter(qtr != "2016 Q4")) look(was_ts2_m) # autocorrelation check autocor(was_ts2_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(was_ts2_m) # assign to standard object for plotting was_ts2_FINAL <- was_ts2_m ``` ## compare events ```{r} AIC(was_ts1_FINAL, was_ts2_FINAL) was_ts1_FINAL ``` ## plot ```{r} was_plot <- gg_cities( d = was_ts1, model = was_ts1_FINAL, event = as.yearqtr("2016 Q3")) was_plot ``` # Wichita ```{r} # view city events protested %>% filter(city == "Wichita") ``` ## John Paul Quintero ```{r} # extract data wic <- its("Wichita", "2015 Q1") # analysis wic_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group + sin1 + sin2 + cos1 + cos2, method = "ML", data = wic %>% filter(qtr != "2015 Q1")) look(wic_m) # non-significant seasonal terms excluded wic_m <- gls( log(rate) ~ level + trend + time + level_assault + trend_assault + time_assault + group, method = "ML", data = wic %>% filter(qtr != "2015 Q1")) look(wic_m) # autocorrelation check autocor(wic_m) # adding ARMA terms ## no residual autocorrelation # compare model AICs # best fit model outputs look(wic_m) # assign to standard object for plotting wic_FINAL <- wic_m ``` ## plot ```{r} wic_plot <- gg_cities( d = wic, model = wic_FINAL, event = as.yearqtr("2015 Q1")) wic_plot ``` # combining plotted results ```{r arrange plots} plots <- ggarrange( alb_plot, asa_plot, atl_plot, aur_plot, aus_plot, bak_plot, bal_plot, br_plot, cha_sc_plot, cha_nc_plot, chi_plot, cin_plot, cle_plot, col_plot, dfwa_plot, den_plot, det_plot, dur_plot, fre_plot, hou_plot, ind_plot, lb_plot, la_plot, lou_plot, mem_plot, mil_plot, msp_plot, no_plot, ny_plot, oak_plot, phi_plot, pho_plot, ral_plot, sac_plot, sa_plot, sd_plot, sf_plot, sea_plot, stl_plot, sto_plot, tam_plot, tul_plot, was_plot, wic_plot, ncol = 4, nrow = 11, align = "hv", common.legend = TRUE, legend = "bottom") ggsave("s1 homicide & assault models.png", plots, height = 16, width = 12) ```