Md Afzal Hossain | Soc 712 |Home Work 6

Part-1:

Re-implement the computations described on pages 5.2-5.18 of the lecture slides using the Zelig 5 syntax

Reading data and making binary variable:

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

Building model using zelig(Slide 5.2):

library(Zelig)
z_tit <- zlogit$new()
z_tit$zelig(survival ~ age + sex*pclass + fare, data = titanic)
summary(z_tit)

Model: 

Call:
z_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

Age effect (Slide 5.3-5.4)

z_tit$setrange(age=min(titanic$age):max(titanic$age))
z_tit$sim()
z_tit$graph()

Fare effect (Slide 5.5-5.6)

z_tit$setrange(fare=min(titanic$fare):max(titanic$fare))
z_tit$sim()
z_tit$graph()

Sex difference | Only simulated first difference(Slide 5.7-5.9)

z_tit$setx(sex="male")
z_tit$setx1(sex="female")
z_tit$sim()
fd <- z_tit$get_qi(xvalue="x1", qi="fd")
summary(fd)

       V1        
 Min.   :0.1259  
 1st Qu.:0.2269  
 Median :0.2542  
 Mean   :0.2556  
 3rd Qu.:0.2841  
 Max.   :0.3905  

z_tit$graph()

Testing the class variations in sex difference using Zelig (Slide 5.10-5.14)

First Class:

z_tit1 <- zlogit$new()
z_tit1$zelig(survival ~ age +sex*pclass + fare, data=titanic)
#z_tit1<-z_tit
z_tit1$setx(sex="male", pclass="1st")
z_tit1$setx1(sex="female", pclass="1st")
z_tit1$sim()
z_tit1$graph()

Second Class:

z_tit2 <- zlogit$new()
z_tit2$zelig(survival ~ age +sex*pclass + fare, data=titanic)
z_tit2$setx(sex="male", pclass="2nd")
z_tit2$setx1(sex="female", pclass="2nd")
z_tit2$sim()
z_tit2$graph()

Third Class:

z_tit3 <- zlogit$new()
z_tit3$zelig(survival ~ age +sex*pclass + fare, data=titanic)
z_tit3$setx(sex="male", pclass="3rd")
z_tit3$setx1(sex="female", pclass="3rd")
z_tit3$sim()
z_tit3$graph()

Putting them in one place and making it long format (Slide 5.15-5.18)

d1 <- z_tit1$get_qi(xvalue="x1", qi="fd")
d2 <- z_tit2$get_qi(xvalue="x1", qi="fd")
d3 <- z_tit3$get_qi(xvalue="x1", qi="fd")
  
dfd <- as.data.frame(cbind(d1, d2, d3))
tidd <- dfd %>% 
  gather(class, simv, 1:3)
head(tidd)

Using Group By

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

Plotting:

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

Part 2: Conducting regression analysis using Zelig on my own data set

Binary Dependent Variable and Data:

The dataset i will use includes all valid felony, misdemeanor, and violation crimes reported to the New York City Police Department (NYPD) for all complete quarters so far in year 2017. I am only interested in the likelihood of whether crime was successfully completed(successfully accomplished by criminal) or attempted(not completed). Therefore, my binary dependent variable will be whether the crime was completed or not.

Prepering Dataset for Binary Dependent Variable:

ny_c<-read.csv("C:/Users/Afzal Hossain/Documents/A.Spring 2019/SOC 712/HW/HW6/NYPD_Complaint_Data_Current__Year_To_Date_.csv")
ny_cr<-select(ny_c, 'CRM_ATPT_CPTD_CD',"VIC_RACE",'VIC_SEX',"VIC_AGE_GROUP",'BORO_NM','LAW_CAT_CD')
dim(ny_cr)

[1] 464065      6

#Recoding Age
ny_cr$VIC_AGE_GROUP<-as.character(ny_cr$VIC_AGE_GROUP)
ny_cr$VIC_AGE_GROUP[ny_cr$VIC_AGE_GROUP=='<18']<-"CHILD"
ny_cr$VIC_AGE_GROUP[ny_cr$VIC_AGE_GROUP=='18-24' |
                  ny_cr$VIC_AGE_GROUP=='25-44']<-'YOUNG'
ny_cr$VIC_AGE_GROUP[ny_cr$VIC_AGE_GROUP=='45-64']<-'MIDDLE AGE'
ny_cr$VIC_AGE_GROUP[ny_cr$VIC_AGE_GROUP=='65+']<-'OLD'
# Recoding Variable for Sex
ny_cr$VIC_SEX<-as.character(ny_cr$VIC_SEX) 
ny_cr$VIC_SEX[ny_cr$VIC_SEX=='E' |
                   ny_cr$VIC_SEX=='D' | ny_cr$VIC_SEX== 'U']<-NA
# Recoding Variable for Race
ny_cr$VIC_RACE<-as.character(ny_cr$VIC_RACE)
ny_cr$VIC_RACE[ny_cr$VIC_RACE=='ASIAN / PACIFIC ISLANDER' | 
                              ny_cr$VIC_RACE=='BLACK HISPANIC'|
                              ny_cr$VIC_RACE =='BLACK' |
                              ny_cr$VIC_RACE=='WHITE HISPANIC' |
                              ny_cr$VIC_RACE=='AMERICAN INDIAN/ALASKAN NATIVE' |
                              ny_cr$VIC_RACE== 'UNKNOWN'] <- "NON WHITE" 
#Changing data type
ny_cr$VIC_RACE<- as.factor(ny_cr$VIC_RACE)
ny_cr$VIC_SEX<- as.factor(ny_cr$VIC_SEX)
ny_cr$BORO_NM<- as.factor(ny_cr$BORO_NM)
ny_cr$LAW_CAT_CD<- as.factor(ny_cr$LAW_CAT_CD)
head(ny_cr)

ny_cr%>% select(VIC_AGE_GROUP)%>% distinct()

Victims age has a lot of negative and unusual number because of human error. Therefore, we will only take our four category that we created.

ny_cr <- subset(ny_cr, VIC_AGE_GROUP == "CHILD" | VIC_AGE_GROUP == "YOUNG" |VIC_AGE_GROUP == "MIDDLE AGE" |VIC_AGE_GROUP == "OLD")
dim(ny_cr)

[1] 332277      6

I lost around ten thousand record but I am not concern since, it is still a large sample.

#Creating Binary variable 1/0 instead of Completed/Attempted
ny_cr <- ny_cr %>% 
  mutate(Crime_Completed =as.integer(CRM_ATPT_CPTD_CD)) %>% 
  select(Crime_Completed,CRM_ATPT_CPTD_CD, everything())
ny_cr <- ny_cr%>% 
          filter(!is.na(VIC_RACE),!is.na(VIC_SEX),!is.na(LAW_CAT_CD),!is.na(CRM_ATPT_CPTD_CD))%>%
                 mutate(Crm_Completed = sjmisc::rec(Crime_Completed, rec = "1=0; 2=1")) %>% 
                          select(Crm_Completed,CRM_ATPT_CPTD_CD,everything())
ny_cr$VIC_AGE_GROUP<- as.factor(ny_cr$VIC_AGE_GROUP)
ny_cr$Crm_Completed<- as.factor(ny_cr$Crm_Completed)
head(ny_cr)

Building logistic model using Zelig:

z.crime <- zelig(Crm_Completed ~ VIC_RACE + VIC_SEX* LAW_CAT_CD+VIC_AGE_GROUP , model = "logit", data = ny_cr, cite = FALSE)

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

summary(z.crime)

Model: 

Call:
z5$zelig(formula = Crm_Completed ~ VIC_RACE + VIC_SEX * LAW_CAT_CD + 
    VIC_AGE_GROUP, data = ny_cr)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.5260   0.1099   0.1183   0.2544   0.3419  

Coefficients:
                               Estimate Std. Error z value Pr(>|z|)
(Intercept)                     3.12091    0.05796  53.846  < 2e-16
VIC_RACEWHITE                  -0.04482    0.03276  -1.368   0.1713
VIC_SEXM                       -0.26563    0.03258  -8.153 3.54e-16
LAW_CAT_CDMISDEMEANOR           1.69260    0.04970  34.056  < 2e-16
LAW_CAT_CDVIOLATION             2.79991    0.10722  26.114  < 2e-16
VIC_AGE_GROUPMIDDLE AGE         0.19673    0.06072   3.240   0.0012
VIC_AGE_GROUPOLD                0.13482    0.07390   1.824   0.0681
VIC_AGE_GROUPYOUNG              0.29341    0.05712   5.137 2.79e-07
VIC_SEXM:LAW_CAT_CDMISDEMEANOR  0.11751    0.06834   1.720   0.0855
VIC_SEXM:LAW_CAT_CDVIOLATION    0.17115    0.17184   0.996   0.3193

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 54163  on 332187  degrees of freedom
Residual deviance: 49223  on 332178  degrees of freedom
AIC: 49243

Number of Fisher Scoring iterations: 8

Next step: Use 'setx' method

I don’t see any statistically significant coefficient which is not surprising to me since, whether a crime completed or just attempted not really depend on age, gender or race. Its more to do with whither someone call police or in some situation, victim is physically superior to the criminal and can defend themself.

But from here, I can say that log odds of crime event completed is less when the victim is Male compare to Female (-.265). Furthermore, when the crime is just a violation compare to Misdemeanor, the log odd are higher that crime will be completed, 2.799 compare to 1.693 and it is true for whither victims are male or female.Age group seems to have no effect, but if victim is old, crime has less log odd of being completed. Now, let’s use simulation-based results interpretation for better understand the result.

Making sense of regression results

First, I will look at sex difference:

x <- setx(z.crime, VIC_SEX = "M")
x1 <- setx(z.crime, VIC_SEX = "F")
c <- sim(z.crime, x = x, x1 = x1)
fd <- c$get_qi(xvalue="x1", qi="fd")
summary(fd)

       V1            
 Min.   :-0.0002070  
 1st Qu.: 0.0007046  
 Median : 0.0009567  
 Mean   : 0.0009613  
 3rd Qu.: 0.0012279  
 Max.   : 0.0020448  

I can see that on average a crime will be completed 0.00096 more if the victims is female.

plot(c)

Looking at race difference:

wx <- setx(z.crime, VIC_RACE = "WHITE")
nwx1 <- setx(z.crime, VIC_RACE = "NON WHITE")
r <- sim(z.crime, x = wx, x1 = nwx1)
fd <- r$get_qi(xvalue ="x1", qi="fd")
summary(fd)

       V1            
 Min.   :-0.0003192  
 1st Qu.: 0.0001266  
 Median : 0.0002658  
 Mean   : 0.0002684  
 3rd Qu.: 0.0004022  
 Max.   : 0.0009140  

If the victims is Non White, on average the crime slightly more likely to be completed(0.00027)

plot(r)

Test the crime variations in sex difference

# Felony
v1x <- setx(z.crime, VIC_SEX = "M", LAW_CAT_CD = "FELONY")
v1x1 <- setx(z.crime, VIC_SEX = "F", LAW_CAT_CD ="FELONY")
v1s <- sim(z.crime, x = v1x, x1 = v1x1)
#Misdemeanor
v2x <- setx(z.crime, VIC_SEX = "M", LAW_CAT_CD = "MISDEMEANOR")
v2x1 <- setx(z.crime, VIC_SEX = "F", LAW_CAT_CD = "MISDEMEANOR")
v2s <- sim(z.crime, x = v2x, x1 = v2x1)
#Violation
v3x <- setx(z.crime, VIC_SEX = "M", LAW_CAT_CD = "VIOLATION")
v3x1 <- setx(z.crime, VIC_SEX = "F", LAW_CAT_CD = "VIOLATION")
v3s <- sim(z.crime, x = v3x, x1 = v3x1)

Plotting Felony Crime:

plot(v1s)

Plotting Misdemeanor Crime:

plot(v2s)

Plotting Violation Crime:

plot(v3s)

Putting all of them in one place:

d1 <- v1s$get_qi(xvalue="x1", qi="fd")
d2 <- v2s$get_qi(xvalue="x1", qi="fd")
d3 <- v3s$get_qi(xvalue="x1", qi="fd")
dfd <- as.data.frame(cbind(d1, d2, d3))
law_categories <- dfd %>% 
  gather(Law, simv, 1:3)
head(law_categories)

Group Different Crime:

law_categories %>% 
  group_by(Law) %>% 
  summarise(mean = mean(simv), sd = sd(simv))

For all three different crime (Felony, Misdemeanor and Violation), if the victims is Female, crime is more likely to be completed, even though it is not statistically significant.

Plotting the Distributions:

ggplot(law_categories, aes(simv)) + geom_histogram() + facet_grid(~Law) +geom_vline(xintercept=mean(law_categories$simv), color="red") +theme_bw()