Analysis of how religiousity contributes to firearm deaths.

religiousity_by_firearm_deaths.R

library(tidyverse)
read_csv("~/Downloads/Religion by Firearm Deaths Data - Firearm Deaths by State.csv",
         skip = 4) %>%
  janitor::clean_names() -> firearm_deaths

read_csv("~/Downloads/Religion by Firearm Deaths Data - Religiousity by State.csv",
         skip = 3) %>%
  janitor::clean_names() %>%
  rename(state = region) %>%
  mutate(religious = 1 - irreligion_percent/100) %>%
  select(state, reliousity = religious) -> religiousity_by_state

firearm_deaths %>%
  select(state, firearm_death_rate_per_100k_per_year = rate) %>%
  left_join(religiousity_by_state, by = "state") %>%
  filter_all(Negate(is.na)) -> d

rlm(firearm_death_rate_per_100k_per_year ~ reliousity, data = d) -> fit
summary(fit)

library(MASS); select <- dplyr::select
d %>%
  ggplot(aes(reliousity, firearm_death_rate_per_100k_per_year)) +
  geom_point(size = 2) +
  # geom_smooth(method = "lm", size = 2) +
  geom_smooth(method = "gam", formula = y ~ s(x, bs = "ps")) +
  ggtitle("Religiosity causes firearm deaths in part (by US state)") +
  labs(y = "Firearm deaths per 100,000 state residents",
       x = "Population identifying as religious") +
  scale_y_continuous(breaks = seq(0, 50, 10)) +
  scale_x_continuous(breaks = seq(0, 1, 0.05),
                     labels = scales::percent) +
  theme_bw(15) +
  expand_limits(y = 0) +
  theme(panel.grid.minor = element_blank(),
        panel.grid.major = element_line(size = 0.1, color = "black"),
        axis.text = element_text(color = "black"))