Examining Mass Shootings and Fatalities

Ernesto Gomez

October 30, 2017

---

Introduction

Recent events involving mass shooters gunning down large groups of people in the United States has catalyzed a national conversation around guns, gun ownership, background checks, and mental health. This analysis is concerned with the question around access to guns and whether access to a larger arsenal on the part of a mass shooter is likelier to result in more fatalities. A multilevel method will be utilized to further analyze the relationship between fatalities and weaponry size and possible variations amongst all 50 states, who have varying rates of gun ownership. In the end, mental illness will be introduced to test for its influence in fatalities. The analysis will begin with an ecological analysis, complete pooling model, a no pool modeling, and conclude with multilevel modeling to ultimately determine the best fit and contributing factors to fatalities.

---

Method

To conduct this analysis, the primary data used was extracted from the Stanford Mass Shootings in America (MSA). Additionally, the percentage of gun ownership across the United States, taken from a heatmap based on General Social Survey data and stratified on the state-level, is also utilized and merged into the dataset.

Variables contained within the dataframe are defined and coded as follows:

• totalvictims - This variable will be the dependent variable, it is a measurement of the total number of fatalities for every incidence of a mass shooting.
• mentalhistory - This variable measures whether a mass shooter had a reported or confirmed history of mental illness, “Yes” for a history of mental illness and “No” for no history.
• mentalill - A recoded 1,0 variable of history of mental illness, 1 for “Yes” and 0 for “No”.
• totalguns - This variable is a measurement of the total number of guns used during a particular shooting.
• statefull - The state in which a mass shooting incidence took place.
• PercHavGun - The percentage of gun ownership from a particular state (derived from the General Social Survey’s Trends in Gun Ownership in the United States) matched by state.

library(nlme)
library(tidyverse)
library(dplyr)
library(sjmisc)
library(magrittr)
library(lmerTest)
library(texreg)
library(visreg)
library(ggplot2)
masshoot <- data.frame(read_csv("/Users/ernesto/Documents/Advanced Analytics/Data/MassShootings.csv")) %>%
  rename("mentalhistory" = History.of.Mental.Illness...General,
    "totalvictims" = Total.Number.of.Victims,
    "state" = State,
     "totalguns" = Total.Number.of.Guns) %>%
  select(state, totalvictims, mentalhistory, totalguns) %>%
  filter(!is.na(totalvictims),
         !is.na(state),
         !is.na(mentalhistory),
         !is.na(totalguns), 
        !totalvictims %in% c("Unknown"),
        !state %in% c("Unknown"),
        !mentalhistory %in% c("Unknown"),
        !totalguns %in% c("Unknown")) %>%
          mutate(mentalill=as.factor(ifelse(mentalhistory=="Yes", 1, 0)),
              totalguns=as.numeric(totalguns))
              #newtotalguns=as.factor(ifelse(totalguns<4, totalguns, 
                                           #ifelse(totalguns>=4, "4+","na"))))
statefull = c("Alaska","Alabama","Arkansas", 
                              "Arizona","California","Colorado",
                              "Connecticut","District of Columbia","Delaware",
                              "Florida" , "Georgia" ,"Hawaii",
                              "Iowa","Idaho","Illinois","Indiana","Kansas",
                              "Kentucky","Louisiana","Massachusetts",
                              "Maryland","Maine" ,"Michigan","Minnesota",
                              "Missouri","Mississippi","Montana",
                              "North Carolina","North Dakota","Nebraska",
                              "New Hampshire" , "New Jersey" ,  "New Mexico",
                              "Nevada","New York","Ohio","Oklahoma" ,"Oregon",
                              "Pennsylvania" , "Rhode Island",
                              "South Carolina", "South Dakota" ,"Tennessee" , 
                              "Texas" , "Utah", "Virginia","Vermont" ,"Washington",
                              "Wisconsin", "West Virginia" , "Wyoming")
PercHavGuns= c(57.8, 51.7, 55.3,
               31.1, 21.3, 34.7,
               16.7, 3.6, 25.5,
               24.5, 40.3, 6.7,
               42.9, 55.3, 20.2, 39.1, 42.1,
               47.7, 44.1, 12.6,
               21.3, 40.5, 38.4, 41.7,
               41.7, 55.3, 57.7,
               41.3, 50.7, 38.6,
               30.0, 12.3, 34.8,
               33.8, 18.0, 32.4, 42.9, 39.8,
               34.7, 12.8,
               42.3, 56.6, 43.9,
               35.9, 43.9, 35.1, 42.0, 33.1,
               44.4, 55.4, 59.7)
stateguns <- data.frame(statefull, PercHavGuns)
mshoot <- merge(stateguns, masshoot, all.x=FALSE, all.y=TRUE, by.x="statefull", by.y="state") %>%
  mutate(statefull=as.factor(statefull))
print(mshoot)

---

Ecological Analysis

The ecological regression seeks to understand whether, at the state level, does the number of guns owned relate to the number of mass shooting incidences in that state? What is found is that simply because a particular state has a higher rate of gun ownership this does not necessarily result in higher incidences of mass shootings. This is somewhat interesting because one might expect a state with more guns to mean easier access, greater ubiquity of access to guns, individuals who own more than one gun, and thus potentially a higher likelihood of more shootings and subsequently deaths.

shootd <- mshoot %>%
  group_by(statefull) %>%
  summarise(gunown=mean(PercHavGuns), totalvictims=sum(totalvictims), incidence=n())
summary(lm(incidence ~ gunown, data = shootd))

Call:
lm(formula = incidence ~ gunown, data = shootd)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.8074 -2.0221 -0.8676  1.1989 20.5414 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  7.08169    1.86923   3.789 0.000488 ***
gunown      -0.07620    0.04687  -1.626 0.111657    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.264 on 41 degrees of freedom
Multiple R-squared:  0.06056,   Adjusted R-squared:  0.03765 
F-statistic: 2.643 on 1 and 41 DF,  p-value: 0.1117

---

Fatalities and Number of Guns Used - A Complete Pooling Model

This individual analysis/complete pooling model examines the number of guns used during an incidence and whether this influences the number of fatalities. The model’s results indicate a highly statistically significant increase in fatalities the more guns used by a shooter during an incident.

cpooling <- lm(totalvictims ~ totalguns, data=mshoot)
summary(cpooling)

Call:
lm(formula = totalvictims ~ totalguns, data = mshoot)

Residuals:
    Min      1Q  Median      3Q     Max 
-19.701  -3.323  -2.206   0.677  54.217 

Coefficients:
            Estimate Std. Error t value           Pr(>|t|)    
(Intercept)   3.1700     0.8145   3.892           0.000127 ***
totalguns     3.1531     0.3990   7.903 0.0000000000000775 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.706 on 260 degrees of freedom
Multiple R-squared:  0.1937,    Adjusted R-squared:  0.1906 
F-statistic: 62.46 on 1 and 260 DF,  p-value: 0.00000000000007745

---

Fatalities and Number of Guns - A No Pooling Model

---

Intercept

When considering the intercept model, which is an indication of fatalities, we see that most intercepts range between 0-18 when there is no consideration of number of guns used.

dcoef <- mshoot %>%
  group_by(statefull) %>%
  do(mod=lm(totalvictims ~ totalguns, data=.))
coef <- dcoef %>% do(data.frame(intc=coef(.$mod)[1]))
ggplot(coef, aes(x=intc)) + geom_histogram(fill="black")

---

Slope

Our slope models show a consistently positive impact, suggesting that the more guns used during a particular shooting resulted in more fatalities.

dcoef1 <- mshoot %>%
  group_by(statefull) %>%
  do(mod=lm(totalvictims ~ totalguns, data=.))
coef1 <- dcoef1 %>% do(data.frame(totalgunc = coef(.$mod)[2]))
ggplot (coef1, aes(x=totalgunc)) + geom_histogram(fill="darkred")

---

Multilevel Modeling

---

Fatalities & Number of Guns by State - A Random Intercept Model

m1_lme <- lme(totalvictims~totalguns, data=mshoot, random = ~1|statefull, method="ML")
summary(m1_lme)

Linear mixed-effects model fit by maximum likelihood
 Data: mshoot 
       AIC      BIC    logLik
  1816.377 1830.651 -904.1887

Random effects:
 Formula: ~1 | statefull
        (Intercept) Residual
StdDev:    1.854707  7.44623

Fixed effects: totalvictims ~ totalguns 
               Value Std.Error  DF  t-value p-value
(Intercept) 3.152600 0.8656005 215 3.642096  0.0003
totalguns   3.076404 0.3953991 215 7.780504  0.0000
 Correlation: 
          (Intr)
totalguns -0.747

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-2.51797073 -0.38017819 -0.19661757  0.09522077  6.91271424 

Number of Observations: 262
Number of Groups: 46 

---

Fatalities and Number of Guns by State - A Random Slope Model

m2_lme <- lme(totalvictims~totalguns, data=mshoot, random = ~totalguns|statefull, method="ML")
summary(m2_lme)

Linear mixed-effects model fit by maximum likelihood
 Data: mshoot 
       AIC      BIC    logLik
  1779.533 1800.943 -883.7663

Random effects:
 Formula: ~totalguns | statefull
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev   Corr  
(Intercept) 3.066593 (Intr)
totalguns   3.203568 -0.987
Residual    6.566262       

Fixed effects: totalvictims ~ totalguns 
               Value Std.Error  DF  t-value p-value
(Intercept) 2.708586 1.0514729 215 2.575992  0.0107
totalguns   3.447483 0.7935152 215 4.344570  0.0000
 Correlation: 
          (Intr)
totalguns -0.892

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-3.5028251 -0.4510911 -0.2105835  0.1330508  5.2161427 

Number of Observations: 262
Number of Groups: 46 

Model(s) Summary

When accounting for the variations/differences between states we continue to see a statistically significant relationship between number of guns used and fatalities, which still increase as the number of guns used similarly increases.

---

Introducing Covariate - The Effect of Mental Illness

The model maintains its significance in regards to number of guns used, and the model also reveal statistical significance and an increase in fatalities when considering mental illness amongst convicted mass shooters.

m3_lme <- lme(totalvictims ~ totalguns+mentalill, data=mshoot, random = ~totalguns|statefull, method="ML")
summary(m3_lme)

Linear mixed-effects model fit by maximum likelihood
 Data: mshoot 
       AIC      BIC   logLik
  1765.498 1790.476 -875.749

Random effects:
 Formula: ~totalguns | statefull
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev   Corr  
(Intercept) 2.721248 (Intr)
totalguns   2.925215 -0.987
Residual    6.393374       

Fixed effects: totalvictims ~ totalguns + mentalill 
               Value Std.Error  DF  t-value p-value
(Intercept) 1.891857 1.0171853 214 1.859894  0.0643
totalguns   3.164054 0.7467198 214 4.237271  0.0000
mentalill1  3.541143 0.8743448 214 4.050054  0.0001
 Correlation: 
           (Intr) ttlgns
totalguns  -0.845       
mentalill1 -0.208 -0.084

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-3.2790555 -0.4322534 -0.1487146  0.2175075  5.2065406 

Number of Observations: 262
Number of Groups: 46 

Mental Illness and Access to Many Guns by State - An Interaction Model

When interacting number of guns used and mental illness – do shooters with a history of mental history kill more with more guns – the model is weakened and loses significance. There may be no interactive effect.

m4_lme <- lme(totalvictims ~ totalguns*mentalill, data=mshoot, random = ~totalguns|statefull, method="ML")
summary(m4_lme)

Linear mixed-effects model fit by maximum likelihood
 Data: mshoot 
       AIC      BIC    logLik
  1765.536 1794.082 -874.7678

Random effects:
 Formula: ~totalguns | statefull
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev   Corr  
(Intercept) 2.684024 (Intr)
totalguns   2.770195 -0.976
Residual    6.390776       

Fixed effects: totalvictims ~ totalguns * mentalill 
                        Value Std.Error  DF   t-value p-value
(Intercept)          2.748995 1.1913864 213 2.3073917  0.0220
totalguns            2.540068 0.8506662 213 2.9859753  0.0032
mentalill1           1.516140 1.6772752 213 0.9039303  0.3671
totalguns:mentalill1 1.342273 0.9392529 213 1.4290860  0.1544
 Correlation: 
                     (Intr) ttlgns mntll1
totalguns            -0.886              
mentalill1           -0.543  0.407       
totalguns:mentalill1  0.528 -0.523 -0.853

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-3.0193090 -0.4341919 -0.1737580  0.2379613  5.0829045 

Number of Observations: 262
Number of Groups: 46 

---

Model Selection

Our third multilevel model proves to be the best fit when models are grouped and viewed for the lowest AIC value.

AIC(cpooling, m1_lme, m2_lme, m3_lme, m4_lme)

Best Fit Models

htmlreg(list(m2_lme, m3_lme, m4_lme))

Statistical models

	Model 1	Model 2	Model 3
(Intercept)	2.71*	1.89	2.75*
	(1.05)	(1.02)	(1.19)
totalguns	3.45***	3.16***	2.54**
	(0.79)	(0.75)	(0.85)
mentalill1		3.54***	1.52
		(0.87)	(1.68)
totalguns:mentalill1			1.34
			(0.94)
AIC	1779.53	1765.50	1765.54
BIC	1800.94	1790.48	1794.08
Log Likelihood	-883.77	-875.75	-874.77
Num. obs.	262	262	262
Num. groups	46	46	46
p < 0.001, p < 0.01, p < 0.05

---

In Conclusion

A graphical representation of the best fit model indicates that more guns used means more fatalities. Additionally, if a shooter is reported to have a history of mental illness they are also likelier to kill more than those with no prior history. The interaction between the two showed no significant effect.

visreg(m3_lme)

Note

It was my intention to also include the covariate in the pooling/no pooling models but encountered an error last minute when adjusting the lm to account for the covariate in the slopes. I tried to understand how to fix it but couldn’t. My mistake.