Background

Import Libraries

library(dplyr)
library(lmtest)
library(ggpmisc)
library(foreach)
library(ggpubr)
library(ggplot2)

Set Plot Theme

theme_set(theme_minimal()+theme(
    plot.title = element_text(size = 12,hjust = 0.5),
    plot.subtitle = element_text(face="bold",
     size = 12)
    ))
update_geom_defaults("col", list(fill = "skyblue"))
update_geom_defaults("point", list(colour = "skyblue3"))
update_geom_defaults("smooth", list(colour = "magenta3", linetype= "dashed" , 
                                    linewidth = .75))

Plot Helper Functions

stat_poly_eq2 <- function(...) {
  stat_poly_eq(formula = y ~ x, coef.digits = 10,rr.digits = 10,p.digits = 10, 
               eq.x.rhs = "~`*`~Age", eq.with.lhs = "`#`~Charges~~`=`~~", ...)
}


plot_start <- function(df,w,...) {
  ggplot(df, aes(x = Age, y = AvNumCharges))+labs(y = "Average # Charges", x="Age (years)",
          title ="Age-Averaged Number\n of Charges per Arrest")+w+
    geom_smooth(method = "lm",formula = y ~ x, ...) +
     stat_poly_eq2(use_label("R2"), label.y = 0.9, label.x = 0.05)
}

 

1) The Data

This project utilizes police arrest data generated from the police department (Record Management System (RMS)) in Scottsdale, Arizona. 

df=read.csv("Police_Arrests.csv") 
names(df)
##  [1] "OBJECTID"           "Arrestee"           "Race"              
##  [4] "Sex"                "Age"                "District"          
##  [7] "Beat"               "Zone"               "ArrestDate"        
## [10] "ArrestTime"         "ArrestNum"          "DrNum"             
## [13] "ArrestType"         "ArizonaStatuteCode" "ChargeDescription" 
## [16] "ArrestLocation"     "CityOfArrestee"     "OfficerSerialNum"

Each row represents one charge, and charges belonging to the same arrest incident share an arrest number.

 

1.1) Data Preprocessing

For each arrest number: combine all charges into an overall count and retain arrest details. 

df=df %>% 
  mutate(
      NumCharges=n(), 
      .by= ArrestNum) %>% 
filter(!duplicated(df$ArrestNum))%>% 
  select(-c("ChargeDescription","ArizonaStatuteCode"))

 

1.2) Data Goal

Fit a linear regression model to estimate the number of charges from other arrest data.

In this analysis the arrestees age will be used to estimate the number of charges. Records with uncommon ages are removed. 

df=df %>% select(c("NumCharges","Age"))%>% 
  filter(Age<=83 & Age > 23)

 

1.3) Data Grouping

Now group the arrest records by age and compute the mean number of charges for each group. 

df2=df %>% summarise(AvNumCharges = mean(NumCharges),.by=Age) 
head(df2)

 

2) Data Exploration

plot_start(df2,geom_point()) + expand_limits(y = 2.5)+stat_poly_eq2(use_label("eq")) + 
            stat_poly_eq2( use_label("P"), label.y = .85, label.x = .05)

 

Heteroskedasticity appears to be present in the model.

Heteroskedasticity happens when the error values from a regression model have unequal variability across the range of another variable. This violates the mathematical assumptions of least ordinary squares regression.

 

2.1) Heteroskedasticity Check 1: Residual Analysis

Check for heteroskedasticity by plotting the residuals of the linear model against the explanatory variable, in this case age. 

df2$ols=lm(AvNumCharges ~ Age, data = df2)
df2$Residuals <- df2$ols$residuals
ggplot(data = df2, aes(y = Residuals, x = Age)) + 
  geom_point() + 
  geom_abline(slope = 0)

 

There appears to be more variation in the residual value for older adults, further suggesting heteroskedasticity is present.

 

2.2) Heteroskedasticity Check 2: Breusch-Pagan Test

To formally detect heteroskedasticity, conduct a Breusch-Pagan test. 

pval=bptest(df2$ols)$p.value
print(pval)
##           BP 
## 9.865546e-05

 

The p-value 9.8655461^{-5} is less than \(\alpha=0.05\), suggesting that heteroskedasticity is present.

 

3) Data Restructuring

This data set can be restructured to minimize heteroskedasticity.

Data binning is the restructing technique chosen due to its ability to average out local charges in variability.

Steps:

  1. Group the records based on arestees age into discrete intervals (bins) then compute the mean number of charges for each bin.

  2. Compute the R-Squared statistic of linear fit and test the model for heteroskedasticity with the Breusch-Pagan test (BP Score).

  3. Repeat for each bin size under consideration.

results <- foreach(i = 1:6) %do% {
  data = df %>% mutate(Age = Age-(Age %% i)) %>%
                summarise(AvNumCharges = mean(NumCharges),.by=Age)
  bpscore=round(bptest(lm(AvNumCharges ~ Age, data = data))$p.value,5)
  list(plt = plot_start(data, geom_col(),se=F)+
         labs(subtitle=paste("(bin width:",i,"years)\nBP Score:", bpscore)) + 
         expand_limits(y = 3),
       bpscore = bpscore)
}

ggarrange(plotlist=lapply(results, function(x) x$plt),
          labels = c("A)", "B)", "C)", "D)", "E)", "F)"),
          ncol = 3, nrow = 2)

 

The bin width is set to 4 years. 

df=df %>% mutate(Age = Age-(Age %% 4)) %>%
          summarise(AvNumCharges = mean(NumCharges),.by=Age)

 

Since the p-value 0.11477 is greater than \(\alpha=0.05\) the evidence no longer suggests that heteroskedasticity is present in the model.

 

4) Data Summary

plot_start(df,geom_col())+ expand_limits(y = 3)+stat_poly_eq2(use_label("eq")) + 
            stat_poly_eq2( use_label("P"), label.y = .85, label.x = .05)

 

When arrest records are grouped by the arestees age into 4 year bins, and the mean number of charges per arrest is computed for each group: