Homework 6 - Titanic, Mass Shooters & Zelig

Ernesto Gomez

October 16, 2017

---

Part 1 of this assignment will first implement the Zelig 5 syntax on the computations presented in slides 5.2-5.18 from Lecture 5. This will use the titanic dataset.

Part 2 will then utilize data from The Stanford Mass Shootings in America (MSA), a non-exhaustive dataset of mass shootings from 1966-2016 compiled by The Stanford Geospatial Center. This collection of data includes variables such as mass shooters’ race, name, gender, cause of death, gun type, number of fatalities, and more.

For the purposes of this analysis, race, gender, and a shooter’s fate at the scene of the crime will be used to perform a logistic regression analyzed under the Zelig package to analyze whether a mass shooter’s race or gender influences their fate at the scene of their crime.

In other words, is one gender or one race more likely to die at the scene of their crime or be apprehended by law enforcement?

---

library(tidyverse)
library(dplyr)
library(sjmisc)
library(tidyr)
library(ggplot2)
library(pander)
library(radiant.data)
library(Zelig)
library(texreg)
library(visreg)

---

Part 1: Titanic and Zelig 5

titanic2 <- titanic %>% 
  mutate(survived1 = as.integer(survived)) %>%
  mutate (age = as.integer(age)) %>%
  mutate(survival = sjmisc::rec(survived1, rec = "2=0; 1=1")) %>%
  select(survived, survival, age, sex, pclass, fare)
head(titanic2)

---

Generating Zelig Model (Titanic)

z5tit<- zlogit$new()
z5tit$zelig(survival ~ age + sex*pclass + fare, data = titanic2)
summary(z5tit)

Model: 

Call:
z5tit$zelig(formula = survival ~ age + sex * pclass + fare, data = titanic2)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.0628  -0.6636  -0.4941   0.4337   2.4940  

Coefficients:
                    Estimate Std. Error z value            Pr(>|z|)
(Intercept)        4.8959649  0.6128145   7.989 0.00000000000000136
age               -0.0386781  0.0067926  -5.694 0.00000001239977237
sexmale           -3.9001038  0.5015680  -7.776 0.00000000000000750
pclass2nd         -1.5922712  0.6024689  -2.643             0.00822
pclass3rd         -4.1381922  0.5601346  -7.388 0.00000000000014922
fare              -0.0009074  0.0020510  -0.442             0.65820
sexmale:pclass2nd -0.0603255  0.6373047  -0.095             0.92459
sexmale:pclass3rd  2.5016703  0.5479814   4.565 0.00000498908018340

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1409.99  on 1042  degrees of freedom
Residual deviance:  931.42  on 1035  degrees of freedom
AIC: 947.42

Number of Fisher Scoring iterations: 5

Next step: Use 'setx' method

---

Age Effect With Zelig 5 (Titanic)

z5tit$setrange(age = min(titanic2$age):max(titanic2$age))
z5tit$sim()
ci.plot(z5tit)

---

Fare Effect With Zelig 5 (Titanic)

z5tit$setrange(fare = min(titanic2$fare):max(titanic2$fare))
z5tit$sim()
ci.plot(z5tit)

---

Gender Difference with Zelig 5 (Titanic)

z5sex <- zlogit$new()
z5sex$zelig(survival ~ age + sex*pclass + fare, data = titanic2)
  z5sex$setx(sex = "male") 
  z5sex$setx1(sex = "female")
  z5sex$sim()
  summary(z5sex)

 sim x :
 -----
ev
          mean         sd       50%      2.5%     97.5%
[1,] 0.1404976 0.01878571 0.1400461 0.1059377 0.1829745
pv
         0     1
[1,] 0.858 0.142

 sim x1 :
 -----
ev
          mean         sd       50%      2.5%     97.5%
[1,] 0.3968055 0.04144396 0.3971523 0.3197929 0.4782572
pv
         0     1
[1,] 0.619 0.381
fd
         mean         sd       50%      2.5%     97.5%
[1,] 0.256308 0.04312428 0.2574681 0.1727245 0.3425787

  plot(z5sex)

---

Class Variations in Gender Differences with Zelig 5 (Titanic)

Gender and First Class

z5.1s <- zlogit$new()
z5.1s$zelig(survival ~ age + sex*pclass + fare, data = titanic2)
z5.1s$setx(sex = "male", pclass = "1st") 
z5.1s$setx1(sex = "female", pclass = "1st")
z5.1s$sim()
plot(z5.1s)

---

Gender and Second Class

z5.2s <- zlogit$new()
z5.2s$zelig(survival ~ age + sex*pclass + fare, data = titanic2)
z5.2s$setx(sex = "male", pclass = "2nd") 
z5.2s$setx1(sex = "female", pclass = "2nd")
z5.2s$sim()
plot(z5.2s)

---

Gender and Third Class

z5.3s <- zlogit$new()
z5.3s$zelig(survival ~ age + sex*pclass + fare, data = titanic2)
z5.3s$setx(sex = "male", pclass = "3rd") 
z5.3s$setx1(sex = "female", pclass = "3rd")
z5.3s$sim()
plot(z5.3s)

---

Grouping Gender Differences Between Classes (Titanic)

d1 <- z5.1s$get_qi(xvalue="x1", qi="fd")
d2 <- z5.2s$get_qi(xvalue="x1", qi="fd")
d3 <- z5.3s$get_qi(xvalue="x1", qi="fd")
  
dfd <- as.data.frame(cbind(d1, d2, d3))
head(dfd)

---

tidd <- dfd %>% 
  gather(class, simv, 1:3)
head(tidd)

---

tidd %>% 
  group_by(class) %>% 
  summarise(mean = mean(simv), sd = sd(simv))

---

Plotting Gender Differences Between Classes (Titanic)

ggplot(tidd, aes(simv)) + geom_histogram() + facet_grid(~class)

---

Part 2 - Mass Shooters Data & Zelig 4

Variables & Data Management

This analysis will use the following variables:

• gender - Measures whether a mass shooter is a Male or Female
• race - Records a mass shooter’s race (white, black, asian, two or more races)
• shooterfate - Records the fate of a shooter at the scene of their crime (killed, comitted sucide, or apprehended “i.e not applicable”)

The new variables, shooterwhite, was generated to represent mass shooters who are white and non-white, while the variable shooterdied represents shooters who died at the scene, either by suicide or by law enforcement, and those who did not die/were apprehended.

masshoot <- data.frame(read_csv("/Users/ernesto/Documents/Advanced Analytics/Data/MassShootings.csv")) %>%
  rename("shooterfate" = Shooter.s.Cause.of.Death,
    "shootergender" = Shooter.Sex,
     "shooterrace" = Shooter.Race) %>%
  select(shootergender, shooterrace, shooterfate) %>%
  filter(!is.na(shooterrace),
         !is.na(shootergender),
         !is.na(shooterfate), 
        !shooterfate %in% c("Unknown"),
        !shooterrace %in% c("Unknown"),
        !shootergender %in% c("Unknown", "Male/Female"),
        shootergender %in% c("Male", "Female")) %>%
      mutate(shootergender=as.factor(shootergender),
             shooterwhite=as.factor(ifelse(shooterrace=="White American or European American", 1, 0)),
             shooterfate = tolower(shooterfate),
             shooterdied=as.factor(ifelse(shooterfate %in% c("killed","suicide"), 1, 0)))
             
print(masshoot)

---

Overview of Data

Mass shooters are overwhelmingly white compared to other races.

ggplot(masshoot, aes(shooterrace) ) +
  geom_bar() + coord_flip()

---

Mass shooters also tend to be men.

ggplot(masshoot, aes(shootergender) ) +
  geom_bar() + coord_flip()

---

The majority of mass shooters are either apprehended (not applicable), however many of them commit suicide.

ggplot(masshoot, aes(shooterfate) ) +
  geom_bar() + coord_flip()

---

Generating Zelig Model

Men and whites more likely to die at the scene of their crime.

shootzelig <- zelig(shooterdied ~ shooterwhite + shootergender, model = "logit", data = masshoot, cite = F)
summary(shootzelig)

Model: 

Call:
z5$zelig(formula = shooterdied ~ shooterwhite + shootergender, 
    data = masshoot)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.2667  -0.9643  -0.9643   1.0907   1.4067  

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)
(Intercept)        -1.0021     0.9454  -1.060  0.28915
shooterwhite1       0.7319     0.2549   2.871  0.00409
shootergenderMale   0.4776     0.9332   0.512  0.60883

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 357.66  on 258  degrees of freedom
Residual deviance: 349.18  on 256  degrees of freedom
AIC: 355.18

Number of Fisher Scoring iterations: 4

Next step: Use 'setx' method

---

Generating Simulations - Racial Variations in Gender Difference

The first simulation finds that mass shooters who are white men are more likely to die at the scene of their crimes than white women who are mass shooters, with a first difference of -.11

a1x <-  setx(shootzelig, shootergender = "Male", shooterwhite=1)
a1x1 <- setx(shootzelig, shootergender= "Female", shooterwhite=1)
a1s <- sim(shootzelig, x = a1x, x1= a1x1)
summary(a1s)

 sim x :
 -----
ev
          mean         sd       50%      2.5%     97.5%
[1,] 0.5512355 0.04417524 0.5502608 0.4638718 0.6339924
pv
         0     1
[1,] 0.459 0.541

 sim x1 :
 -----
ev
          mean        sd       50%      2.5%     97.5%
[1,] 0.4399224 0.1939713 0.4302544 0.1083263 0.8134935
pv
         0     1
[1,] 0.559 0.441
fd
          mean       sd        50%       2.5%     97.5%
[1,] -0.111313 0.200032 -0.1180531 -0.4606997 0.2731931

plot(a1s)

---

The second simulation shows that mass shooters who are not white but male are also more likely to die at the scene of their crimes in comparison to non-white female shooters, with a first difference of -.07

Evidently, mass shooters who are men are more likely to die at the scene than women, whether they are white or non-white.

a2x <- setx(shootzelig, shootergender="Male", shooterwhite=0)
a2x1 <- setx(shootzelig, shootergender= "Female", shooterwhite=0)
a2s <- sim(shootzelig, x=a2x, x1=a2x1)
summary(a2s)

 sim x :
 -----
ev
          mean         sd       50%      2.5%     97.5%
[1,] 0.3710795 0.04167486 0.3684918 0.2903538 0.4608056
pv
         0     1
[1,] 0.639 0.361

 sim x1 :
 -----
ev
          mean        sd       50%      2.5%     97.5%
[1,] 0.3068707 0.1769439 0.2712629 0.0600497 0.7181584
pv
         0     1
[1,] 0.699 0.301
fd
            mean        sd         50%       2.5%     97.5%
[1,] -0.06420876 0.1751594 -0.09017618 -0.3255986 0.3404416

plot(a2s)

Combining Simulations Into Dataframe

fd1 <- a1s$get_qi(xvalue="x1", qi="fd")
fd2 <- a2s$get_qi(xvalue="x1", qi="fd")
newfd <- as.data.frame(cbind(fd1, fd2))
colnames(newfd)[1] <- "whiteM/F"
colnames(newfd)[2] <- "nonwhiteM/F"
head(newfd)

Group Summary Stats

Among those who are white (V1), male and female mass shooters are less likely to survive (killed or by suicide), with a gender difference of .11, compared to non-white male and female mass shooters (V2).

fdshooter%>%
  group_by(v1)%>%
  summarize(mean=mean(v2), sd=sd(v2))

Plotting Results

The figures show that white male and female shooters are less likely to survive as the bars skew more to the left, compared to positive skewness of non-white mass shooters who are male and female.

ggplot(fdshooter, aes(v2)) + geom_histogram() +facet_grid(~v1)

Analysis

In the end, the results shows that a gender difference exists in the fate of mass shooters, with male mass shooters, regardless of race, being more likely to die than female mass shooters. Whites are also more likely to die. Racial variations in gender differences show that white male and female shooters are less likely to survive compared to non-white female and male mass shooters.

It should be noted that much of this analysis hinges on one reality indicated by the introductory data overview: the majority of recorded mass shooters are white and/or male. Naturally, because of this, it is not surprising to see that these groups least likely to survive after a mass shooting (either killed by law enforcement or they take their own lives).

Needless to say, this alone is interesting in a society where mass shootings are becoming more commonplace and the conversation around (domestic) terrorism remains largely centered on so-called perceived “social ills” of particular communities of color (i.e. Muslims).