Md Afzal Hossain | SOC 712 | Home Work 4

Binary Dependent Variable Models

The dataset I will use is stop and frisk data from \(2018\). Every time a police officer stops a person in \(NYC\), the officer is require to fill out a form recording the details of the stop. The forms were filled out by hand and manually entered into an \(NYPD\) database. When the forms became electronic, database are available for download for public . Each row of the data set represent a stop by \(NYPD\) officer. Dataset also contain variable such as the age, race, weight, height of the person stopped, if a person was frisked (when an officer pass the hands over someone in a search for hidden weapons, drugs, or other items and the location of the stop.

I am interested to find likelihood of someone being arrested based on different variables. Therefore, my binary dependent variable will be betting arrested or not.

Prepering Dataset for Binary Dependent Variable:

library(readxl)
# Reading the excel data
sqf<- read_excel("C:/Users/Afzal Hossain/Documents/A.Spring 2019/SOC 712/HW/HW4/sqf_2018.xlsx")
#Keeping the variable that i am interested in.
sq<-select(sqf, 'STOP_FRISK_ID',"MONTH2",'SUSPECT_ARRESTED_FLAG','FRISKED_FLAG','SEARCHED_FLAG','WEAPON_FOUND_FLAG','SUSPECT_REPORTED_AGE':'SUSPECT_BODY_BUILD_TYPE','STOP_LOCATION_BORO_NAME')
# Lets look the Data
head(sq)

We need to handle the missing value that is mention as (null) in categorical variable that we will use later in our model and only consider the recors with complete status in our dataset. I will also recode the race variable in order to easy understanding.

# making the the variable character temporarily
sq$SUSPECT_SEX <- as.character(sq$SUSPECT_SEX)
sq$SUSPECT_SEX[sq$SUSPECT_SEX == "(null)"] <- NA
sq$SUSPECT_RACE_DESCRIPTION <- as.character(sq$SUSPECT_RACE_DESCRIPTION)
sq$SUSPECT_RACE_DESCRIPTION[sq$SUSPECT_RACE_DESCRIPTION == "(null)"]<- NA
sq$SUSPECT_BODY_BUILD_TYPE <- as.character(sq$SUSPECT_BODY_BUILD_TYPE)
sq$SUSPECT_BODY_BUILD_TYPE[sq$SUSPECT_BODY_BUILD_TYPE== "(null)"]<- NA
sq<-sq[complete.cases(sq), ]
# Recoding Variable for bodu type
sq$SUSPECT_BODY_BUILD_TYPE[sq$SUSPECT_BODY_BUILD_TYPE=='EA' | 
                             sq$SUSPECT_BODY_BUILD_TYPE=='HEA' |
                             sq$SUSPECT_BODY_BUILD_TYPE=='THN'|
                             sq$SUSPECT_BODY_BUILD_TYPE=='U' |
                             sq$SUSPECT_BODY_BUILD_TYPE=='XXX'] <- "BIG/Athletic" 
sq$SUSPECT_BODY_BUILD_TYPE[sq$SUSPECT_BODY_BUILD_TYPE=='MED' |
                             sq$SUSPECT_BODY_BUILD_TYPE=='THN'] <- "Small" 
# Recoding Variable for Race
sq$SUSPECT_RACE_DESCRIPTION[sq$SUSPECT_RACE_DESCRIPTION=='ASIAN / PACIFIC ISLANDER' | 
                              sq$SUSPECT_RACE_DESCRIPTION=='BLACK HISPANIC'|
                              sq$SUSPECT_RACE_DESCRIPTION=='WHITE HISPANIC' |
                              sq$SUSPECT_RACE_DESCRIPTION=='AMERICAN INDIAN/ALASKAN NATIVE'] <- "LATINO/OTHER" 
# Taking records that are only complete
sq<-sq[complete.cases(sq), ]
dim(sq)

[1] 10776    13

I can see that it reduce the data size by deleting incomplete dataset.

In order to build our model we need to convert our data to apocopate data format and i will create new variable arrast, Frisk and weapon for the from the existing binary variable which contain \(Y/N\) values.

# Making apropate data type
sq$SUSPECT_ARRESTED_FLAG<- as.factor(sq$SUSPECT_ARRESTED_FLAG)
sq$FRISKED_FLAG<- as.factor(sq$FRISKED_FLAG)
sq$WEAPON_FOUND_FLAG<- as.factor(sq$WEAPON_FOUND_FLAG)
sq$SUSPECT_REPORTED_AGE<- as.double(sq$SUSPECT_REPORTED_AGE)
sq$SUSPECT_SEX<- as.factor(sq$SUSPECT_SEX)
sq$SUSPECT_RACE_DESCRIPTION<- as.factor(sq$SUSPECT_RACE_DESCRIPTION)
sq$SUSPECT_BODY_BUILD_TYPE<- as.factor(sq$SUSPECT_BODY_BUILD_TYPE)
sq$SUSPECT_HEIGHT<- as.double(sq$SUSPECT_HEIGHT)
sq$SUSPECT_WEIGHT<- as.double(sq$SUSPECT_WEIGHT)
# Creating New Variables
sq <- sq %>% 
  mutate(arrast = as.integer(SUSPECT_ARRESTED_FLAG),Frisk = as.integer(FRISKED_FLAG), weapon =as.integer(WEAPON_FOUND_FLAG))
sq <- sq %>% 
  select(arrast, Frisk, weapon,SUSPECT_SEX,SUSPECT_REPORTED_AGE,SUSPECT_RACE_DESCRIPTION, everything())
head(sq)

I will now recode the variable to \(0/1\) instead of \(Y/N\) for my newly created variables which are arrast, Frisk and weapon. Also Filter any missing value.

sq <- sq %>% filter(!is.na(SUSPECT_REPORTED_AGE),!is.na(SUSPECT_HEIGHT),!is.na(SUSPECT_WEIGHT),!is.na(SUSPECT_SEX))%>%
  mutate(arrasted = sjmisc::rec(arrast, rec = "1=0; 2=1"), Frisked = sjmisc::rec(Frisk,   rec = "1=0; 2=1"),
         weapon_found = sjmisc::rec(weapon, rec = "1=0; 2=1")) %>% 
             select(arrasted, SUSPECT_ARRESTED_FLAG, Frisked,FRISKED_FLAG,   
                    weapon_found,WEAPON_FOUND_FLAG,SUSPECT_SEX,SUSPECT_REPORTED_AGE,
                          SUSPECT_RACE_DESCRIPTION, everything()) %>% 
                                select(-arrast)
head(sq)

Logistic Regression

Lets create our First simple model where we predict using our bineary dependent variable arrasted and a independent continuas variable Age of the suspect stop by officer, SUSPECT_REPORTED_AGE

m0 <- glm(arrasted ~ SUSPECT_REPORTED_AGE, family = binomial, data = sq)
summary(m0)

Call:
glm(formula = arrasted ~ SUSPECT_REPORTED_AGE, family = binomial, 
    data = sq)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.9719  -0.8460  -0.8250   1.5098   1.6037  

Coefficients:
                      Estimate Std. Error z value Pr(>|z|)    
(Intercept)          -1.021811   0.056976 -17.934  < 2e-16 ***
SUSPECT_REPORTED_AGE  0.005945   0.001768   3.363 0.000771 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 11963  on 9788  degrees of freedom
Residual deviance: 11952  on 9787  degrees of freedom
AIC: 11956

Number of Fisher Scoring iterations: 4

From \(m_0\), the coefficient is positive,so we can say that as people get older, chances of them being arrested in crease. From the \(Y-intercept\), I can say that when age is zero, the log odd of being arrested is \(-1.028211\). From the slope value, I can say that for every one-unit increase of age, the log odds of being arrested increases by \(0.005945\). Furthermore, from the \(Z\) \(value\) we can say that it is more than two standard deviation way from zero therefore it is statistically significant and this is confirm by the small \(p-value\). Therefore, suspect age is significant for predicting arrest.

Improving Model:

I will add suspect gender and race information to our second model and expect to find statistically not significant since, a police office are not allow to consider race and gender to arrest someone.

m1 <- glm(arrasted ~ SUSPECT_SEX  + SUSPECT_REPORTED_AGE + SUSPECT_RACE_DESCRIPTION , family = binomial, data = sq)
summary(m1)

Call:
glm(formula = arrasted ~ SUSPECT_SEX + SUSPECT_REPORTED_AGE + 
    SUSPECT_RACE_DESCRIPTION, family = binomial, data = sq)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.9980  -0.8538  -0.8180   1.5046   1.6264  

Coefficients:
                                      Estimate Std. Error z value Pr(>|z|)    
(Intercept)                          -0.960084   0.090326 -10.629  < 2e-16 ***
SUSPECT_SEXMALE                      -0.112396   0.075536  -1.488 0.136756    
SUSPECT_REPORTED_AGE                  0.005977   0.001775   3.367 0.000759 ***
SUSPECT_RACE_DESCRIPTIONLATINO/OTHER  0.114725   0.047819   2.399 0.016433 *  
SUSPECT_RACE_DESCRIPTIONWHITE         0.003273   0.077097   0.042 0.966135    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 11963  on 9788  degrees of freedom
Residual deviance: 11944  on 9784  degrees of freedom
AIC: 11954

Number of Fisher Scoring iterations: 4

We can see that age remain significant variable. However, we can see only been lotion have some significant when cones being arrested, but no other race and gender are statistically significant.

Lets Add Intersection and More:

Let me add more variable and intersection between weapon found and been frisked(searched for drug, weapon etc.)

m2 <- glm(arrasted ~ SUSPECT_SEX +SUSPECT_REPORTED_AGE + weapon_found*Frisked +SUSPECT_RACE_DESCRIPTION + SUSPECT_HEIGHT +SUSPECT_BODY_BUILD_TYPE, family = binomial, data = sq)
summary(m2)

Call:
glm(formula = arrasted ~ SUSPECT_SEX + SUSPECT_REPORTED_AGE + 
    weapon_found * Frisked + SUSPECT_RACE_DESCRIPTION + SUSPECT_HEIGHT + 
    SUSPECT_BODY_BUILD_TYPE, family = binomial, data = sq)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9317  -0.7796  -0.7237   0.7718   1.8328  

Coefficients:
                                      Estimate Std. Error z value Pr(>|z|)    
(Intercept)                          -1.290108   0.350854  -3.677 0.000236 ***
SUSPECT_SEXMALE                      -0.226927   0.080349  -2.824 0.004739 ** 
SUSPECT_REPORTED_AGE                  0.006213   0.001896   3.277 0.001050 ** 
weapon_found                          1.876116   0.195850   9.579  < 2e-16 ***
Frisked                              -0.182285   0.050219  -3.630 0.000284 ***
SUSPECT_RACE_DESCRIPTIONLATINO/OTHER  0.139205   0.050820   2.739 0.006159 ** 
SUSPECT_RACE_DESCRIPTIONWHITE         0.080195   0.080718   0.994 0.320454    
SUSPECT_HEIGHT                        0.060340   0.062626   0.964 0.335293    
SUSPECT_BODY_BUILD_TYPESmall         -0.168502   0.049008  -3.438 0.000585 ***
weapon_found:Frisked                  0.468597   0.215828   2.171 0.029919 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 11963  on 9788  degrees of freedom
Residual deviance: 11017  on 9779  degrees of freedom
AIC: 11037

Number of Fisher Scoring iterations: 4

I can see that weapon, frisked significant coefficient that is expected since if weapon found, chance of getting arrested are high. Furthermore, intersection of weapon found and frisk and body type are very important coefficient as well. The age remain significant, and gender become an important coefficient. Therefore, log odds of been arrested change based those coefficient. It is interesting to see that Male has negative log odds of been arrested.

Model Selection

We will use likelihood ratio test and information criteria for our best model selection

Likelihood Ratio Test:

anova(m0, m1, m2, test = "Chisq")

Analysis of Deviance Table

Model 1: arrasted ~ SUSPECT_REPORTED_AGE
Model 2: arrasted ~ SUSPECT_SEX + SUSPECT_REPORTED_AGE + SUSPECT_RACE_DESCRIPTION
Model 3: arrasted ~ SUSPECT_SEX + SUSPECT_REPORTED_AGE + weapon_found * 
    Frisked + SUSPECT_RACE_DESCRIPTION + SUSPECT_HEIGHT + SUSPECT_BODY_BUILD_TYPE
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)    
1      9787      11952                         
2      9784      11944  3      8.3  0.04019 *  
3      9779      11018  5    926.4  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Form this test, we can see that model 3 which is \(m_2\) has a significantly improved based on very small P-value. Now we will use Information criteria to find which model is better.

Information Criteria:

Bayesian information criterion (BIC): \[BIC = -2L + k\times log(n)\]

Akaike information criterion (AIC): \[AIC = 2(L-k)\]

library(texreg)
#htmlreg(list(m0, m1, m2), doctype = FALSE)

Statistical models

	Model 1	Model 2	Model 3
(Intercept)	-1.02***	-0.96***	-1.29***
	(0.06)	(0.09)	(0.35)
SUSPECT_REPORTED_AGE	0.01***	0.01***	0.01**
	(0.00)	(0.00)	(0.00)
SUSPECT_SEXMALE		-0.11	-0.23**
		(0.08)	(0.08)
SUSPECT_RACE_DESCRIPTIONLatino/Other		0.11*	0.14**
		(0.05)	(0.05)
SUSPECT_RACE_DESCRIPTIONWHITE		0.00	0.08
		(0.08)	(0.08)
weapon_found			1.88***
			(0.20)
Frisked			-0.18***
			(0.05)
SUSPECT_HEIGHT			0.06
			(0.06)
SUSPECT_BODY_BUILD_TYPESmall			-0.17***
			(0.05)
weapon_found:Frisked			0.47*
			(0.22)
AIC	11956.17	11953.86	11037.47
BIC	11970.54	11989.81	11109.36
Log Likelihood	-5976.08	-5971.93	-5508.73
Deviance	11952.17	11943.86	11017.47
Num. obs.	9789	9789	9789
p < 0.001, p < 0.01, p < 0.05

From both \(AIC\) and \(BIC\), we can see that models get better as their value decrease with more variable and complexity. The best performing model based on both of these metrics is Model 3 which is model \(m_2\)

Visualizations Our Best Model:

library(visreg)
visreg(m2, "SUSPECT_REPORTED_AGE",scale="response",line=list(col="magenta2"),fill=list(col="burlywood2"),
       xlab="Suspect Age") +  theme_bw()

Conditions used in construction of plot
SUSPECT_SEX: MALE
weapon_found: 0
Frisked: 1
SUSPECT_RACE_DESCRIPTION: BLACK
SUSPECT_HEIGHT: 5.7
SUSPECT_BODY_BUILD_TYPE: BIG/Athletic
NULL

We can see that, as age increase probability of being arrested slightly increased.

visreg(m2, "SUSPECT_SEX", by = "weapon_found", scale = "response",line=list(col="magenta2"),
                             fill=list(col="navyblue"), xlab="Suspect Sex")

For above, we can see that when weapon found and not found, the probability of getting arrest for Female is slightly high. This may explain earlier result where we saw male has negative log odd of being arrested. This result surprise me, but the difference is very insignificant and given that male sample is larger than female may contribute to this slight different.

visreg(m2, "SUSPECT_RACE_DESCRIPTION", by = "Frisked",scale ="response",line=list(col="deeppink1"),
                             fill=list(col="darkcyan"), xlab="Suspect Race")

From this plot, we see that when comes to weather different race was frisked (searched for drug, weapon etc.) or not, Black people are less likely to get arrested compare to white.But they are the one stop most compare to other race almost \(56\%\) of stop were Black. Therefore, we might want to look into this more to find weather there is a racial bias exist or not. Latinos/ other race has highest probability of being arrested in both scenario. May be they tend to have burgs or illegal substances compare to other demography.

visreg(m2, "SUSPECT_BODY_BUILD_TYPE", by = "Frisked", scale = "response",line=list(col="darkgoldenrod1"),
                             fill=list(col="firebrick4"), xlab="Suspect Body Type")

This plot clearly show that whether they were frisked or not people with big body are more likely to get arrested. Maybe because office feel threaten by their body type.

Checking Our Model Performance:

Now let’s check how our model actually did on our data. Therefore, we will draw a graph that will show the predicted probabilities of each arrest with their actual arrest status.

First, we create a new data frame that contains the probabilities of being arrested along with actual arrest status.Then we sort the data frame from low to high probabilities. Then we add a new column to the data frame that has the rank of each sample, from low probability to high probability.

library(ggplot2)
#Create a new data frame that contains the probabilities of being arrested along with actual arrest status
pred_data<-data.frame( predict_of_arrested = m2$fitted.values, Arrested = sq$SUSPECT_ARRESTED_FLAG)
#Sort the data frame from low to high probabilities
pred_data <- pred_data[order(pred_data$predict_of_arrested,decreasing = FALSE),]
#add a new column to the data frame that has the `rank` of each sample
pred_data$rank <- 1:nrow(pred_data)
ggplot(data= pred_data,aes(x=rank,y=predict_of_arrested))+
  geom_point(aes(color =Arrested),alpha =1,shape=5, stroke =2)+
  xlab("Index") + ylab("Probability of Getting Arrested") +theme_minimal()

From our graph, we can see that most of the people how got arrested (one in turquoise/blue) are predicted to have a high probability of getting arrest. Other hand, most of the people how did not get arrested (the ones in red/salmon) are predicted to have low probability of getting arrested.

Which mean our mode did a decent job of predicting arrest. However,it is not perfect. We can see that some of people how were not arrested are shown high probability of being arrested(red/salmon color in top). In future we might want to add more variable like neighborhood in our model since police stop usually convert to arrest in high crime area more frequently. In addition,in future we might want to split our data to test and train data set to actually test our model to a dataset it never saw.