72  Modello di Poisson

In questo capitolo imparerai a
  • utilizzare brms per adattare ai dati un modello di Poisson.
Prerequisiti
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)

72.1 Introduzione

In questo capitolo percorriamo l’intero processo che porta dall’idea di contare un evento raro – le sparatorie fatali da parte della polizia statunitense – alla stima del suo tasso medio annuo attraverso un modello di Poisson implementato con il pacchetto brms. Per contestualizzare il fenomeno si suggerisce di leggere l’articolo Racial Disparities in Police Use of Deadly Force Against Unarmed Individuals Persist After Appropriately Benchmarking Shooting Data on Violent Crime Rates (Ross et al., 2021). Non è obbligatorio per seguire il capitolo, ma fornisce lo sfondo sociale e metodologico dei dati che stiamo per analizzare.

72.2 Perché un modello di Poisson?

Quando il fenomeno di interesse è un conteggio – per esempio il numero di incidenti, diagnosi o, come nel nostro caso, sparatorie in un anno – la distribuzione di Poisson è spesso una scelta naturale. Questa distribuzione è definita da un solo parametro, \(\lambda\), che rappresenta la media (e varianza) del conteggio. In altre parole, se conosci \(\lambda\) conosci già la forma completa della distribuzione.

Per non appesantire la lettura ricordiamo qui solo l’essenziale: se \(Y\) è il numero di eventi osservati in un certo intervallo di tempo, dire che

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

significa che la probabilità di osservare esattamente \(k\) eventi è

\[ P(Y = k) = \frac{e^{-\lambda}\, \lambda^{k}}{k!}, \quad k = 0,1,2,\dots \]

72.3 Come brms codifica \(\lambda\)

Il pacchetto brms, così come la maggior parte dei software di regressione, non stima \(\lambda\) direttamente. Per garantire che il tasso rimanga positivo adotta un link logaritmico: all’interno del modello viene quindi stimata la quantità

\[ \eta = \log(\lambda). \]

Nel caso più semplice, senza predittori, brms utilizza un solo coefficiente, l’intercetta b_Intercept, che è proprio l’equivalente di \(\eta\). Una volta ottenuti i campioni posteriori di b_Intercept è sufficiente applicare l’esponenziale per tornare sul piano di \(\lambda\):

lambda <- exp(b_Intercept)  # trasforma log-lambda in lambda

Ogni volta che in questo capitolo nomineremo “campioni di \(\lambda\)” sottintenderemo che abbiamo già eseguito questa trasformazione.

72.4 La domanda di ricerca

Grazie all’archivio pubblico del Washington Post disponiamo di tutti i casi di sparatorie fatali accadute negli Stati Uniti dal 2015 in poi. L’interesse è stimare quante se ne verificano in media in un anno e descrivere l’incertezza associata a tale stima.

72.5 Importiamo e prepariamo i dati

url <- "https://raw.githubusercontent.com/washingtonpost/data-police-shootings/master/v2/fatal-police-shootings-data.csv"
raw <- read.csv(url, stringsAsFactors = FALSE)
raw$date <- as.Date(raw$date)
raw$year <- lubridate::year(raw$date)

# Escludiamo il 2025 perché l’anno è ancora in corso e i dati sarebbero incompleti
shootings <- subset(raw, year < 2025)

df <- shootings %>%
  dplyr::count(year, name = "events")
head(df)
#>   year events
#> 1 2015    995
#> 2 2016    959
#> 3 2017    984
#> 4 2018    992
#> 5 2019    993
#> 6 2020   1021

A questo punto abbiamo una tabella df con due colonne: year, che va dal 2015 al 2024, ed events, che contiene il numero di sparatorie registrate in ciascun anno.

72.6 Specificare la Distribuzione a Priori

Prima di osservare i dati vogliamo dichiarare che, secondo la nostra conoscenza precedente, un intervallo plausibile per \(\lambda\) va grosso modo da 400 a 900 casi l’anno, con media attorno a 600. Per ottenere una distribuzione lognormale con queste caratteristiche possiamo lavorare sulla scala logaritmica e scegliere normal(6.4, 0.3). Il valore 6.4 è infatti il logaritmo naturale di 600; la deviazione 0.3 produce l’ampiezza desiderata dell’intervallo. In brms la specifica è molto compatta:

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

72.7 Adattare il modello in brms

Il cuore dell’analisi si riduce a poche righe, perché il linguaggio di brms è pensato per assomigliare alla formula syntax di lm:

m0 <- brm(
  events ~ 1,                  # ~1 indica solo l’intercetta
  family = poisson(),          # distribuzione di errore
  data   = df,
  prior  = prior_lambda,
  iter   = 3000,
  warmup = 1000,
  chains = 4,
  seed   = 123,
  backend = "cmdstanr"
)

Il comando produce più di quattromila campioni posteriori (\(1000\) di warm‑up per ciascuna delle quattro catene + \(2000\) validi) dell’intercetta logaritmica. Le diagnosi di convergenza – \(\hat R\) vicino a 1, effective sample size buona – sono riportate da summary(m0).

72.8 Dalla scala log a \(\lambda\)

Per passare dalla stima di b_Intercept alla stima di \(\lambda\) basta eseguire l’esponenziale sui campioni posteriori. Usando tidybayes l’operazione è quasi in linguaggio naturale:

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

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  1043.  1023.  1062.   0.94 median qi

L’output ci dice che il valore più credibile di \(\lambda\) è intorno a 1043 casi l’anno, con un intervallo di credibilità al 94 % che va all’incirca da 1023 a 1062.

72.9 Visualizzare la distribuzione a posteriori

Un grafico spesso vale più di mille numeri. Con ggplot2 e il tema che abbiamo impostato all’inizio la figura è pronta in tre righe:

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à"
  )

L’area più scura al centro dell’half‑eye mette in evidenza l’intervallo più denso del 50 %; l’intera “pena” laterale dell’arco rappresenta invece il 94 %.

72.10 Un’estensione: confronto tra gruppi

Finora abbiamo trattato tutte le sparatorie fatali come se provenissero da un’unica popolazione. Una domanda molto concreta, invece, è se il tasso di vittime disarmate differisce tra persone classificate come bianche (codice W nel dataset) e tutte le altre. Per rispondere ci basta duplicare lo schema già usato: al posto di un solo tasso medio (λ) ne stimiamo due, uno per ciascun gruppo.

72.10.1 Costruire il dataset

Il primo passo è filtrare le righe che ci interessano (vittime disarmate) e poi contare, per ogni anno, quante di queste vittime appartengono a ciascun gruppo. Il codice qui sotto svolge tutto il lavoro di preparazione in un’unica catena di operazioni:

url <- "https://raw.githubusercontent.com/washingtonpost/data-police-shootings/master/v2/fatal-police-shootings-data.csv"

df_groups <- read.csv(url, stringsAsFactors = FALSE) |>
  dplyr::mutate(
    date  = as.Date(date),
    year  = lubridate::year(date),
    group = dplyr::if_else(race == "W", "White", "NonWhite")
  ) |>
  dplyr::filter(year < 2025, armed_with == "unarmed") |>
  dplyr::count(year, group, name = "events")

In df_groups ogni riga è l’abbinamento tra un anno e un gruppo, con il relativo conteggio di vittime disarmate.

72.10.2 Specificare il modello

La formula events ~ 0 + group dice a brms di rinunciare a un’intercetta comune e di stimarne una diversa per ogni valore di group. Poiché l’intercetta è, in scala logaritmica, il nostro parametro centrale, otteniamo due log‑tassi distinti.

Per la prior partiamo dall’idea che, in ciascun gruppo, potremmo aspettarci circa trenta vittime disarmate l’anno ma con incertezza ampia. Lavorando nella scala log questo si traduce in una distribuzione Normale con media 3.4 e deviazione 0.3 applicata indistintamente alle due intercette:

prior_race <- prior(normal(3.4, 0.3), class = "b")

Il modello completo è ora immediato:

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

72.10.3 Dal log‑tasso al tasso annuale

Dopo aver controllato che tutti gli indicatori di convergenza (in particolare R‑hat) siano a posto, trasformiamo i campioni posteriori con l’esponenziale così da tornare alla scala dei tassi veri e propri:

post <- as_draws_df(m_groups) |>
  dplyr::transmute(
    lambda_White    = exp(b_groupWhite),
    lambda_NonWhite = exp(b_groupNonWhite),
    diff_lambda     = lambda_NonWhite - lambda_White
  )

Con tidybayes riassumiamo in una riga le quantità d’interesse, usando ad esempio un intervallo di credibilità al 94 %:

post |>
  median_qi(lambda_White, lambda_NonWhite, diff_lambda, .width = 0.94)
#> # A tibble: 1 × 12
#>   lambda_White lambda_White.lower lambda_White.upper lambda_NonWhite
#>          <dbl>              <dbl>              <dbl>           <dbl>
#> 1         22.5               19.9               25.4            34.1
#>   lambda_NonWhite.lower lambda_NonWhite.upper diff_lambda diff_lambda.lower
#>                   <dbl>                 <dbl>       <dbl>             <dbl>
#> 1                  30.9                  37.6        11.6              7.21
#>   diff_lambda.upper .width .point .interval
#>               <dbl>  <dbl> <chr>  <chr>    
#> 1              16.0   0.94 median qi

Ne risulta che il gruppo NonWhite registra in media circa 12 vittime disarmate in più l’anno rispetto al gruppo White; l’intero intervallo di credibilità rimane sopra lo zero, perciò l’ipotesi di un tasso maggiore fra le persone non bianche è fortemente supportata dai dati.

Con due soli cambiamenti – la variabile group nella formula e una prior ragionevole per ciascun gruppo – abbiamo esteso il modello di Poisson a una comparazione fra categorie. Tutto il resto resta identico: link logaritmico, diagnostica delle catene, trasformazione a valle dei campioni. Una volta compreso questo meccanismo, aggiungere ulteriori gruppi o predittori diventa un esercizio di routine.

72.11 Riflessioni Conclusive

Il modello di Poisson con prior informativa e link logaritmico è un punto di partenza potente e relativamente semplice per modellare conteggi. L’implementazione in brms riduce la complessità sintattica al minimo, ma non per questo dobbiamo tralasciare di capire cosa accade dietro le quinte:

  • \(\lambda\) è il vero protagonista, ma viene stimato indirettamente tramite la sua trasformazione logaritmica.
  • La scelta della prior su b_Intercept va fatta nella scala log, pensando però a cosa significa nella scala originale.
  • Una volta compreso il passaggio dall’intercetta logaritmica al tasso \(\lambda\), è possibile arricchire il modello con più predittori, sapendo esattamente come trasformare il modello.

Esercizi

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).

  1. Considera che il numero di gol segnati da una squadra segua una Poisson con parametro \(\lambda\).
  2. 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).
  3. Aggiorna la distribuzione a posteriori conoscendo i gol segnati (5 per la Spagna e 3 per la Francia in una singola partita).
  4. 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.5
#> 
#> 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          data.table_1.17.2   
#> [22] labeling_0.4.3       bridgesampling_1.1-2 htmlwidgets_1.6.4   
#> [25] curl_6.2.2           pkgbuild_1.4.7       mnormt_2.1.1        
#> [28] RColorBrewer_1.1-3   cmdstanr_0.9.0       abind_1.4-8         
#> [31] withr_3.0.2          stats4_4.5.0         grid_4.5.0          
#> [34] colorspace_2.1-1     inline_0.3.21        xtable_1.8-4        
#> [37] emmeans_1.11.1       cli_3.6.5            mvtnorm_1.3-3       
#> [40] rmarkdown_2.29       generics_0.1.4       RcppParallel_5.1.10 
#> [43] rstudioapi_0.17.1    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] hms_1.1.3            arrayhelpers_1.1-0   ggdist_3.3.3        
#> [55] glue_1.8.0           codetools_0.2-20     ps_1.9.1            
#> [58] distributional_0.5.0 stringi_1.8.7        gtable_0.3.6        
#> [61] QuickJSR_1.7.0       pillar_1.10.2        htmltools_0.5.8.1   
#> [64] Brobdingnag_1.2-9    R6_2.6.1             rprojroot_2.0.4     
#> [67] evaluate_1.0.3       lattice_0.22-7       backports_1.5.0     
#> [70] rstantools_2.4.0     coda_0.19-4.1        nlme_3.1-168        
#> [73] checkmate_2.3.2      xfun_0.52            pkgconfig_2.0.3

Bibliografia

Downey, A. B. (2021). Think Bayes. " O’Reilly Media, Inc.".
Ross, C. T., Winterhalder, B., & McElreath, R. (2021). Racial disparities in police use of deadly force against unarmed individuals persist after appropriately benchmarking shooting data on violent crime rates. Social Psychological and Personality Science, 12(3), 323–332.