72  Confronto tra due proporzioni indipendenti

In questo capitolo imparerai a

• condurre un confronto bayesiano tra le proporzioni di due gruppi indipendenti utilizzando la funzione brm() del pacchetto brms.

Prerequisiti

• Leggere l’articolo “Children’s arithmetic skills do not transfer between applied and academic mathematics” (Banerjee et al., 2025).

Preparazione del Notebook

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

# Load packages
if (!requireNamespace("pacman")) install.packages("pacman")
pacman::p_load(brms, cmdstanr, posterior, brms, bayestestR, insight)

72.1 Introduzione

Supponiamo di voler capire se due gruppi di persone hanno la stessa probabilità di “successo” in una certa attività. Per esempio, vogliamo sapere se due trattamenti diversi portano alla stessa percentuale di guarigione, oppure se studenti che seguono due metodi di studio differenti superano un esame con la stessa frequenza.

Quando il risultato per ogni persona è un valore binario — successo o insuccesso, sì o no, guarito o non guarito — possiamo usare un modello statistico chiamato regressione logistica.

In particolare, se i due gruppi sono indipendenti e distinti (cioè ogni persona appartiene a uno solo dei due gruppi), possiamo usare una versione semplice della regressione logistica con una sola variabile esplicativa binaria (una “dummy”).

Adottando un approccio bayesiano, possiamo:

1. rendere esplicita l’incertezza sulle nostre ipotesi iniziali tramite le distribuzioni a priori;
2. descrivere l’intera gamma di risultati plausibili, ottenendo una distribuzione a posteriori dei parametri invece di un singolo valore “stimato”.

72.2 La struttura dei dati

Consideriamo i dati che raccogliamo:

• ogni partecipante è identificato da un indice \(i = 1, 2, \dots, N\);

• per ciascuno osserviamo un esito binario:

  \[ y_i = \begin{cases} 1 & \text{se c’è un successo},\\ 0 & \text{se c’è un insuccesso}. \end{cases} \]

• ogni partecipante appartiene a uno e un solo gruppo. Chiamiamo gruppo 0 il primo gruppo (ad esempio, il gruppo di controllo) e gruppo 1 il secondo gruppo (ad esempio, il gruppo che riceve un trattamento sperimentale).

Per rappresentare questa appartenenza al gruppo, definiamo una variabile indicatrice:

\[ D_i = \begin{cases} 0 & \text{se il partecipante è nel gruppo 0},\\ 1 & \text{se il partecipante è nel gruppo 1}. \end{cases} \]

Questa variabile ci permette di costruire un modello che distingue i due gruppi.

72.3 Un modello statistico per dati binari

Poiché il risultato \(y\_i\) può essere solo 0 o 1, possiamo descriverlo con una distribuzione di Bernoulli, che rappresenta proprio questo tipo di variabili:

\[ y_i \sim \text{Bernoulli}(p_i), \]

dove \(p\_i\) è la probabilità che il partecipante \(i\) ottenga un successo (cioè che \(y_i = 1\)).

A questo punto potremmo pensare di modellare direttamente \(p_i\) come una funzione di \(D_i\). Ma c’è un problema tecnico importante.

72.3.1 Perché non modelliamo direttamente la probabilità \(p_i\)?

Le probabilità devono sempre stare tra 0 e 1. Ma se usassimo un modello lineare classico (come \(p\_i = \alpha + \gamma D_i\)), potremmo ottenere dei valori fuori da questo intervallo — ad esempio, una “probabilità” negativa o maggiore di 1, che non ha senso.

Per evitare questo, si usa una trasformazione matematica chiamata logit, che ha due proprietà molto utili:

• accetta in ingresso solo valori tra 0 e 1 (cioè le probabilità),
• restituisce un numero reale qualsiasi, da \(-\infty\) a \(+\infty\).

La trasformazione logit è definita così:

\[ \text{logit}(p_i) = \log\left(\frac{p_i}{1 - p_i}\right) . \]

Questa quantità si chiama log-odds e rappresenta il logaritmo del rapporto tra la probabilità di successo e quella di insuccesso.

Esempio:

• se \(p_i = 0.5\), allora \(\text{logit}(p_i) = \log(1) = 0\);
• se \(p_i = 0.8\), allora \(\text{logit}(p_i) = \log(4) \approx 1.39\);
• se \(p_i = 0.2\), allora \(\text{logit}(p_i) = \log(0.25) \approx -1.39\).

Grazie a questa trasformazione, possiamo costruire un modello lineare senza rischiare di ottenere valori fuori dall’intervallo [0,1].

72.4 Il modello di regressione logistica

Mettiamo insieme tutti i pezzi. Il nostro modello diventa:

\[ \begin{aligned} y_i &\sim \text{Bernoulli}(p_i), \\ \text{logit}(p_i) &= \alpha + \gamma D_i. \end{aligned} \]

Vediamo cosa significano i due parametri del modello:

• \(\alpha\) è il log-odds di successo per il gruppo 0 (quando \(D_i = 0\)).

  • La probabilità di successo corrispondente si ottiene con la funzione logistica:

    \[ p_0 = \frac{e^{\alpha}}{1 + e^{\alpha}} = \text{logistic}(\alpha) . \]

• \(\gamma\) rappresenta la differenza nei log-odds tra il gruppo 1 e il gruppo 0.

  • Il log-odds nel gruppo 1 è quindi $+ $, e la probabilità è:

    \[ p_1 = \frac{e^{\alpha + \gamma}}{1 + e^{\alpha + \gamma}} = \text{logistic}(\alpha + \gamma) . \]

72.4.1 Interpretazione pratica

• Se \(\gamma > 0\), il gruppo 1 ha una probabilità di successo più alta rispetto al gruppo 0.
• Se \(\gamma < 0\), il gruppo 1 ha una probabilità più bassa.
• Se \(\gamma = 0\), i due gruppi hanno la stessa probabilità di successo.

72.4.2 Vantaggi dell’approccio bayesiano

In un’analisi bayesiana, non ci limitiamo a stimare un singolo valore per \(\gamma\). Invece, otteniamo una distribuzione completa che rappresenta tutte le ipotesi plausibili sui valori di \(\gamma\), tenendo conto:

• della variabilità nei dati,
• delle nostre ipotesi iniziali (le distribuzioni a priori),
• e delle informazioni che emergono dai dati osservati.

Questo ci permette di rispondere a domande come:

• Quanto è probabile che \(\gamma\) sia maggiore di 0?
• Qual è l’intervallo più credibile in cui può trovarsi \(\gamma\) con il 95% di probabilità?
• Qual è la probabilità che la differenza tra i gruppi sia sostanziale, non solo presente?

72.4.3 Inferenza bayesiana

1. Scelta delle prior

   • Un’opzione comune è usare prior debolmente informative, ad esempio \(\alpha\sim\mathcal N(0,\,2.5)\) e \(\gamma\sim\mathcal N(0,\,2.5)\). Queste varianze larghe lasciano che i dati “parlino”, ma impediscono che le probabilità si avvicinino troppo a 0 o 1 senza evidenza.

2. Calcolo della distribuzione a posteriori

   • Si usa normalmente l’algoritmo MCMC (per es. No‑U‑Turn Sampler di Stan).
   • Otteniamo campioni \(\{\alpha^{(s)},\gamma^{(s)}\}_{s=1}^S\).

3. Quantità derivate di interesse

   • Probabilità nei due gruppi: \(p_0^{(s)}=\operatorname{logistic}(\alpha^{(s)})\), \(p_1^{(s)}=\operatorname{logistic}(\alpha^{(s)}+\gamma^{(s)})\).
   • Differenza di probabilità: \(\Delta^{(s)} = p_1^{(s)}-p_0^{(s)}\). Mostra quanto, in media, il gruppo 1 supera (o non supera) il gruppo 0 in termini di proporzione.
   • Rapporto di odds: \(\text{OR}^{(s)} = e^{\gamma^{(s)}}\). Se \(\text{OR}=2\) significa che gli odds di successo nel gruppo 1 sono il doppio di quelli nel gruppo 0.

4. Sintesi dei risultati

   • Media (o mediana) a posteriori per \(p_0, p_1, \Delta, \text{OR}\).
   • Intervalli di credibilità al 95 %: tagliamo il 2.5 % di campioni in ciascuna coda.
   • Probabilità che \(\Delta>0\) o che \(\text{OR}>1\): basta contare la frazione di campioni corrispondenti.

72.4.4 Perché tutto questo funziona?

• Il logit “apre” l’intervallo (0, 1) rendendo possibile usare un modello lineare.

• La regressione logistica con una dummy è l’esatto equivalente, in termini di parametri, al test bayesiano sulle due proporzioni, ma:

  • consente estensioni (più covariate, effetti casuali, interazioni);
  • permette di riportare risultati direttamente interpretabili (differenza di probabilità, OR, predicted probabilities).

• L’approccio bayesiano produce output che si leggono come “date le nostre ipotesi preliminari e i dati, la plausibilità che la vera differenza di probabilità stia in questo intervallo è 95 %”, concetto spesso più intuitivo del valore‑p.

72.5 Inferenza sulle proporzioni

Il confronto tra le proporzioni di due gruppi indipendenti può essere affrontato sia con un approccio frequentista, basato sulla distribuzione campionaria, sia con un approccio bayesiano. Vediamo i due approcci in dettaglio.

72.6 Approccio Frequentista

Quando vogliamo confrontare le probabilità di successo in due gruppi distinti, possiamo analizzare la differenza tra le proporzioni di successo osservate. L’approccio frequentista affronta il problema studiando la distribuzione campionaria di questa differenza.

72.6.1 Modello di riferimento

Supponiamo di avere due gruppi indipendenti. In ciascun gruppo osserviamo esiti binari, come “successo” (1) o “insuccesso” (0), per ogni partecipante. Formalmente:

\[ Y_1 \sim \text{Bernoulli}(p_1), \quad Y_2 \sim \text{Bernoulli}(p_2), \]

dove \(p_1\) e \(p_2\) sono le probabilità di successo nella popolazione del primo e del secondo gruppo, rispettivamente.

72.6.2 Obiettivo dell’inferenza

Siamo interessati a stimare e fare inferenza sulla differenza tra le due proporzioni:

\[ \Delta = p_1 - p_2. \]

Poiché non conosciamo \(p_1\) e \(p_2\), li stimiamo usando le proporzioni campionarie:

\[ \hat{p}_1 = \frac{X_1}{n_1}, \quad \hat{p}_2 = \frac{X_2}{n_2}, \quad \text{e dunque} \quad \hat{\Delta} = \hat{p}_1 - \hat{p}_2. \]

72.6.3 Proprietà della distribuzione campionaria

Per capire se la differenza osservata tra \(\hat{p}_1\) e \(\hat{p}_2\) è attribuibile al caso oppure riflette una differenza reale tra i gruppi, analizziamo la distribuzione campionaria di \(\hat{\Delta}\).

72.6.3.1 Valore atteso

Il valore atteso della differenza stimata è:

\[ E(\hat{p}_1 - \hat{p}_2) = p_1 - p_2, \]

cioè, in media, la stima è corretta (è uno stimatore non distorto).

Dimostrazione

Siano \(X_1\) il numero di successi in un campione di dimensione \(n_1\) estratto dalla prima popolazione, e \(X_2\) il numero di successi in un campione di dimensione \(n_2\) estratto dalla seconda popolazione. Assumiamo che i campioni siano indipendenti.

Gli stimatori delle proporzioni delle popolazioni \(p_1\) e \(p_2\) sono le proporzioni campionarie, definite come:

\[\hat{p}_1 = \frac{X_1}{n_1}\]

\[\hat{p}_2 = \frac{X_2}{n_2}\]

Sappiamo che \(X_1\) segue una distribuzione binomiale con parametri \(n_1\) e \(p_1\), e \(X_2\) segue una distribuzione binomiale con parametri \(n_2\) e \(p_2\).

Il valore atteso di una variabile casuale binomiale \(X \sim B(n, p)\) è \(E(X) = np\). Quindi:

\[E(X_1) = n_1 p_1\]

\[E(X_2) = n_2 p_2\]

Ora, consideriamo il valore atteso di \(\hat{p}_1\):

\[E(\hat{p}_1) = E\left(\frac{X_1}{n_1}\right)\]

Per la proprietà di linearità del valore atteso (\(E(cX) = cE(X)\) dove \(c\) è una costante):

\[E(\hat{p}_1) = \frac{1}{n_1} E(X_1)\]

Sostituendo il valore di \(E(X_1)\):

\[E(\hat{p}_1) = \frac{1}{n_1} (n_1 p_1) = p_1\]

Questo dimostra che \(\hat{p}_1\) è uno stimatore non distorto per \(p_1\).

In modo analogo, per \(\hat{p}_2\):

\[E(\hat{p}_2) = E\left(\frac{X_2}{n_2}\right)\]

\[E(\hat{p}_2) = \frac{1}{n_2} E(X_2)\]

Sostituendo il valore di \(E(X_2)\):

\[E(\hat{p}_2) = \frac{1}{n_2} (n_2 p_2) = p_2\]

Questo dimostra che \(\hat{p}_2\) è uno stimatore non distorto per \(p_2\).

Ora vogliamo trovare il valore atteso della differenza tra gli stimatori, \(E(\hat{p}_1 - \hat{p}_2)\). Per la proprietà di linearità del valore atteso (\(E(X-Y) = E(X) - E(Y)\)):

\[E(\hat{p}_1 - \hat{p}_2) = E(\hat{p}_1) - E(\hat{p}_2)\]

Sostituendo i valori che abbiamo trovato per \(E(\hat{p}_1)\) e \(E(\hat{p}_2)\):

\[E(\hat{p}_1 - \hat{p}_2) = p_1 - p_2\]

Quindi, abbiamo dimostrato che il valore atteso della differenza tra le proporzioni campionarie è uguale alla differenza tra le proporzioni delle popolazioni. Questo significa che \(\hat{p}_1 - \hat{p}_2\) è uno stimatore non distorto per \(p_1 - p_2\).

72.6.3.2 Varianza

Assumendo che i campioni siano indipendenti, la varianza della differenza stimata è la somma delle varianze delle due proporzioni:

\[ \operatorname{Var}(\hat{p}_1 - \hat{p}_2) = \frac{p_1(1 - p_1)}{n_1} + \frac{p_2(1 - p_2)}{n_2}. \]

Questa formula ci dice quanto può variare la differenza stimata da un campione all’altro.

Poiché \(p_1\) e \(p_2\) sono ignoti, nella pratica li sostituiamo con le proporzioni osservate \(\hat{p}_1\) e \(\hat{p}_2\) per stimare la varianza.

Dimostrazione

Come visto in precedenza, gli stimatori delle proporzioni delle popolazioni sono:

\[ \hat{p}_1 = \frac{X_1}{n_1}, \]

\[ \hat{p}_2 = \frac{X_2}{n_2}. \]

dove \(X_1 \sim B(n_1, p_1)\) e \(X_2 \sim B(n_2, p_2)\) sono il numero di successi nei rispettivi campioni.

Vogliamo calcolare \(\operatorname{Var}(\hat{p}_1 - \hat{p}_2)\).

Per la proprietà della varianza della differenza di due variabili casuali indipendenti, \(\operatorname{Var}(Y - Z) = \operatorname{Var}(Y) + \operatorname{Var}(Z)\), poiché i campioni (e quindi \(X_1\) e \(X_2\), e di conseguenza \(\hat{p}_1\) e \(\hat{p}_2\)) sono assunti indipendenti.

Quindi, possiamo scrivere:

\[ \operatorname{Var}(\hat{p}_1 - \hat{p}_2) = \operatorname{Var}(\hat{p}_1) + \operatorname{Var}(\hat{p}_2).\]

Ora, dobbiamo trovare la varianza di una proporzione campionaria. Consideriamo \(\hat{p} = X/n\), dove \(X \sim B(n, p)\). La varianza di una variabile casuale binomiale \(X \sim B(n, p)\) è \(\operatorname{Var}(X) = np(1-p)\). Usando la proprietà della varianza \(\operatorname{Var}(cX) = c^2\operatorname{Var}(X)\), dove \(c\) è una costante:

\[ \operatorname{Var}(\hat{p}) = \operatorname{Var}\left(\frac{X}{n}\right) = \left(\frac{1}{n}\right)^2 \operatorname{Var}(X). \]

Sostituendo la varianza di \(X\):

\[ \operatorname{Var}(\hat{p}) = \frac{1}{n^2} (np(1-p)) = \frac{p(1-p)}{n}. \]

Applichiamo questa formula per trovare \(\operatorname{Var}(\hat{p}_1)\) e \(\operatorname{Var}(\hat{p}_2)\):

\[ \operatorname{Var}(\hat{p}_1) = \frac{p_1(1 - p_1)}{n_1}, \]

\[ \operatorname{Var}(\hat{p}_2) = \frac{p_2(1 - p_2)}{n_2}. \]

Ora sostituiamo queste varianze nell’espressione per \(\operatorname{Var}(\hat{p}_1 - \hat{p}_2)\):

\[ \operatorname{Var}(\hat{p}_1 - \hat{p}_2) = \frac{p_1(1 - p_1)}{n_1} + \frac{p_2(1 - p_2)}{n_2}. \]

Questo completa la dimostrazione. La formula ottenuta mostra che la varianza della differenza tra le proporzioni campionarie è la somma delle varianze individuali delle proporzioni campionarie, il che ha senso dato l’assunto di indipendenza dei campioni.

Come menzionato nel testo, poiché i veri valori delle proporzioni di popolazione \(p_1\) e \(p_2\) sono generalmente ignoti, per stimare la varianza nella pratica si sostituiscono \(p_1\) e \(p_2\) con le loro stime campionarie \(\hat{p}_1\) e \(\hat{p}_2\). Lo stimatore della varianza diventa quindi:

\[ \widehat{\operatorname{Var}}(\hat{p}_1 - \hat{p}_2) = \frac{\hat{p}_1(1 - \hat{p}_1)}{n_1} + \frac{\hat{p}_2(1 - \hat{p}_2)}{n_2}. \]

72.6.4 Approssimazione normale

Quando i campioni sono sufficientemente grandi, possiamo applicare il Teorema del Limite Centrale, che ci assicura che la distribuzione della differenza \(\hat{p}_1 - \hat{p}_2\) si avvicina a una distribuzione normale:

\[ \hat{p}_1 - \hat{p}_2 \sim \mathcal{N}\left(p_1 - p_2,\ \sqrt{\frac{p_1(1 - p_1)}{n_1} + \frac{p_2(1 - p_2)}{n_2}}\right). \]

Questa approssimazione è alla base della costruzione di:

• intervalli di confidenza per \(p_1 - p_2\),
• test di ipotesi per verificare se la differenza tra le proporzioni è nulla.

72.7 Un esempio illustrativo: studio su bambini lavoratori e scolari

Per mettere in pratica quanto detto, consideriamo un caso reale tratto dallo studio di Banerjee et al. (2025), che ha confrontato le abilità matematiche di:

• bambini lavoratori nei mercati di Kolkata e Delhi;
• bambini scolarizzati che non lavorano.

L’obiettivo era valutare se le competenze sviluppate nel lavoro quotidiano (come dare il resto, sommare prezzi, ecc.) si trasferiscono al contesto scolastico e viceversa.

72.7.1 Risultati principali

• I bambini lavoratori si sono dimostrati molto abili nel risolvere problemi matematici concreti, ma hanno faticato con problemi presentati in forma astratta (come quelli scolastici).
• Al contrario, i bambini scolarizzati se la sono cavata meglio con problemi astratti, ma sono risultati poco efficaci nei problemi concreti del mercato.

Vediamo i dati raccolti per due tipi di problemi:

72.7.2 Problemi astratti

Gruppo	Successi	Prove totali	Proporzione
Bambini lavoratori	670	1488	0.45
Bambini scolarizzati	320	542	0.59

72.7.3 Problemi di mercato

Gruppo	Successi	Prove totali	Proporzione
Bambini lavoratori	134	373	0.36
Bambini scolarizzati	3	271	0.01

72.8 Applicazione dell’approccio frequentista

72.8.1 Confronto per i problemi astratti

Vogliamo verificare se la differenza tra 0.45 (lavoratori) e 0.59 (scolarizzati) è spiegabile dal caso.

Ipotesi.

• \(H_0\): \(p_1 = p_2\) (nessuna differenza tra i gruppi);
• \(H_1\): \(p_1 \ne p_2\) (esiste una differenza).

Calcolo.

1. Proporzione combinata (pooled):

   \[ \hat{p} = \frac{670 + 320}{1488 + 542} = \frac{990}{2030} \approx 0.487. \]

2. Varianza stimata della differenza:

   \[ \text{Var} = \hat{p}(1 - \hat{p}) \left( \frac{1}{1488} + \frac{1}{542} \right) \approx 0.000629. \]

3. Deviazione standard: $ $.

4. Statistica z:

   \[ z = \frac{0.45 - 0.59}{0.0251} \approx -5.58. \]

Il valore z è molto lontano da 0: la probabilità di osservare una tale differenza per caso (p-value) è inferiore a 0.0001. Possiamo quindi rifiutare l’ipotesi nulla.

72.8.2 Confronto per i problemi di mercato

Dati:

• Lavoratori: 134 su 373 (\(\hat{p}_1 = 0.36\))
• Scolarizzati: 3 su 271 (\(\hat{p}_2 = 0.01\))

Ripetiamo i passaggi:

1. \(\hat{p} = \frac{137}{644} \approx 0.213\)

2. \(\text{Var} \approx 0.001067\)

3. Deviazione standard: \(\sqrt{0.001067} \approx 0.0327\)

4. Statistica z:

   \[ z = \frac{0.36 - 0.01}{0.0327} \approx 10.70 \]

Anche qui, il p-value è praticamente zero: le differenze sono molto più grandi di quanto ci si aspetti per puro caso.

In sintesi, l’approccio frequentista ci consente di:

• stimare la differenza tra le proporzioni,
• quantificare l’incertezza (varianza e intervallo di confidenza),
• testare ipotesi sul fatto che la differenza sia zero o meno.

In entrambi i confronti (problemi astratti e problemi concreti), abbiamo trovato evidenze chiare di una differenza tra i due gruppi.

72.8.3 Svolgimento con R

Di seguito mostriamo due modalità per replicare i calcoli in R: (1) passo passo usando le formule manuali, (2) usando la funzione prop.test() di R. Verranno illustrate entrambe le analisi: quella per i problemi astratti e quella per i problemi matematici di mercato.

72.8.3.1 Problemi astratti

Dati

• Bambini lavoratori (gruppo 1): 670 successi su 1488 prove.
• Bambini scolarizzati (gruppo 2): 320 successi su 542 prove.

# Dati
X1 <- 670
n1 <- 1488
X2 <- 320
n2 <- 542

# Proporzioni campionarie
p1 <- X1 / n1
p2 <- X2 / n2

# Proporzione pooled (combinata)
p_pool <- (X1 + X2) / (n1 + n2)

# Differenza fra le proporzioni
diff_p <- p1 - p2

# Calcolo della varianza della differenza (usando la formula per due proporzioni)
var_diff <- p_pool * (1 - p_pool) * (1/n1 + 1/n2)

# Deviazione standard della differenza
sd_diff <- sqrt(var_diff)

# Statistica z
z_value <- diff_p / sd_diff

# p-value (test a due code)
p_value <- 2 * pnorm(abs(z_value), lower.tail = FALSE)

# Visualizzazione risultati
z_value
#> [1] -5.588
p_value
#> [1] 2.295e-08

72.8.3.2 Analisi con funzione prop.test()

Per ottenere direttamente il test di confronto di due proporzioni, possiamo usare la funzione prop.test di R. Attenzione che, di default, prop.test effettua una correzione per la continuità (Yates), che in questo contesto disattiviamo per confrontare i risultati con i calcoli manuali (impostando correct = FALSE).

# Confronto due proporzioni senza correzione per la continuità
test_astratti <- prop.test(
  x = c(X1, X2),
  n = c(n1, n2),
  alternative = "two.sided",
  correct = FALSE
)

test_astratti
#> 
#>  2-sample test for equality of proportions without continuity
#>  correction
#> 
#> data:  c(X1, X2) out of c(n1, n2)
#> X-squared = 31, df = 1, p-value = 2e-08
#> alternative hypothesis: two.sided
#> 95 percent confidence interval:
#>  -0.18864 -0.09163
#> sample estimates:
#> prop 1 prop 2 
#> 0.4503 0.5904

72.8.3.3 Problemi matematici di mercato

# Dati
X1_m <- 134
n1_m <- 373
X2_m <- 3
n2_m <- 271

# Proporzioni campionarie
p1_m <- X1_m / n1_m
p2_m <- X2_m / n2_m

# Proporzione pooled
p_pool_m <- (X1_m + X2_m) / (n1_m + n2_m)

# Differenza di proporzioni
diff_p_m <- p1_m - p2_m

# Varianza della differenza (usando la proporzione pooled)
var_diff_m <- p_pool_m * (1 - p_pool_m) * (1/n1_m + 1/n2_m)

# Deviazione standard
sd_diff_m <- sqrt(var_diff_m)

# Statistica z
z_value_m <- diff_p_m / sd_diff_m

# p-value (test a due code)
p_value_m <- 2 * pnorm(abs(z_value_m), lower.tail = FALSE)

# Visualizzazione risultati
z_value_m
#> [1] 10.66
p_value_m
#> [1] 1.581e-26

Analogamente a quanto fatto per i problemi astratti:

# Confronto due proporzioni
test_mercato <- prop.test(
  x = c(X1_m, X2_m),
  n = c(n1_m, n2_m),
  alternative = "two.sided",
  correct = FALSE
)

test_mercato
#> 
#>  2-sample test for equality of proportions without continuity
#>  correction
#> 
#> data:  c(X1_m, X2_m) out of c(n1_m, n2_m)
#> X-squared = 114, df = 1, p-value <2e-16
#> alternative hypothesis: two.sided
#> 95 percent confidence interval:
#>  0.2979 0.3984
#> sample estimates:
#>  prop 1  prop 2 
#> 0.35925 0.01107

72.8.3.3.1 Interpretazione dei risultati

1. Problemi astratti
   • Statistica z (calcolo manuale): circa -5.58
     
   • p-value: molto piccolo (<< 0.001)
     
   • Conclusione: rifiutiamo \(H_0\) e concludiamo che la differenza tra le proporzioni di successo dei bambini lavoratori e scolarizzati sono degne di nota.
2. Problemi matematici di mercato
   • Statistica z (calcolo manuale): circa 10.70
     
   • p-value: estremamente piccolo (<< 0.001)
     
   • Conclusione: rifiutiamo \(H_0\) e concludiamo che anche in questo caso le differenze tra le proporzioni di successo dei due gruppi sono degne di nota.

72.9 Approccio Bayesiano

Come indicato nell’introduzione di questo capitolo, è possibile applicare l’approccio bayesiano al problema dell’inferenza sulla differenza tra due proporzioni indipendenti utilizzando un modello di regressione con una variabile indicatrice (dummy) per distinguere i due gruppi. Per fare un esempio, consideriamo qui i problemi astratti. Utilizziamo il pacchetto brms per stimare il modello, con distribuzioni a priori debolmente informative specificate di default:

X1 <- 670
n1 <- 1488
X2 <- 320
n2 <- 542

dat_a <- data.frame(
  count = c(X1, X2),
  tot = c(n1, n2),
  group = c("working", "non-working")
)
dat_a

fit_a <- brm(count | trials(tot) ~ group, data = dat_a, family = binomial())

summary(fit_a)
#>  Family: binomial 
#>   Links: mu = logit 
#> Formula: count | trials(tot) ~ group 
#>    Data: dat_a (Number of observations: 2) 
#>   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
#>          total post-warmup draws = 4000
#> 
#> Regression Coefficients:
#>              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept        0.37      0.09     0.19     0.54 1.00     1561     2015
#> groupworking    -0.57      0.10    -0.77    -0.36 1.00     1873     1865
#> 
#> 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).

# Extract posterior draws of the Intercept and groupworking coefficients
post <- posterior_samples(fit_a) 
# or posterior_draws(fit_a) in newer versions

# Probability for reference group = logistic(Intercept)
post$p_ref <- plogis(post$b_Intercept)

# Probability for working group = logistic(Intercept + groupworking)
post$p_work <- plogis(post$b_Intercept + post$b_groupworking)

# Now summarize
summary_ref  <- quantile(post$p_ref,  probs = c(0.025, 0.5, 0.975))
summary_work <- quantile(post$p_work, probs = c(0.025, 0.5, 0.975))

summary_ref
#>   2.5%    50%  97.5% 
#> 0.5471 0.5907 0.6321
summary_work
#>   2.5%    50%  97.5% 
#> 0.4243 0.4504 0.4762

# 1. Extract samples
post <- as_draws_df(fit_a) 

# 2. Summaries on logit 
# b_groupworking
b_groupworking_CI <- quantile(post$b_groupworking, probs = c(0.025, 0.5, 0.975))
print("b_groupworking (log-odds) [2.5%, 50%, 97.5%]:")
#> [1] "b_groupworking (log-odds) [2.5%, 50%, 97.5%]:"
b_groupworking_CI
#>    2.5%     50%   97.5% 
#> -0.7707 -0.5670 -0.3612

# 3. Summaries on OR
post$OR_groupworking <- exp(post$b_groupworking)
OR_groupworking_CI <- quantile(post$OR_groupworking, probs = c(0.025, 0.5, 0.975))
print("Odds Ratio for 'working' vs. reference [2.5%, 50%, 97.5%]:")
#> [1] "Odds Ratio for 'working' vs. reference [2.5%, 50%, 97.5%]:"
OR_groupworking_CI
#>   2.5%    50%  97.5% 
#> 0.4627 0.5672 0.6968

# 4. Summaries on prob. scale

# Probability reference group
post$p_ref  <- plogis(post$b_Intercept)

# Probability working group
post$p_work <- plogis(post$b_Intercept + post$b_groupworking)

# Difference on probability scale
post$diff_p_work_ref <- post$p_work - post$p_ref

summary_ref  <- quantile(post$p_ref,  probs = c(0.025, 0.5, 0.975))
summary_work <- quantile(post$p_work, probs = c(0.025, 0.5, 0.975))
diff_p_CI    <- quantile(post$diff_p_work_ref, probs = c(0.025, 0.5, 0.975))

cat("\nReference group probability [2.5%, 50%, 97.5%]:\n")
#> 
#> Reference group probability [2.5%, 50%, 97.5%]:
print(summary_ref)
#>   2.5%    50%  97.5% 
#> 0.5471 0.5907 0.6321

cat("\nWorking group probability [2.5%, 50%, 97.5%]:\n")
#> 
#> Working group probability [2.5%, 50%, 97.5%]:
print(summary_work)
#>   2.5%    50%  97.5% 
#> 0.4243 0.4504 0.4762

cat("\nDifference in probability (Working - Reference) [2.5%, 50%, 97.5%]:\n")
#> 
#> Difference in probability (Working - Reference) [2.5%, 50%, 97.5%]:
print(diff_p_CI)
#>     2.5%      50%    97.5% 
#> -0.18948 -0.14056 -0.08989

Si noti come, usando prior debolmente informativi, i risultati ottenuti con i due approcci (frequentista e bayesiano) siano praticamente equivalenti.

Come abbiamo osservato nel caso del confronto tra medie, l’approccio bayesiano è più utile perché:

1. non si basa su un’ipotesi nulla idealizzata che sappiamo non essere vera;
2. stima direttamente la nostra incertezza rispetto alla differenza tra le proporzioni, fornendo una risposta concreta e interpretabile in termini di probabilità;
3. permette di integrare le evidenze dei dati con la conoscenza preesistente, per ottenere inferenze più realistiche e informate.

72.9.1 Intervallo di credibilità

Per ottenere l’intervallo di credibilità (Highest Density Interval, HDI) sulla scala delle probabilità (e non su quella logit), è necessario trasformare manualmente i draw a livello di probabilità e poi calcolare l’HDI su quei valori trasformati. In altre parole:

1. estraiamo i draw posteriori di b_Intercept e b_groupworking;
   
2. trasformiamo i valori con la funzione logit-inversa (\(\operatorname{logistic}(x) = 1/(1+e^{-x})\)) per ottenere le probabilità;
   
3. se ci concentriamo sul solo effetto sulla scala della probabilità (e.g. differenza fra i due gruppi), calcoliamo la differenza tra la probabilità del gruppo “working” e quella del gruppo “reference” per ciascun draw;
   
4. applichiamo hdi() su queste grandezze trasformate.

Di seguito è fornito il codice R con il workflow completo.

1. Estrazione draw posteriori.

post <- as_draws_df(fit_a)

2. Calcolo delle probabilità per ciascun draw. La variabile ‘group’ abbia due livelli:

   • “working” (effetto => b_groupworking)
   • “non-working” (riferimento => b_Intercept)

post <- post %>%
  dplyr::mutate(
    p_ref      = plogis(b_Intercept),                       # probabilità (referenza)
    p_working  = plogis(b_Intercept + b_groupworking),      # prob. gruppo working
    diff_working_ref = p_working - p_ref                    # differenza
  )

3. Calcolo dell’HDI sull’effetto (o sulle probabilità).

1. HDI per la probabilità del gruppo “working”.

hdi_working_prob <- hdi(post$p_working, ci = 0.95)
hdi_working_prob
#> 95% HDI: [0.42, 0.48]

2. HDI per la probabilità del gruppo “non-working” (riferimento).

hdi_ref_prob <- hdi(post$p_ref, ci = 0.95)
hdi_ref_prob
#> 95% HDI: [0.55, 0.64]

3. HDI della differenza fra le due probabilità.

hdi_diff <- hdi(post$diff_working_ref, ci = 0.89)
hdi_diff
#> 89% HDI: [-0.18, -0.10]

Interpretazione

• hdi_working_prob fornisce l’HDI al 95% (o al livello che specifichi) della probabilità di “successo” del gruppo “working”.
  
• hdi_ref_prob fa lo stesso per il gruppo di riferimento.
  
• hdi_diff restituisce l’HDI della differenza in probabilità tra “working” e “reference” (\(p_{\text{working}} - p_{\text{reference}}\)).

In questo modo ottieniamo l’intervallo di credibilità (HDI) sulla scala delle probabilità.

72.9.2 Distribuzione Predittiva a Posteriori

Effettuiamo un controllo grafico per confrontare i dati osservati con quelli predetti dal modello:

pp_check(fit_a)

72.10 Riflessioni Conclusive

In questo capitolo abbiamo esplorato il confronto tra due proporzioni adottando sia l’approccio frequentista sia quello bayesiano. L’obiettivo principale era valutare se la proporzione di successi in un gruppo differisse da quella osservata nell’altro gruppo, e con quale grado di incertezza.

Nella procedura classica frequentista (test per due proporzioni), si calcola una statistica di test (ad esempio, un test z) e il relativo p-value, ossia la probabilità di osservare un risultato così estremo (o più) assumendo che le due proporzioni reali siano uguali (ipotesi nulla).

• Intervallo di confidenza: È possibile costruire un intervallo di confidenza (IC) per la differenza tra le due proporzioni. L’interpretazione frequentista di tale IC, però, si basa su un’ipotetica ripetizione di campionamenti ed è focalizzata sull’eventuale rifiuto o meno dell’ipotesi che la differenza sia zero.
  
• Limiti interpretativi: L’approccio frequentista si fonda sul concetto di ipotesi nulla “nessuna differenza” e non fornisce una probabilità diretta di quanto la differenza vera sia maggiore o minore di un certo valore, limitandosi a indicare se i dati sono inusuali qualora la differenza fosse zero.

Grazie alla regola di Bayes, combiniamo informazioni a priori (sul probabile valore delle proporzioni) con i dati osservati, per ottenere una distribuzione a posteriori della differenza tra le due proporzioni. Questa distribuzione descrive i valori plausibili della differenza, insieme alle relative credibilità (probabilità).

• Credible Interval o Highest Density Interval (HDI): Al posto di un intervallo di confidenza, l’approccio bayesiano fornisce un intervallo di credibilità. Ad esempio, un 95% HDI indica i valori della differenza tra le proporzioni che cumulativamente contengono il 95% della probabilità a posteriori. È un costrutto immediatamente interpretabile: “Abbiamo una probabilità del 95% che la differenza vera cada all’interno di questo intervallo”.
  
• Flessibilità e interpretazione diretta: L’approccio bayesiano permette di rispondere in modo più naturale a domande come: “Qual è la probabilità che la differenza fra le due proporzioni sia maggiore di 0?” oppure “Qual è la probabilità che la proporzione di un gruppo superi quella dell’altro di almeno una certa soglia rilevante?”.

Confronto tra i due approcci:

1. Interpretazione dei risultati: Il p-value frequentista ci dice quanto il dato sia “improbabile” sotto l’ipotesi di uguaglianza delle proporzioni; il Bayesianesimo risponde direttamente a quanto è plausibile ogni possibile valore di differenza.
   
2. Centralità dell’ipotesi nulla: Nel frequentismo, l’ipotesi nulla (differenza = 0) è centrale. Nel modello bayesiano, è invece possibile assegnare direttamente probabilità alla differenza e alla sua distanza da zero, evitando un focus eccessivo sull’uguaglianza perfetta delle due proporzioni.
   
3. Ruolo dei priors: L’uso di priors (non informativi o informativi) può influire sulle stime bayesiane quando i dati sono scarsi, rendendo evidente la necessità di scelte trasparenti e ben motivate. Tuttavia, con campioni ampi, l’influenza dei priors tende a ridursi e la stima a posteriori è dominata dai dati.
   
4. Completezza dell’inferenza: L’approccio bayesiano consente di integrare nuove informazioni e di aggiornare la distribuzione a posteriori man mano che arrivano dati aggiuntivi. Al contrario, l’approccio frequentista non fornisce un meccanismo diretto di “aggiornamento” delle stime alla luce di nuovi dati.

In sintesi,

• approccio frequentista: Si tratta di una metodologia consolidata e standard nella ricerca; fornisce risultati in termini di p-value e IC, ma l’interpretazione del p-value e dell’IC resta legata a procedure di campionamento ipotetico.
  
• approccio bayesiano: Offre una maniera più intuitiva di quantificare l’incertezza, assegnando probabilità dirette ai possibili valori di differenza fra le due proporzioni. Consente di formulare domande più specifiche (es. la probabilità che la differenza superi un valore definito) e di integrare in modo naturale informazioni a priori.

Nella pratica della ricerca, l’approccio frequentista rimane diffuso. Tuttavia, l’inferenza bayesiana fornisce un quadro interpretativo più ricco e flessibile. Utilizzare entrambi i metodi, quando appropriato, può potenziare l’analisi e la comprensione dei dati, permettendo di trarre conclusioni più robuste e trasparenti.

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] insight_1.2.0     bayestestR_0.15.3 posterior_1.6.1   cmdstanr_0.9.0   
#>  [5] brms_2.22.0       Rcpp_1.0.14       thematic_0.1.6    MetBrewer_0.2.0  
#>  [9] ggokabeito_0.1.0  see_0.11.0        gridExtra_2.3     patchwork_1.3.0  
#> [13] bayesplot_1.12.0  psych_2.5.3       scales_1.4.0      markdown_2.0     
#> [17] knitr_1.50        lubridate_1.9.4   forcats_1.0.0     stringr_1.5.1    
#> [21] dplyr_1.1.4       purrr_1.0.4       readr_2.1.5       tidyr_1.3.1      
#> [25] tibble_3.2.1      ggplot2_3.5.2     tidyverse_2.0.0   rio_1.2.3        
#> [29] here_1.0.1       
#> 
#> loaded via a namespace (and not attached):
#>  [1] tidyselect_1.2.1     farver_2.1.2         loo_2.8.0           
#>  [4] fastmap_1.2.0        tensorA_0.36.2.1     pacman_0.5.1        
#>  [7] digest_0.6.37        timechange_0.3.0     lifecycle_1.0.4     
#> [10] StanHeaders_2.32.10  processx_3.8.6       magrittr_2.0.3      
#> [13] compiler_4.5.0       rlang_1.1.6          tools_4.5.0         
#> [16] yaml_2.3.10          labeling_0.4.3       bridgesampling_1.1-2
#> [19] htmlwidgets_1.6.4    curl_6.2.2           pkgbuild_1.4.7      
#> [22] mnormt_2.1.1         plyr_1.8.9           RColorBrewer_1.1-3  
#> [25] abind_1.4-8          withr_3.0.2          grid_4.5.0          
#> [28] stats4_4.5.0         colorspace_2.1-1     inline_0.3.21       
#> [31] cli_3.6.5            mvtnorm_1.3-3        rmarkdown_2.29      
#> [34] generics_0.1.3       RcppParallel_5.1.10  rstudioapi_0.17.1   
#> [37] reshape2_1.4.4       tzdb_0.5.0           rstan_2.32.7        
#> [40] parallel_4.5.0       matrixStats_1.5.0    vctrs_0.6.5         
#> [43] V8_6.0.3             Matrix_1.7-3         jsonlite_2.0.0      
#> [46] callr_3.7.6          hms_1.1.3            glue_1.8.0          
#> [49] codetools_0.2-20     ps_1.9.1             distributional_0.5.0
#> [52] stringi_1.8.7        gtable_0.3.6         QuickJSR_1.7.0      
#> [55] pillar_1.10.2        htmltools_0.5.8.1    Brobdingnag_1.2-9   
#> [58] R6_2.6.1             rprojroot_2.0.4      evaluate_1.0.3      
#> [61] lattice_0.22-7       backports_1.5.0      rstantools_2.4.0    
#> [64] coda_0.19-4.1        nlme_3.1-168         checkmate_2.3.2     
#> [67] xfun_0.52            pkgconfig_2.0.3

Bibliografia

Banerjee, A. V., Bhattacharjee, S., Chattopadhyay, R., Duflo, E., Ganimian, A. J., Rajah, K., & Spelke, E. S. (2025). Children’s arithmetic skills do not transfer between applied and academic mathematics. Nature, 1–9.