Soc 712: Homework 6

Jennifer Ganeles
3/24/19

Part 1: Exploring the Zelig Package (Zelig 4 vs. Zelig 5 Syntax)

Using the Titanic dataset from the radiant.data package, I re-implemented the regression analyses conducted in class (slides 5.2-5.18) using the Zelig 5 syntax.

Preparing the Data: Recoding and Selecting Variables

library(radiant.data)
data(titanic)
titanic=titanic %>%
  mutate(surv = as.integer(survived)) %>% 
  mutate(survival = sjmisc::rec(surv, rec = "2=0; 1=1")) %>% 
  select(pclass, survived, survival, everything())
head(titanic)

Slide 5.2: Interaction Model

#Zelig 4: 
#z_tit <- zelig(survival ~ age + sex*pclass + fare, model = "logit", data = titanic, cite = F)
#Zelig 5:
z5_tit <- zlogit$new()
z5_tit$zelig(survival ~ age + sex*pclass + fare, data = titanic)
summary(z5_tit)

Model: 

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

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.0634  -0.6641  -0.4943   0.4336   2.4941  

Coefficients:
                    Estimate Std. Error z value Pr(>|z|)
(Intercept)        4.8978215  0.6131092   7.988 1.37e-15
age               -0.0387245  0.0068044  -5.691 1.26e-08
sexmale           -3.8996177  0.5015659  -7.775 7.55e-15
pclass2nd         -1.5923247  0.6024844  -2.643  0.00822
pclass3rd         -4.1382715  0.5601819  -7.387 1.50e-13
fare              -0.0009058  0.0020509  -0.442  0.65874
sexmale:pclass2nd -0.0600076  0.6372949  -0.094  0.92498
sexmale:pclass3rd  2.5019110  0.5479901   4.566 4.98e-06

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1409.99  on 1042  degrees of freedom
Residual deviance:  931.45  on 1035  degrees of freedom
AIC: 947.45

Number of Fisher Scoring iterations: 5

Next step: Use 'setx' method

Slides 5.3-5.4: Age Effect

library(Zelig)
#Zelig 4: 
#a.range = min(titanic$age):max(titanic$age)
#x <- setx(z_tit, age = a.range)
#s <- sim(z_tit, x = x)
#ci.plot(s)
#Zelig 5:
a.range=min(titanic$age):max(titanic$age)
z5_tit$setrange(age=a.range)
z5_tit$sim()
z5_tit$graph()

Looking at the full age range, we find that infants have a higher expected value of survivial (median~.35) than people who are older, but uncertainty is higher for infants than for those who are older (we are more certain in low survival rate when older)

Slides 5.5-5.6: Fare Effect

#Zelig 4:
#f.range = min(titanic$fare):max(titanic$fare)
#x <- setx(z_tit, fare = f.range)
#s <- sim(z_tit, x = x)
#ci.plot(s)
#Zelig 5:
f.range = min(titanic$fare):max(titanic$fare)
z5_tit$setrange(fare=f.range)
z5_tit$sim()
z5_tit$graph()

Looking at the full fare range, we find that the expected value of surviving is almost constant as fare increases (with a very minor decrease in survival and a large increase in uncertainty as fare increases). This suggests that fare does not really influence survival chance (after controlling for class).

Slides 5.7-5.9: Sex Difference

#Zelig 4:
#x <- setx(z_tit, sex = "male")
#x1 <- setx(z_tit, sex = "female")
#s <- sim(z_tit, x = x, x1 = x1)
#fd <- s$get_qi(xvalue="x1", qi="fd")
#summary(fd)
#plot(s)
#Zelig 5:
z5_tit$setx(sex="male")
z5_tit$setx1(sex="female")
z5_tit$sim()
fd<-z5_tit$get_qi(xvalue="x1", qi="fd")
summary(fd)

       V1        
 Min.   :0.1364  
 1st Qu.:0.2291  
 Median :0.2556  
 Mean   :0.2578  
 3rd Qu.:0.2859  
 Max.   :0.3928  

Looking at the simulated first difference between males and females in expected survival, we find that on average, females have a .26 higher chance of surviving than males (can be as low as .14 and as high as .39)

z5_tit$graph()

Slides 5.10-5.15: Class Variation in Sex Difference

#Zelig 4:
#c1x <- setx(z_tit, sex = "male", pclass = "1st")
#c1x1 <- setx(z_tit, sex = "female", pclass = "1st")
#c1s <- sim(z_tit, x = c1x, x1 = c1x1)
#Zelig 5:
z5.1 <- zlogit$new()
z5.1$zelig(survival ~ age +sex*pclass + fare, data=titanic)
z5.1$setx(sex="male", pclass="1st")
z5.1$setx1(sex="female", pclass="1st")
z5.1$sim()
z5.1$graph()

Slide 5.12:
For first class, the difference in expected value of survival between females and males is about .52, suggesting that first class females have a .52 higher chance of surviving than first class males.

#Zelig 4:
#c2x <- setx(z_tit, sex = "male", pclass = "2nd")
#c2x1 <- setx(z_tit, sex = "female", pclass = "2nd")
#c2s <- sim(z_tit, x = c2x, x1 = c2x1)
#Zelig 5:
z5.2 <- zlogit$new()
z5.2$zelig(survival ~ age +sex*pclass + fare, data=titanic)
z5.2$setx(sex="male", pclass="2nd")
z5.2$setx1(sex="female", pclass="2nd")
z5.2$sim()
z5.2$graph()

Slide 5.13:
For second class, the simulated first difference between male and female survival increases to about .75, suggesting that second class females have a .75 higher chance of survial than second class males.

#Zelig 4:
#c3x <- setx(z_tit, sex = "male", pclass = "3rd")
#c3x1 <- setx(z_tit, sex = "female", pclass = "3rd")
#c3s <- sim(z_tit, x = c3x, x1 = c3x1)
#Zelig 5:
z5.3 <- zlogit$new()
z5.3$zelig(survival ~ age +sex*pclass + fare, data=titanic)
z5.3$setx(sex="male", pclass="3rd")
z5.3$setx1(sex="female", pclass="3rd")
z5.3$sim()
z5.3$graph()

Slide 5.14:
For third class, the simulated first difference between female and male survival drops to about .26, suggesting that third class women have a .26 higher chance of survival than third class men. (As we move down the class ranks, women’s expected survival decreases.)

Slide 5.15: Putting Them in One Place

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

Slide 5.16: Changing from Wide to Long Format

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

Slide 5.17: Group by Class

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

Slide 5.18: Plotting

Looking at the expected first difference between female and male passengers across all 3 classes, we can see the first difference is highest in second class and lowest in third class. Male and female passengers have closer expected survival rates in third class.

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

Part 2: Exploring the Zelig Package Using My Own Dataset (Marvel Characters as Living or Deceased)

Using a FiveThirtyEight Comic Characters Dataset found on Kaggle, this week’s homework seeks to determine how variables such as alignment (good, bad, or neutral), gender (male or female), identity (public or secret), and number of appearances (as of 2014) explain the likelihood of whether a Marvel character would be alive in the comics or not. Since my dependent variable is dichotomous, I conducted a logit regression analysis.

Preparing the Data: Recoding and Selecting Variables

My dependent variable was recoded as binary (alive_binary), where 1=living and 0=deceased. My identity variable was also recoded so that “Public Identity”" included “No Dual Identity” and “Known to Authorities.” For my gender variable, “Agender” and “Genderfluid” characters were recoded as “Other.”

library(readr)
library(dplyr)
library(magrittr)
Marvel<-read_csv("/Users/jenniferganeles/Downloads/marvel-wikia-data.csv")
marvel_heroes<-
  mutate (Marvel, alive_binary=
           recode(ALIVE, "Living Characters"=1,"Deceased Characters"=0),
          Gender=
            recode(SEX, "Agender Characters"="Other", "Genderfluid Characters"="Other"),
          Identity=
            recode(ID, "No Dual Identity"="Public Identity", "Known to Authorities Identity"="Public Identity"),
          alignment = factor(ALIGN),
          gender = factor(Gender),
          identity=factor(Identity))%>%
    filter(!is.na(APPEARANCES),
           !is.na(alignment))%>%
  
  select(ALIVE, alive_binary, alignment, gender, identity, APPEARANCES)
head(marvel_heroes)

Running Three Regression Models and Determining the Best Model Fit

Both AIC and BIC agree that Model 2 is the best model, which contains all selected variables but without an interaction term.

M1<-zlogit$new()
M1$zelig(alive_binary~ identity+gender, data=marvel_heroes)
M2<-zlogit$new()
M2$zelig(alive_binary~ identity+gender+alignment+APPEARANCES, data=marvel_heroes)
M3<-zlogit$new()
M3$zelig(alive_binary~ identity+gender*APPEARANCES+alignment, data=marvel_heroes)
texreg::htmlreg(list(M1, M2, M3),doctype = FALSE)

Statistical models

	Model 1	Model 2	Model 3
(Intercept)	1.38***	0.94***	0.94***
	(0.05)	(0.07)	(0.07)
identitySecret Identity	-0.13**	-0.02	-0.02
	(0.05)	(0.05)	(0.05)
genderMale Characters	-0.37***	-0.23***	-0.23***
	(0.06)	(0.06)	(0.06)
genderOther	0.54	0.69	0.42
	(0.54)	(0.54)	(0.65)
alignmentGood Characters		0.63***	0.63***
		(0.05)	(0.05)
alignmentNeutral Characters		0.31***	0.31***
		(0.07)	(0.07)
APPEARANCES		0.00**	0.00
		(0.00)	(0.00)
genderMale Characters:APPEARANCES			-0.00
			(0.00)
genderOther:APPEARANCES			0.03
			(0.06)
AIC	11511.20	11354.41	11357.65
BIC	11540.06	11404.92	11422.58
Log Likelihood	-5751.60	-5670.21	-5669.82
Deviance	11503.20	11340.41	11339.65
Num. obs.	10047	10047	10047
p < 0.001, p < 0.01, p < 0.05

*The above table displays the regression results in terms of log odds, but I want to produce more understandable coefficients using simulation…

Bad vs. Good Characters:

Looking at the simulated first difference between bad and good characters, my analysis shows that on average, good characters have a .12 higher chance of being alive than bad characters (first difference can be as low as .09 and as high as .15).

M2$setx(alignment="Bad Characters")
M2$setx1(alignment="Good Characters")
M2$sim()
fd <- M2$get_qi(xvalue="x1", qi="fd")
summary(fd)

       V1         
 Min.   :0.08721  
 1st Qu.:0.11400  
 Median :0.12094  
 Mean   :0.12105  
 3rd Qu.:0.12795  
 Max.   :0.15065  

M2$graph()

Bad vs. Neutral Characters:

Looking at the simulated first difference between bad and neutral characters, my analysis shows that on average, neutral characters have a .06 higher chance of being alive than bad characters (first difference can be as low as .02 and as high as .10).

M2$setx(alignment="Bad Characters")
M2$setx1(alignment="Neutral Characters")
M2$sim()
fd <- M2$get_qi(xvalue="x1", qi="fd")
summary(fd)

       V1         
 Min.   :0.02291  
 1st Qu.:0.05597  
 Median :0.06499  
 Mean   :0.06473  
 3rd Qu.:0.07356  
 Max.   :0.10435  

M2$graph()

Male vs. Female Characters:

Looking at the first difference between male and female characters, my analysis shows that on average, female characters have a .05 higher chance of being alive than male characters (first difference can be as low as .01 and as high as .09).

M2$setx(gender="Male Characters")
M2$setx1(gender="Female Characters")
M2$sim()
fd <- M2$get_qi(xvalue="x1", qi="fd")
summary(fd)

       V1         
 Min.   :0.01015  
 1st Qu.:0.04152  
 Median :0.04910  
 Mean   :0.04912  
 3rd Qu.:0.05680  
 Max.   :0.08541  

M2$graph()

Secret vs. Public Identity:

Here, the analysis shows that on average, characters with public identities have a .004 higher chance of being alive than characters with secret identities. However, this difference is not significant (as the AIC/BIC table showed above) and the first difference can be as high as .036 and as low as -.029 (meaning, characters with public identities can also have a lower percentage of being alive). Therefore, whether a character’s identity is secret or public does not seem to have much influence on their chances of being alive.

M2$setx(identity="Secret Identity")
M2$setx1(identity="Public Identity")
M2$sim()
fd <- M2$get_qi(xvalue="x1", qi="fd")
summary(fd)

       V1           
 Min.   :-0.028765  
 1st Qu.:-0.002651  
 Median : 0.004200  
 Mean   : 0.004119  
 3rd Qu.: 0.011159  
 Max.   : 0.035572  

M2$graph()

Number of Appearances Effect:

Looking at a range between 0 and 1000 appearances, I found that the expected value of being alive increases as the number of appearances in the comic increases. This makes sense since characters would likely appear in the comics as living (rather than as deceased). However, the uncertainty of being alive also increases as number of appearances increase.

M2$setrange(APPEARANCES=0:1000)
M2$sim()

|================================================================|100% ~0 s remaining     

M2$graph()

Conclusion:

After running a logistic regression using Zelig, I was able to determine the chances of a Marvel character being alive or not by exploring independent variables such as alignemnt, gender, identity, and appearance number. I found that characters who are good or neutral are more likely to be alive in the comics than those who are bad, and characters who are female are more likely to be alive than those who are male. Whether a character’s identity is secret or public does not have much influence on whether he/she is alive in the comics. Lastly, characters have a higher chance of being alive as they appear more often in the comics; however, as number of appearances increases, the uncertainty of being alive increases as well.