---
title: "Analyzing the Stop, Question, and Frisk Program: A Statistical Approach to Equal Protection"
author: "by Mauricio Claudio"
output: 
  flexdashboard::flex_dashboard:
    theme: simplex
    orientation: columns
    vertical_layout: fill
    source_code: embed
---

---
title: "Analyzing the Stop, Question, and Frisk Program: A Statistical Approach to Equal Protection"
author: ''
output: flexdashboard::flex_dashboard
---



```{r setup, include=FALSE}
library(ggstream)
library(sjPlot)
library(DT)
library(flexdashboard)
library(vembedr)
library(glmtoolbox)
library(knitr)
library(plotly)
library(performance)
library(summarytools)
library(tidyverse)
library(emmeans)
library(parcats)
library(easyalluvial)

raw.dat = read.csv("sqf-2023.csv")
```


```{r message=FALSE, warning=FALSE}
NYPDfrisk2018 = raw.dat |>
   dplyr::select(SUSPECT_REPORTED_AGE,
                 FRISKED_FLAG,
                 SUSPECT_SEX,
                 SUSPECT_RACE_DESCRIPTION,
                 STOP_LOCATION_BORO_NAME) |>
   filter(SUSPECT_RACE_DESCRIPTION!='(null)',
          SUSPECT_SEX!='(null)',
          SUSPECT_RACE_DESCRIPTION !='MIDDLE EASTERN/SOUTHWEST ASIAN',
          SUSPECT_RACE_DESCRIPTION !='AMERICAN INDIAN/ALASKAN NATIVE',
          SUSPECT_RACE_DESCRIPTION !='ASIAN / PACIFIC ISLANDER',
          STOP_LOCATION_BORO_NAME != 'STATEN ISLAND') |>
   mutate(frisked = factor(ifelse(FRISKED_FLAG=='Y','yes','no')),
          boro = factor(STOP_LOCATION_BORO_NAME,
                        c('MANHATTAN','BRONX','BROOKLYN','QUEENS')),
          age = as.integer(SUSPECT_REPORTED_AGE),
          sex = factor(SUSPECT_SEX),
          race = factor(SUSPECT_RACE_DESCRIPTION,
                        c('WHITE', 'WHITE HISPANIC', 'BLACK HISPANIC', 'BLACK'))) |>
   dplyr::select(race, sex, boro, age, frisked) |>
   filter(age <=60,
          age >=16) |>
   drop_na()
```



**Introduction**
=====================================  


Column {data-width=400}
-------------------------------------

### {data-padding=8 data-height=400}

#### **Summary**

In 2022, New York City Mayor E. Adams reinstated the controversial Stop, Question and Frisk program, allowing police officers to stop individuals based on reasonable suspicion and conduct body searches if necessary. The program has faced criticism, particularly after a 2013 court ruling deemed it discriminatory against minorities. Despite this, it remains in effect. This statistical analysis of the program reveals significant disparities in frisk probabilities across different race, sex, age and borough groups. It highlights the need for robust data and nuanced analysis to better understand its impact on equal protection, offering a more informed approach to the ongoing debate.

___

### {data-padding=8}


```{r}
embed_url('https://youtu.be/2_hpXMExtJ4?si=HdzelVEMgzTH5BbN') |>
   use_bs_responsive()
```



Column {data-width=600}
-------------------------------------

### {data-padding=8}

#### **Introduction**  

In 2022, shortly after taking office, New York City Mayor Eric Adams reinstated the New York Police Department's (NYPD) Stop, Question, and Frisk program. This policy allows officers to temporarily stop individuals based on reasonable suspicion. If the officer further believes the individual is armed and dangerous, a body pat-down, or frisk, is permitted. However, the program has a controversial history. In 2013, a U.S. District Court ruled that the NYPD's implementation violated equal protection constitutional provisions due to its disproportionate impact on minorities.

Despite ongoing criticism, the policy remains in use and continues to fuel debate. At the core of this debate is the issue of equal protection. Yet, advocacy efforts and media coverage often rely on inconsistent estimates of frisk occurrences. Discussions are frequently hindered by vague references to victims as "persons of color" and a lack of robust, comparative data that would allow for meaningful inferences about the affected population.

To address this gap, we conduct a statistical analysis to provide more rigorous and nuanced estimates of how the program is applied across 32 distinct race | sex | borough cohorts. Our analysis seeks to answer two key questions:  

1.	What is the probability of a frisk?

2.	How do the odds of a frisk compare to those of a white subject in Manhattan?
  
Unlike simple descriptive calculations of averages and counts, our approach estimates the relative differences and uncertainties in frisk probabilities. It also accounts for variations not only by race but also by sex, age, and borough. This method allows for more robust inferences about the broader New York City population and directly addresses concerns about equal protection.

___

Researched and engineered at **Decision Analytics Hub**  
by [Mauricio Claudio](https://www.linkedin.com/in/mauricio-claudio-550a6944/){target="_blank"} | <DecisionAnalyticsHub@proton.me>





**Data**
=====================================


Column {data-width=380}
-----------------------------------------------------------------------

### {data-padding=8}

#### **Data**

The 2023 [Stop, Question, and Frisk data](https://www.nyc.gov/site/nypd/stats/reports-analysis/stopfrisk.page){target="_blank"} used in this analysis comes from the New York Police Department website. It includes 12,035 recorded instances of police stops, each documenting key characteristics of the individual and the stop location. Each instance consists of:  
  
* Three categorical variables: *race*, *sex*, and *boro* (borough).  

* One numerical variable: *age* (the individual's age).  

* One binary outcome variable: *frisked* (indicating whether a frisk occurred).  

To maintain statistical reliability, observations for the *race* categories MIDDLE EASTERN/SOUTHWEST ASIAN, AMERICAN INDIAN/ALASKAN NATIVE, and ASIAN/PACIFIC ISLANDER, as well as stops in *boro* STATEN ISLAND, were excluded due to their small sample sizes. Additionally, only individuals aged 16 to 60 were included in the analysis.  

The graphs are interactive - hover over them for greater detail on the data!

___


Column {data-width=620 .tabset}
-----------------------------------------------------------------------

### Summary

#### Summary table
  
```{r}
### Print data summary
print(dfSummary(NYPDfrisk2018,
                graph.col=TRUE,
                graph.magnif = 1,
                plain.ascii = FALSE,
                headings = FALSE,
                display.labels = FALSE,
                silent = TRUE,
                valid.col = TRUE,
                na.col = TRUE),
      method = "render")
```
  
___
  
  
#### Grouped totals
  
```{r}
ggdat = NYPDfrisk2018 |>
  group_by(race, sex, boro, frisked) |>
  count()

ggdat |>
  mutate(PercentTotal = ( (n / nrow(NYPDfrisk2018)) *100),
          PercentTotal = round(PercentTotal, 2)) |>
  arrange(desc(frisked), desc(n)) |>
  rename(Borough=boro, Sex=sex, Race=race, Frisk=frisked, Total=n) |>
  dplyr::select(Race, Sex, Borough, Frisk, Total, PercentTotal) |>
  kable(digits=5)

```



### Stops

```{r}
library(dplyr)
library(highcharter)

hcdat = NYPDfrisk2018 |> group_by(race, sex, boro) |> count()
hc = hchart(hcdat, "packedbubble", hcaes(name = race:sex, value = n, group = boro))

hc %>%
  hc_title(
    text = "Police Stops by Borough, Race and Sex",
    margin = 20,
    align = "center",
    style = list(color = "#000000", useHTML = TRUE)
    ) |>
  hc_tooltip(
    useHTML = TRUE,
    pointFormat = "<b>{point.name}:</b> {point.value}"
  ) %>%
  hc_plotOptions(
    packedbubble = list(
      maxSize = "150%",
      zMin = 0,
      layoutAlgorithm = list(
        gravitationalConstant =  0.05,
        splitSeries =  TRUE, # TRUE to group points
        seriesInteraction = TRUE,
        dragBetweenSeries = FALSE,
        parentNodeLimit = TRUE
      ),
      dataLabels = list(
        enabled = TRUE,
        format = "{point.name}",
        filter = list(
          property = "y",
          operator = ">",
          value = 5
        ),
        style = list(
          color = "black",
          textOutline = "none",
          fontWeight = "normal"
        )
      )
    )
  )
```



### Stop outcome

```{r}
ggplotly(
   ggplot(ggdat, aes(y=reorder(race:boro, n), x=n, fill=frisked)) +
   geom_col() + theme_minimal(base_size=10) +
   facet_grid(~ factor(sex, levels = c('MALE','FEMALE')), scales='free_x') +
      scale_fill_manual(values=c('darkgreen','red3')) +
   labs(title='Stop outcome by Sex, Race & Borough',
        y='', x='Total', fill='Frisked') +
     theme(strip.text=element_text(face="bold", size=10))
)
```



### Stop outcome, rel.

```{r}
ggplotly(
   ggplot(ggdat, aes(y=reorder(race:boro, n), x=n, fill=frisked)) +
   geom_col(position='fill') + theme_minimal(base_size=10) +
   facet_grid(~ factor(sex, levels = c('MALE','FEMALE')), scales='free_x') +
      scale_fill_manual(values=c('darkgreen','red3')) +
   labs(title='Stop outcome by Sex, Race & Borough, relative',
        y='', x='Total', fill='Frisked')
)
```



### Stop outcome: Age

```{r}
ggplotly(
NYPDfrisk2018 |> group_by(frisked, age, race) |>
   rename(Frisked=frisked) |> count() |>
   ggplot(aes(x=age, y=n, fill=race)) +
   geom_col() + theme_minimal() +
   facet_wrap(~Frisked, scales='free', labeller='label_both') +
   scale_fill_manual(values = c("khaki3", "orange2", 'red2','red4')) +
   labs(title='Stop outcome by Age', x='Age', y='Total', fill='Race') +
   theme(strip.text=element_text(face="bold", color="darkgreen", hjust=0, size=10))
)
```



### Stop outcome: Age, rel.

```{r}
ggplotly(
NYPDfrisk2018 |> group_by(frisked, age, race) |>
   rename(Frisked=frisked) |> count() |>
   ggplot(aes(x=age, y=n, fill=race)) +
   geom_col(position="fill") +
   theme_minimal() +
   facet_grid(~Frisked, scales='free', labeller='label_both') +
   scale_fill_manual(values = c("khaki3", "orange2", 'red2','red4')) +
   labs(title='Stop outcome by Age, relative', x='Age', y='Total', fill='Race') +
   theme(strip.text=element_text(face="bold", color="darkgreen", hjust=0, size=10))
)
```



### Alluvial

```{r fig.height=5, fig.width=9, message=FALSE, warning=FALSE}
p <- alluvial_wide(NYPDfrisk2018 |> dplyr::select(frisked, race, sex, boro, age),
                   max_variables = 5, bin_labels = 'min_max',
                   fill_by = 'first_variable',
                   col_vector_flow = c('#01700c', '#cf0202'))

parcats(p, marginal_histograms = FALSE, data_input = NYPDfrisk2018,
        hoveron = 'category', hoverinfo = 'count', arrangement='freeform',
        sortpaths = 'forward', labelfont = list(size=24), tickfont = list(size=12))
```






```{r eval=FALSE, include=FALSE}
### Download
NYPDfrisk2018 |>
   datatable(extensions = c('Buttons','ColReorder'),
             options = list(pageLength = 100,
                            dom = 'Blrtip',
                            buttons = c('csv'),
                            colReorder = TRUE,
                            deferRender = TRUE))
```





**Method**
=====================================


Column {data-width=380}
-----------------------------------------------------------------------

### {data-padding=8}

#### **Model**

The data is analyzed using a binary logistic regression model consisting of four interaction terms: *race:boro*, *boro:sex*, *race:age*, and *boro:age*. This model was selected over alternative models based on its lower AIC and BIC scores, as well as its superior performance in posterior prediction and the Hosmer-Lemeshow goodness-of-fit test. Deviance analysis and joint tests further validate the model's fit and confirm the significance of these interaction terms.

Mean probabilities and odds-ratio contrasts are estimated with 95% confidence at four-year intervals of the age variable. The contrasts are adjusted using the Dunnett method.

The analysis is conducted in `R` (v4.1.0) and `RStudio` (v2023.12.1+402). Mean probabilities and odds-ratios are estimated with the `emmeans` library (v1.10.1). 


___


Column {.tabset data-width=620}
-----------------------------------------------------------------------

### Summary

#### Model

$$\mathbf{Y = \alpha\beta + \beta\gamma + \alpha\delta + \beta\delta}$$

$$where\ Y = frisked
\\ {\bf \alpha} = race
\\ {\bf \beta} = boro
\\ {\bf \gamma} = sex
\\ {\bf \delta} = age$$

___



#### Deviance
```{r echo=FALSE, message=FALSE, warning=FALSE}
mod1= glm(frisked ~ race:boro + race:age + sex:boro + boro:age,
           data = NYPDfrisk2018, family=binomial(link="logit"))

car::Anova(mod1)
```

___

  
#### Joint tests

```{r message=FALSE, warning=FALSE}
joint_tests(mod1)
```

___


### Diagnostics


#### Posterior Predictions
```{r fig.align='center', fig.height=5, fig.width=9}
check_predictions(mod1)
```
  
___  
___


#### Hosmer-Lemeshow test
```{r fig.align='center'}
hltest(mod1)
```



####

```{r fig.align='center'}
###
hl = data.frame(hltest(mod1, verbose=FALSE))

ggplotly(
   ggplot(hl, aes(x=hm.Observed, y=hm.Expected)) +
      geom_point(size=6, color='red3') + theme_minimal(base_size = 9) +
      geom_abline(slope = 1, color='black', alpha=0.2, size=2) +
      labs(color='Bins', x='Observed', y='Expected',
           title='Hosmer-Lemeshow goodness-of-fit test')
)
```
  
___





**Results**
=====================================


Column {data-width=380}
-----------------------------------------------------------------------

### {data-padding=8}

#### **Results**

The results are displayed on the plots and tables on the right. In the plots, the point is the estimate and the line is the 95% confidence interval around the estimate. Probabilities indistinct from random chance and relative odds at parity appear as points with muted color. ``Play`` cycles through *age* or it can be selected on the slider.   
  
**Frisk probability**  
The probability of a frisk shows great variability among the 32 *race | sex | borough* cohorts and is greatly dependent on age. At the lower age range of 16 years, the range of probabilities ranges from a high of 0.92 for a black male in the Bronx to a low of 0.16 for a white female in Brooklyn. The probability of six cohorts is indistinct from the random chance probability of 0.5.

At the mean age of 28 years, frisk probabilities continue to display great variability among cohorts, with a maximum of 0.85 for black male in the Bronx and a minimum of 0.14 for white females in Brooklyn, while the probabilities of nine cohorts are indistinct from random chance.  
  
At the upper age range of 60 years, frisk probabilities range from 0.62 for a white male in the Bronx to 0.10 for a white female in Brooklyn, with the probability of four highest cohorts indistinct from random chance. Generally, individual cohort probabilities and probability differences among cohorts both decrease with increasing age.
  
  
**Relative odds**  
The relative of odds of frisk among individual *race | borough* cohorts against the *white | Manhattan* reference cohort displays substantial variability and is highly dependent on age. At the lower age range of 16 years, a black in the Bronx has 14 times greater odds of a frisk than the reference cohort while the odds of a white Hispanic are 2.4 times greater. 

At the mean age of 28 years, relative odds display substantial variability among cohorts, with blacks in the Bronx having 8.7 times and white Hispanics 8.4 greater odds of a frisk than the reference cohort. By contrast, white Hispanics in Queens and Brooklyn have equal odds to the reference cohort.   
  
At the upper age range of 60 years, a white in the Bronx has 3.7 times greater odds than the reference cohorts while a white Hispanic in the same borough has 3.2 times greater odds. All other relative odds at this age are indistinct from parity. Generally, the individual cohort odds of a frisk and the odds differences among cohorts both decrease with increasing age.  
  
___


Column {data-width=620 .tabset}
-----------------------------------------------------------------------

### Frisk probability


```{r message=FALSE, warning=FALSE}
###
age.dat = data.frame(emmeans(object=mod1,
        specs = ~race:sex:boro,
        by='age', at = list(age = seq(16,60,4)),
        level=0.95, type='response')) |> dplyr::select(-SE, -df) |>
   mutate(ctr = factor(paste(race, sex, 'in', boro)),
          chance = ifelse(asymp.LCL > 0.5 | asymp.UCL < 0.5, 1, 0.5))

###
ggplotly(
   ggplot(age.dat, aes(x=prob, y=reorder(ctr, prob), frame=age)) +
      geom_pointrange(aes(xmin=asymp.LCL, xmax=asymp.UCL,
                          color=sex, alpha=chance), size=2) +
      labs(title='Frisk probability by Race, Sex, Borough and Age',
           y='', x='Probability', alpha='') +
      geom_vline(xintercept = 0.50, linetype='dotted') +
      theme_minimal(base_size = 9) + scale_x_continuous(breaks=seq(0,1,0.1)) +
      scale_color_manual(values = c("FEMALE"="red4", "MALE" = 'dodgerblue4'))
) |> animation_slider(currentvalue=list(prefix="Age: ", font=list(color="darkred")))
```




  
___
  
  
####

```{r echo=FALSE, warning=FALSE}
age.dat |> dplyr::select(age, race, sex, boro, prob, asymp.LCL, asymp.UCL) |>
   arrange(age, desc(prob)) |>
   kable(digits=3)
```
  
___



### Relative odds

```{r message=FALSE, warning=FALSE}
odds.dat = emmeans(mod1, "trt.vs.ctrl" ~ race:boro,
             by='age', at = list(age = seq(16,60,4)),
             level=0.95, type='response', infer=TRUE)

odds.dat = data.frame(odds.dat$contrasts) |>
   dplyr::select(-SE, -df, -null, -z.ratio) |>
   mutate(age=as.integer(as.character(age)),
          significant=ifelse(p.value<=0.05, 1, 0))

#glimpse(odds.dat)

ggplotly(
   ggplot(odds.dat,
          aes(x=odds.ratio, y=reorder(contrast, odds.ratio), frame=age,
              alpha=significant)) +
      geom_pointrange(aes(xmin=asymp.LCL, xmax=asymp.UCL),
                      size=4, color='red4') +
      labs(title='Relative odds by Race, Borough and Age',
           y='', x='Odds') +
      geom_vline(xintercept = 1, linetype='dotted') +
      theme_minimal(base_size = 10) +
      scale_x_continuous(breaks=seq(0,30,10))
) |> animation_slider(currentvalue=list(prefix="Age: ", font=list(color="darkred")))
```
  
___
  
  
#### 

```{r}
### 5.
odds.dat |>
  arrange(age, desc(odds.ratio)) |>
  dplyr::select(-significant) |>
  kable(digit=3) #|>
  # datatable(extensions = c('Buttons'),
  #           options = list(pageLength = 15,
  #                           dom = 'Blrtip',
  #                           buttons = c('copy', 'pdf', 'csv'))) |>
  #  formatRound(columns=c('odds.ratio',
  #                        'asymp.LCL',
  #                        'asymp.UCL',
  #                        'p.value'), digits=5)
```






**Conclusion**
=====================================


Column {data-width=380}
-----------------------------------------------------------------------

### {data-padding=8}


#### **Conclusion**
  
<br>  

*"All animals are equal, but some animals are more equal that others."*  
  
<br>

The analysis uncovers significant disparities in frisk probabilities and odds when compared to the white-in-Manhattan reference group, with variations across race, sex, borough, and age. Additionally, it highlights substantial uncertainties in the data-an aspect often overlooked in public discourse. By incorporating both these estimates and their uncertainties, the public can take a more informed and measured approach in evaluating whether and to what extent Stop, Question, and Frisk deviates from the principles of equal protection. 
  
___





Column {data-width=620}
-----------------------------------------------------------------------

### 
![](img1.jpeg)