73 Modello di Poisson

In questo capitolo imparerai a

utilizzare brms per adattare ai dati un modello di Poisson.

Prerequisiti

Leggere Racial Disparities in Police Use of Deadly Force Against Unarmed Individuals Persist After Appropriately Benchmarking Shooting Data on Violent Crime Rates per ottenere una panoramica approfondita su questo fenomeno e sul relativo ambito di ricerca.

Preparazione del Notebook

here::here("code", "_common.R") |> 
  source()

# Load packages
if (!requireNamespace("pacman")) install.packages("pacman")
pacman::p_load(HDInterval, lubridate, brms, bayesplot, tidybayes)

73.1 Introduzione

Nel capitolo precedente abbiamo visto come, in un contesto bayesiano, la distribuzione a posteriori del parametro di una distribuzione di Poisson possa essere calcolata usando una prior Gamma. In questo capitolo, mostreremo come eseguire un’analisi simile in R usando il pacchetto brms, focalizzandoci sull’interpretazione diretta del tasso \(\lambda\) di sparatorie fatali per anno negli Stati Uniti.

73.2 Domanda della ricerca

Analizziamo i dati raccolti dal Washington Post su tutte le sparatorie fatali da parte della polizia statunitense dal 2015. Il nostro obiettivo è stimare il tasso medio annuo \(\lambda\) di queste sparatorie e quantificare l’incertezza associata.

73.3 Preparazione dei dati

Importiamo i dati direttamente da GitHub:

url <- "https://raw.githubusercontent.com/washingtonpost/data-police-shootings/master/v2/fatal-police-shootings-data.csv"
fps_dat <- read.csv(url, stringsAsFactors = FALSE)
fps_dat$date <- as.Date(fps_dat$date)
fps_dat$year <- as.numeric(format(fps_dat$date, "%Y"))
fps <- subset(fps_dat, year != 2025)  # Rimuove dati parziali del 2025

Creiamo un dataset aggregato per anno:

year_counts <- table(fps$year)
df <- data.frame(
  year   = as.numeric(names(year_counts)),
  events = as.vector(year_counts)
)

73.4 Specifica del modello

Assumiamo che il numero di sparatorie \(Y\) in ciascun anno segua una distribuzione di Poisson con tasso costante \(\lambda\):

\[ Y \sim \text{Poisson}(\lambda) \]

Vogliamo stimare \(\lambda\), che rappresenta il numero medio di eventi per anno.

73.5 Scelta della distribuzione a priori

Invece di ragionare in termini di \(\log(\lambda)\), possiamo scegliere una distribuzione lognormale che, trasformata, abbia una media plausibile per \(\lambda\). Supponiamo di credere, a priori, che il numero medio annuo di sparatorie sia circa 600, con incertezza di circa ±200.

Una distribuzione lognormale approssima bene questa idea. Una normal(6.4, 0.3) sul logaritmo corrisponde a una distribuzione su \(\lambda\) centrata intorno a 600, con valori plausibili tra circa 400 e 900.

Scriviamo:

prior_lambda <- c(
  prior(normal(6.4, 0.3), class = "Intercept")
)

Anche se specifichiamo la prior su \(\log(\lambda)\) (per compatibilità con brms), interpreteremo i risultati direttamente sulla scala di \(\lambda\).

73.6 Stima di \(\lambda\) con `brms`

Costruiamo un modello molto semplice, senza effetti per anno:

m0 <- brm(
  formula = events ~ 1,
  family  = poisson(),
  data    = df,
  prior   = prior_lambda,
  iter    = 3000, warmup = 1000, chains = 4, seed = 123
)

Otteniamo un riepilogo:

summary(m0)
#>  Family: poisson 
#>   Links: mu = log 
#> Formula: events ~ 1 
#>    Data: df (Number of observations: 10) 
#>   Draws: 4 chains, each with iter = 3000; warmup = 1000; thin = 1;
#>          total post-warmup draws = 8000
#> 
#> Regression Coefficients:
#>           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept     6.95      0.01     6.93     6.97 1.00     3290     3457
#> 
#> Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).

73.7 Interpretazione dei risultati

Per ottenere direttamente i valori del tasso \(\lambda\) (non il log), possiamo trasformare i campioni posteriori dell’intercetta:

posterior_lambda <- as_draws(m0) %>%
  spread_draws(b_Intercept) %>%
  mutate(lambda = exp(b_Intercept))

Otteniamo la stima e l’intervallo di credibilità al 94%:

posterior_lambda %>%
  median_qi(lambda, .width = 0.94)
#> # A tibble: 1 × 6
#>   lambda .lower .upper .width .point .interval
#>    <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
#> 1  1042.  1023.  1062.   0.94 median qi

73.8 Visualizzazione

posterior_lambda %>%
  ggplot(aes(x = lambda)) +
  stat_halfeye(fill = "skyblue", alpha = 0.6) +
  labs(
    title = "Distribuzione a posteriori del tasso λ",
    x = "Tasso annuo di sparatorie fatali (λ)",
    y = "Densità"
  )

73.9 Estensione: confronto tra gruppi

Supponiamo di voler confrontare il tasso di sparatorie fatali di vittime disarmate tra due gruppi: persone bianche e non bianche.

73.9.1 Filtraggio e aggregazione

fps_unarmed <- subset(fps, armed_with == "unarmed")
white_df <- subset(fps_unarmed, race == "W")
non_white_df <- subset(fps_unarmed, race != "W")

count_white <- table(white_df$year)
count_non_white <- table(non_white_df$year)

df_white <- data.frame(year = as.numeric(names(count_white)), events = as.vector(count_white), group = "White")
df_nonwhite <- data.frame(year = as.numeric(names(count_non_white)), events = as.vector(count_non_white), group = "NonWhite")
df_groups <- rbind(df_white, df_nonwhite)

73.9.2 Prior plausibile

Supponiamo una media a priori di circa 30 eventi l’anno con deviazione 10. Una lognormale approssimativa corrisponde a normal(3.4, 0.3) sulla scala logaritmica.

prior_race <- c(
  prior(normal(3.4, 0.3), class = "b", coef = "groupWhite"),
  prior(normal(3.4, 0.3), class = "b", coef = "groupNonWhite")
)

73.9.3 Modello con due gruppi

m_groups <- brm(
  formula = events ~ 0 + group,
  family  = poisson(),
  data    = df_groups,
  prior   = prior_race,
  iter = 3000, warmup = 1000, chains = 4, seed = 123
)

Otteniamo i tassi stimati:

post <- posterior_samples(m_groups)
lambda_white <- exp(post$b_groupWhite)
lambda_nonwhite <- exp(post$b_groupNonWhite)
diff_lambda <- lambda_nonwhite - lambda_white

Intervallo credibile:

median_qi(diff_lambda, .width = 0.94)
#>       y  ymin  ymax .width .point .interval
#> 1 11.69 7.245 16.12   0.94 median        qi

73.10 Conclusioni

Il tasso medio annuo stimato di vittime disarmate uccise dalla polizia è più alto per il gruppo non caucasico rispetto a quello caucasico.
La differenza media stimata è di circa 11.5 casi l’anno, con un intervallo di credibilità al 94% che va da circa 6.8 a 16.4.
I tassi stimati per ciascun gruppo sono:
- Non caucasici: media ≈ 35.6 (CI: 32–39.3)
- Caucasici: media ≈ 24.1 (CI: 21.1–27.3)

Questi risultati forniscono una base solida per affermare che il tasso di sparatorie fatali di vittime disarmate è più alto tra le persone non caucasiche.

Esercizi

Problemi

Nella finale olimpica di calcio 2024, la Spagna ha sconfitto la Francia per 5 a 3. Supponiamo di voler calcolare la probabilità di superiorità della Spagna rispetto alla Francia utilizzando un modello coniugato Gamma-Poisson (o l’approssimazione brms con prior lognormale).

Considera che il numero di gol segnati da una squadra segua una Poisson con parametro \(\lambda\).
Specifica un prior su \(\lambda\) per entrambe le squadre, ad esempio \(\alpha=1\) e \(\beta=1\) nella parametrizzazione Gamma classica (oppure una Normal(0,1.4) sull’intercetta, in modo da avere una media a posteriori analoga).
Aggiorna la distribuzione a posteriori conoscendo i gol segnati (5 per la Spagna e 3 per la Francia in una singola partita).
Calcola la probabilità che \(\lambda_{\text{Spagna}} > \lambda_{\text{Francia}}\).

(Ispirato a “The World Cup Problem”, (Downey, 2021).)

Suggerimento: puoi risolvere il problema in modo analitico (Gamma-Poisson con un solo conteggio) oppure puoi usare brms costruendo un dataframe:

df_soccer <- data.frame(
  team = c("Spain", "France"),
  goals = c(5, 3)
)

Modello: goals ~ 0 + team, family=poisson().
Prior su b_teamSpain e b_teamFrance.
Infine, estrai i draws e calcola la probabilità \(\Pr(\exp(b_{\text{Spain}}) > \exp(b_{\text{France}}))\).

Informazioni sull’Ambiente di Sviluppo

sessionInfo()
#> R version 4.5.0 (2025-04-11)
#> Platform: aarch64-apple-darwin20
#> Running under: macOS Sequoia 15.4.1
#> 
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib 
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
#> 
#> locale:
#> [1] C/UTF-8/C/C/C/C
#> 
#> time zone: Europe/Rome
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#>  [1] tidybayes_3.0.7  brms_2.22.0      Rcpp_1.0.14      HDInterval_0.2.4
#>  [5] thematic_0.1.6   MetBrewer_0.2.0  ggokabeito_0.1.0 see_0.11.0      
#>  [9] gridExtra_2.3    patchwork_1.3.0  bayesplot_1.12.0 psych_2.5.3     
#> [13] scales_1.4.0     markdown_2.0     knitr_1.50       lubridate_1.9.4 
#> [17] forcats_1.0.0    stringr_1.5.1    dplyr_1.1.4      purrr_1.0.4     
#> [21] readr_2.1.5      tidyr_1.3.1      tibble_3.2.1     ggplot2_3.5.2   
#> [25] tidyverse_2.0.0  rio_1.2.3        here_1.0.1      
#> 
#> loaded via a namespace (and not attached):
#>  [1] tidyselect_1.2.1     svUnit_1.0.6         farver_2.1.2        
#>  [4] loo_2.8.0            fastmap_1.2.0        tensorA_0.36.2.1    
#>  [7] pacman_0.5.1         digest_0.6.37        timechange_0.3.0    
#> [10] estimability_1.5.1   lifecycle_1.0.4      StanHeaders_2.32.10 
#> [13] processx_3.8.6       magrittr_2.0.3       posterior_1.6.1     
#> [16] compiler_4.5.0       rlang_1.1.6          tools_4.5.0         
#> [19] utf8_1.2.5           yaml_2.3.10          labeling_0.4.3      
#> [22] bridgesampling_1.1-2 htmlwidgets_1.6.4    curl_6.2.2          
#> [25] pkgbuild_1.4.7       mnormt_2.1.1         plyr_1.8.9          
#> [28] RColorBrewer_1.1-3   abind_1.4-8          withr_3.0.2         
#> [31] grid_4.5.0           stats4_4.5.0         colorspace_2.1-1    
#> [34] inline_0.3.21        xtable_1.8-4         emmeans_1.11.1      
#> [37] cli_3.6.5            mvtnorm_1.3-3        rmarkdown_2.29      
#> [40] generics_0.1.3       RcppParallel_5.1.10  rstudioapi_0.17.1   
#> [43] reshape2_1.4.4       tzdb_0.5.0           rstan_2.32.7        
#> [46] parallel_4.5.0       matrixStats_1.5.0    vctrs_0.6.5         
#> [49] V8_6.0.3             Matrix_1.7-3         jsonlite_2.0.0      
#> [52] callr_3.7.6          hms_1.1.3            arrayhelpers_1.1-0  
#> [55] ggdist_3.3.3         glue_1.8.0           ps_1.9.1            
#> [58] codetools_0.2-20     distributional_0.5.0 stringi_1.8.7       
#> [61] gtable_0.3.6         QuickJSR_1.7.0       pillar_1.10.2       
#> [64] htmltools_0.5.8.1    Brobdingnag_1.2-9    R6_2.6.1            
#> [67] rprojroot_2.0.4      evaluate_1.0.3       lattice_0.22-7      
#> [70] backports_1.5.0      rstantools_2.4.0     coda_0.19-4.1       
#> [73] nlme_3.1-168         checkmate_2.3.2      xfun_0.52           
#> [76] pkgconfig_2.0.3

Bibliografia

Downey, A. B. (2021). Think Bayes. " O’Reilly Media, Inc.".

73.1 Introduzione

73.2 Domanda della ricerca

73.3 Preparazione dei dati

73.4 Specifica del modello

73.5 Scelta della distribuzione a priori

73.6 Stima di \(\lambda\) con brms

73.7 Interpretazione dei risultati

73.8 Visualizzazione

73.9 Estensione: confronto tra gruppi

73.9.1 Filtraggio e aggregazione

73.9.2 Prior plausibile

73.9.3 Modello con due gruppi

73.10 Conclusioni

Esercizi

Informazioni sull’Ambiente di Sviluppo

Bibliografia

73.6 Stima di \(\lambda\) con `brms`