Confronto tra due gruppi

L’obiettivo di questo capitolo è di ampliare la discussione del capitolo bayes_one_mean_notebook, affrontando il confronto tra le medie di due gruppi indipendenti. Per cominciare, carichiamo le librerie necessarie.

Preparazione del Notebook

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import statistics as st
from scipy.stats import gaussian_kde
import pymc as pm
import pymc.sampling_jax
import arviz as az
import warnings

warnings.filterwarnings("ignore", category=FutureWarning)

/Users/corrado/opt/anaconda3/envs/pymc_env/lib/python3.11/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

%config InlineBackend.figure_format = 'retina'
RANDOM_SEED = 8927
rng = np.random.default_rng(RANDOM_SEED)
az.style.use("arviz-darkgrid")

Stima bayesiana e test dell’ipotesi nulla

Spesso, ci troviamo ad affrontare la necessità di confrontare due gruppi di dati. Potrebbe interessarci sapere se la media di un gruppo è maggiore o diversa rispetto a quella di un altro gruppo. Per effettuare tale confronto, è fondamentale utilizzare un modello statistico, poiché le vere differenze tra i gruppi sono spesso accompagnate da rumore di misurazione o fluttuazioni casuali del fenomeno in esame. Questo rende difficile trarre conclusioni basandosi unicamente sulle differenze calcolate dai dati osservati.

Il metodo tradizionale per confrontare statisticamente due o più gruppi è quello di utilizzare un test statistico. Questo approccio prevede l’individuazione di un’ipotesi nulla, che solitamente afferma che non ci sono differenze tra i gruppi, e l’utilizzo di una statistica test per determinare se i dati osservati sono plausibili sotto questa ipotesi. L’ipotesi nulla viene rifiutata quando la statistica test calcolata supera una soglia predefinita.

Tuttavia, i test di ipotesi possono essere complessi e i risultati spesso soggetti a interpretazioni errate. La scelta delle specifiche del test statistico (ad esempio, quale test utilizzare, quale ipotesi nulla testare, quale livello di significatività adottare) è spesso arbitraria e basata su convenzioni piuttosto che sulla specificità del problema o delle decisioni da prendere (Johnson, 1999). Inoltre, i risultati forniti dai test sono spesso indiretti, incompleti e tendono a sovrastimare le evidenze contro l’ipotesi nulla (Goodman, 1999).

Un approccio più informativo ed efficace per il confronto tra gruppi è quello basato sulla stima invece che sul test dell’ipotesi nulla, ed è guidato dalla probabilità bayesiana anziché dalla frequentista. In pratica, invece di testare se ci sono differenze tra i gruppi, si cerca di ottenere una stima di quanto siano effettivamente diversi. Questo approccio è intrinsecamente più informativo. Inoltre, viene inclusa una stima dell’incertezza associata a tale differenza, che tiene conto sia dell’incertezza dovuta alla nostra mancanza di conoscenza dei parametri del modello (incertezza epistemica) sia dell’incertezza causata dalla variabilità intrinseca del sistema (incertezza aleatoria).

Un esempio illustrativo

In questo esempio, l’obiettivo è stimare la differenza tra le medie del quoziente di intelligenza dei bambini di due gruppi distinti in base al livello di scolarità della madre. Il primo gruppo include i bambini la cui madre non ha completato le scuole superiori, mentre il secondo gruppo comprende quelli la cui madre ha ottenuto il diploma superiore. Per questo, useremo i dati kidiq e un modello bayesiano al fine di ottenere una stima affidabile della differenza tra le medie dei due gruppi nella popolazione.

I dati utilizzati sono forniti da Gelman e Hill (2007) e costituiscono un sottocampione estratto dal National Longitudinal Survey of Youth.

Leggiamo i dati.

url = "http://www.stat.columbia.edu/~gelman/arm/examples/child.iq/kidiq.dta"
df = pd.read_stata(url)
df.head()

	kid_score	mom_hs	mom_iq	mom_work	mom_age
0	65	1.0	121.117529	4	27
1	98	1.0	89.361882	4	25
2	85	1.0	115.443165	4	27
3	83	1.0	99.449639	3	25
4	115	1.0	92.745710	4	27

Il dataset contiene le seguenti colonne:

• “kid_score”: il quoziente intellettivo (QI) dei bambini. È una misura dell’intelligenza del bambino.

• “mom_hs”: una variabile binaria che indica se la madre del bambino ha completato o meno la scuola superiore. Può assumere i valori 0 o 1, dove 0 rappresenta “no” (la madre non ha completato la scuola superiore) e 1 rappresenta “sì” (la madre ha completato la scuola superiore).

Ci sono 93 bambini la cui madre non ha completato la scuola superiore e 341 bambini la cui madre ha ottenuto il diploma di scuola superiore.

df.groupby(["mom_hs"]).size()

mom_hs
0.0     93
1.0    341
dtype: int64

Le statistiche descrittive si ottengono nel modo seguente.

summary_stats = [np.mean, st.stdev]
df.groupby(["mom_hs"]).aggregate(summary_stats)

	kid_score		mom_iq		mom_work		mom_age
	mean	stdev	mean	stdev	mean	stdev	mean	stdev
mom_hs
0.0	77.548387	22.573800	91.889152	12.630498	2.322581	1.226175	21.677419	2.727323
1.0	89.319648	19.049483	102.212049	14.848414	3.052786	1.120727	23.087977	2.617453

I bambini la cui madre ha completato le superiori tendono ad avere un QI maggiore di 11.8 punti rispetto ai bambini la cui madre non ha concluso le superiori.

89.319648 - 77.548387

11.771260999999996

Creiamo due vettori che contengono il QI dei bambini dei due gruppi.

# Vector of kid_score when mom_hs is 1
kid_score_mom_hs_1 = df[df["mom_hs"] == 1]["kid_score"]

# Vector of kid_score when mom_hs is 0
kid_score_mom_hs_0 = df[df["mom_hs"] == 0]["kid_score"]

Esaminiamo la distribuzione del QI nei due gruppi.

# Create a KDE for each group
kde_hs_1 = gaussian_kde(kid_score_mom_hs_1)
kde_hs_0 = gaussian_kde(kid_score_mom_hs_0)

# Generate x values for the plot
x_vals = np.linspace(
    min(np.min(kid_score_mom_hs_1), np.min(kid_score_mom_hs_0)),
    max(np.max(kid_score_mom_hs_1), np.max(kid_score_mom_hs_0)),
    100,
)

# Calculate KDE densities for the x values
density_hs_1 = kde_hs_1(x_vals)
density_hs_0 = kde_hs_0(x_vals)

# Create the plot
plt.title("KDE Plot of Kid IQ for Mothers \nwith and without High School Completion")
plt.xlabel("Kid IQ")
plt.ylabel("Density")

plt.plot(x_vals, density_hs_1, label="HS=1")
plt.plot(x_vals, density_hs_0, color="C3", label="HS=0")

# Show the legend
plt.legend();

Dimensione dell’effetto

Nel caso presente, la differenza tra le medie dei due gruppi è di 11.8 punti sulla scala del QI, e potrebbe sembrare un risultato rilevante, considerando che la metrica del QI è facilmente interpretabile. Tuttavia, è importante notare che il test utilizzato in questo studio non è il WISC, che ha una distribuzione normale con media 100 e deviazione standard 15, ma il test PIAT.

In generale, è difficile comprendere il significato di una differenza tra le medie di due gruppi quando viene presentata solo come valore assoluto, soprattutto quando le varianze dei gruppi sono diverse. Per ottenere una misura più informativa, è necessario considerare sia la differenza tra le medie dei gruppi che l’incertezza associata a queste stime delle medie della popolazione. L’indice statistico che soddisfa questo scopo è noto come “dimensione dell’effetto” (effect size).

La dimensione dell’effetto è una misura della forza dell’associazione osservata, che tiene conto sia della grandezza della differenza tra i gruppi attesi che dell’incertezza sui dati. Tra gli indici più comunemente utilizzati per quantificare la dimensione dell’effetto, vi è l’indice \(d\) di Cohen.

Nel caso di due medie, questo indice è dato da:

\[ d={\frac {{\bar {x}}_{1}-{\bar {x}}_{2}}{s}}, \]

laddove

\[ s={\sqrt {\frac {(n_{1}-1)s_{1}^{2}+(n_{2}-1)s_{2}^{2}}{n_{1}+n_{2}-2}}} \]

e la varianza di ciascun gruppo è calcolata come

\[ s_{1}^{2}={\frac {1}{n_{1}-1}}\sum _{i=1}^{n_{1}}(x_{1,i}-{\bar {x}}_{1})^{2}. \]

Solitamente, l’indice \(d\) di Cohen si interpreta usando la metrica seguente:

Dimensione dell’effetto	\(d\)
Very small	0.01
Small	0.20
Medim	0.50
Large	0.80
Very large	1.20
Huge	2.0

Per una trattazione bayesiana della stima della dimensione dell’effetto, si veda [Kru14].

Modello bayesiano

Il modello bayesiano per il confronto tra le medie di due gruppi indipendenti comprende la definizione della verosimiglianza per i dati di ciascun gruppo e la descrizione delle distribuzioni a priori dei parametri rilevanti. Inoltre, in questo caso, abbiamo incluso anche la stima della dimensione dell’effetto, che ci permette di valutare la forza dell’associazione osservata tra i gruppi, tenendo conto dell’incertezza sui dati.

with pm.Model() as model:
    # Prior per la deviazione standard dei gruppi
    group1_std = pm.HalfCauchy("group1_std", beta=30)
    group0_std = pm.HalfCauchy("group0_std", beta=30)

    # Prior per la media dei gruppi
    group1_mean = pm.Normal("group1_mean", mu=df.kid_score.mean(), sigma=20)
    group0_mean = pm.Normal("group0_mean", mu=df.kid_score.mean(), sigma=20)

    # Likelihood per i dati osservati nei gruppi
    mom_hs_1 = pm.Normal(
        "mom_hs_1", mu=group1_mean, sigma=group1_std, observed=kid_score_mom_hs_1
    )
    mom_hs_0 = pm.Normal(
        "mom_hs_0", mu=group0_mean, sigma=group0_std, observed=kid_score_mom_hs_0
    )

    # Calcolare la grandezza dell'effetto
    diff_of_means = pm.Deterministic("difference of means", group1_mean - group0_mean)
    diff_of_stds = pm.Deterministic("difference of stds", group1_std - group0_std)
    effect_size = pm.Deterministic(
        "effect size", diff_of_means / np.sqrt((group1_std**2 + group0_std**2) / 2)
    )

Vediamo passo dopo passo cosa rappresenta ogni parte del modello:

1. Prior per la deviazione standard dei gruppi (group1_std e group0_std): Queste sono distribuzioni Half-Cauchy che rappresentano le deviazioni standard dei punteggi dei bambini nei due gruppi. La scelta della distribuzione Half-Cauchy viene guidata dalla considerazione che le deviazioni standard non possono essere negative. Il parametro “beta” controlla la scala della distribuzione e influisce sull’incertezza associata alla deviazione standard.

2. Prior per la media dei gruppi (group1_mean e group0_mean): Queste sono distribuzioni normali che rappresentano le medie dei punteggi dei bambini nei due gruppi. La media delle distribuzioni normali è inizializzata con il valore medio dei punteggi dei bambini nell’intero dataset (“df.kid_score.mean()”), mentre il parametro “sigma” rappresenta la deviazione standard della distribuzione normale. Questo implica che si assume una distribuzione delle medie attorno alla media dei punteggi di tutti i bambini nel dataset, con una certa incertezza data dalla deviazione standard “sigma”.

3. Likelihood per i dati osservati nei gruppi (mom_hs_1 e mom_hs_0): Queste sono distribuzioni normali osservate che modellano i dati reali dei punteggi dei bambini nei due gruppi, a seconda che le madri abbiano completato o meno la scuola superiore. Queste distribuzioni rappresentano le stime dei punteggi dei bambini nei rispettivi gruppi, basate sulle medie e le deviazioni standard campionate dai dati. L’argomento “observed” indica che stiamo fornendo i dati osservati per le due distribuzioni normali, che vengono confrontati con i valori campionati dai modelli per calcolare la verosimiglianza dei dati rispetto alle stime dei punteggi dei bambini nei gruppi.

Nel modello, dopo aver definito le distribuzioni a priori e la verosimiglianza, viene calcolata la dimensione dell’effetto (effect size). Per calcolare la dimensione dell’effetto, vengono eseguite tre operazioni:

1. Differenza delle medie (diff_of_means): Viene calcolata la differenza tra le medie dei due gruppi (group1_mean - group0_mean). Questa grandezza indica quanto i punteggi del QI dei bambini nel gruppo delle madri con istruzione superiore differiscono dai punteggi dei bambini nel gruppo delle madri senza istruzione superiore.

2. Differenza delle deviazioni standard (diff_of_stds): Viene calcolata la differenza tra le deviazioni standard dei due gruppi (group1_std - group0_std).

3. La dimensione dell’effetto è calcolata dividendo la differenza delle medie (diff_of_means) per la radice quadrata della media delle deviazioni standard al quadrato dei due gruppi (np.sqrt((group1_std**2 + group0_std**2) / 2)). Questo calcolo tiene conto sia della grandezza della differenza tra le medie che della variabilità dei dati, fornendo così un indice comprensibile della dimensione dell’effetto.

pm.model_to_graphviz(model)

Distribuzioni predittive a priori

Prima di eseguire il campionamento, esaminiamo la distribuzione a priori predittiva per valutare se le scelte delle distribuzioni a priori per i parametri del modello siano adeguate.

with model:
    prior_samples = pm.sample_prior_predictive(30)

Sampling: [group0_mean, group0_std, group1_mean, group1_std, mom_hs_0, mom_hs_1]

Osserviamo che le distribuzioni a priori predittive per le medie dei due gruppi appaiono ragionevoli, in quanto contengono i valori delle medie campionate e includono una gamma di valori plausibili, senza essere eccessivamente estese né troppo limitate.

ax = az.plot_dist(prior_samples.prior_predictive["mom_hs_0"])
ax.set_xlim(-100, 300);

ax = az.plot_dist(prior_samples.prior_predictive["mom_hs_1"])
ax.set_xlim(-300, 500);

Campionamento

Eseguiamo ora il campionamento.

with model:
    idata = pm.sampling_jax.sample_numpyro_nuts()

Show code cell output Hide code cell output

Compiling...

Compilation time = 0:00:01.135585

Sampling...

  0%|                                                                                  | 0/2000 [00:00<?, ?it/s]

Compiling.. :   0%|                                                                    | 0/2000 [00:00<?, ?it/s]

  0%|                                                                                  | 0/2000 [00:00<?, ?it/s]

Compiling.. :   0%|                                                                    | 0/2000 [00:00<?, ?it/s]

  0%|                                                                                  | 0/2000 [00:00<?, ?it/s]

Compiling.. :   0%|                                                                    | 0/2000 [00:00<?, ?it/s]

  0%|                                                                                  | 0/2000 [00:00<?, ?it/s]

Compiling.. :   0%|                                                                    | 0/2000 [00:00<?, ?it/s]

Running chain 3:   0%|                                                                 | 0/2000 [00:01<?, ?it/s]

Running chain 0:   0%|                                                                 | 0/2000 [00:01<?, ?it/s]

Running chain 1:   0%|                                                                 | 0/2000 [00:01<?, ?it/s]

Running chain 2:   0%|                                                                 | 0/2000 [00:01<?, ?it/s]

Running chain 0: 100%|████████████████████████████████████████████████████| 2000/2000 [00:01<00:00, 1249.77it/s]

Running chain 1: 100%|████████████████████████████████████████████████████| 2000/2000 [00:01<00:00, 1250.77it/s]

Running chain 2: 100%|████████████████████████████████████████████████████| 2000/2000 [00:01<00:00, 1251.66it/s]

Running chain 3: 100%|████████████████████████████████████████████████████| 2000/2000 [00:01<00:00, 1252.61it/s]

Sampling time = 0:00:01.844539

Transforming variables...

Transformation time = 0:00:00.108068

Distribuzioni a posteriori

Esaminiamo le distribuzioni a posteriori e le diagnostiche del campionamento.

az.plot_trace(idata, combined=True, divergences="bottom", kind="rank_vlines")
plt.tight_layout();

/var/folders/cl/wwjrsxdd5tz7y9jr82nd5hrw0000gn/T/ipykernel_78954/1817941948.py:2: UserWarning: The figure layout has changed to tight
  plt.tight_layout();

Abbiamo ottenuto due distribuzioni a posteriori per ciascuno dei due parametri \(\mu\) e \(\sigma\) per ciascun grouppo.

Un sommario numerico delle distribuzioni a posteriori si ottiene nel modo seguente.

az.summary(idata, var_names=["group1_mean", "group0_mean", "group1_std", "group0_std"])

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
group1_mean	89.312	1.008	87.413	91.191	0.014	0.010	5061.0	3648.0	1.0
group0_mean	77.629	2.308	73.372	81.871	0.036	0.025	4206.0	3016.0	1.0
group1_std	19.103	0.726	17.779	20.426	0.012	0.008	3857.0	2946.0	1.0
group0_std	22.820	1.750	19.582	25.991	0.027	0.020	4265.0	2853.0	1.0

Il nostro obiettivo è comprendere se le medie dei due gruppi sono diverse, e l’incertezza associata alla stima a posteriori del parametro diff_of_means è fondamentale per rispondere a questa domanda. Se l’intervallo di credibilità associato a diff_of_means non include lo 0, allora possiamo concludere con un certo grado di sicurezza che le medie dei due gruppi sono diverse. In altre parole, se l’intervallo di credibilità non contiene lo 0, allora ci sono prove convincenti che le medie dei due gruppi sono diverse.

az.summary(idata, var_names=["difference of means", "difference of stds", "effect size"])

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
difference of means	11.683	2.506	7.145	16.443	0.037	0.026	4534.0	3073.0	1.0
difference of stds	-3.716	1.888	-7.309	-0.339	0.029	0.022	4198.0	2912.0	1.0
effect size	0.556	0.120	0.341	0.793	0.002	0.001	4667.0	3111.0	1.0

Nel caso presente, l’intervallo di credibilità al 94% non include lo 0. Pertanto, possiamo concludere, con un livello di sicurezza soggettivo del 94%, che il QI dei bambini le cui madri hanno completato le scuole superiori tende ad essere più elevato rispetto a quello dei bambini le cui madri non hanno completato le scuole superiori.

az.plot_posterior(idata, var_names="difference of means", ref_val=0, figsize=(6, 3));

La figura seguente mostra la distribuzione a posteriori della grandezza dell’effetto.

az.plot_posterior(idata, var_names="effect size", ref_val=0, figsize=(6, 3));

Possiamo dunque concludere che, per ciò che concerne l’effetto della scolarità della madre sul quoziente di intelligenza del bambino, la dimensione dell’effetto è “media”.

Verifica di ipotesi bayesiana

Come “bonus” di questa analisi statistica, possiamo esaminare l’approccio bayesiano equivalente al test di ipotesi tradizionale.

Dopo aver ottenuto un campione dalla distribuzione a posteriori del parametro di interesse \(\mu\) per ciascun gruppo, possiamo porci la domanda: qual è la probabilità che il QI di un bambino in un gruppo sia maggiore di quello di un bambino nell’altro gruppo? Per rispondere a questa domanda, useremo campioni casuali dalle distribuzioni a posteriori dei parametri. Dovremo confrontare le coppie di valori campionati dalle due distribuzioni a posteriori del parametro di interesse e calcolare la media di tali confronti. Questo ci fornirà un’indicazione sulla probabilità che il QI di un bambino in un gruppo sia maggiore di quello di un bambino nell’altro gruppo, basandoci sulla distribuzione a posteriori dei parametri stimati dal modello.

Per recuperare i valori numerici della distribuzione a posteriori del parametro \(\mu\) dall’oggetto idata, possiamo utilizzare le funzioni di ArviZ. Esaminiamo le proprietà dell’oggetto idata.

idata

arviz.InferenceData

• posterior

  <xarray.Dataset>
  Dimensions:              (chain: 4, draw: 1000)
  Coordinates:
    * chain                (chain) int64 0 1 2 3
    * draw                 (draw) int64 0 1 2 3 4 5 6 ... 994 995 996 997 998 999
  Data variables:
      group1_mean          (chain, draw) float64 89.48 89.94 87.98 ... 89.02 89.51
      group0_mean          (chain, draw) float64 78.64 79.61 78.58 ... 74.19 76.95
      group1_std           (chain, draw) float64 18.85 19.7 19.42 ... 19.63 20.5
      group0_std           (chain, draw) float64 21.67 24.27 22.53 ... 23.27 22.86
      difference of means  (chain, draw) float64 10.84 10.33 9.393 ... 14.83 12.56
      difference of stds   (chain, draw) float64 -2.824 -4.57 ... -3.641 -2.357
      effect size          (chain, draw) float64 0.5337 0.4674 ... 0.6889 0.5785
  Attributes:
      created_at:     2024-01-26T17:59:51.448577
      arviz_version:  0.17.0

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 1000
  • Coordinates: (2)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 5 ... 995 996 997 998 999

      array([  0,   1,   2, ..., 997, 998, 999])

  • Data variables: (7)
    • group1_mean
      (chain, draw)
      float64
      89.48 89.94 87.98 ... 89.02 89.51

      array([[89.47993858, 89.93599918, 87.97560246, ..., 87.04364353,
              91.4285284 , 87.84108104],
             [89.87818295, 88.76829027, 90.2915316 , ..., 90.27555962,
              89.33418335, 88.61000362],
             [89.2117991 , 89.03178963, 88.42159145, ..., 88.32004416,
              88.31587056, 88.49490234],
             [88.31011802, 90.20199239, 88.26128404, ..., 89.43541314,
              89.01814444, 89.50647438]])

    • group0_mean
      (chain, draw)
      float64
      78.64 79.61 78.58 ... 74.19 76.95

      array([[78.64154323, 79.60616406, 78.58285897, ..., 77.11758731,
              76.67104164, 75.73123668],
             [81.76121427, 73.62522105, 81.37889455, ..., 76.87532534,
              77.47445558, 77.82979318],
             [71.100319  , 76.00777222, 78.43372643, ..., 77.91002655,
              82.23866895, 82.85562885],
             [80.02085046, 75.19408973, 77.5560896 , ..., 78.33432467,
              74.1864163 , 76.94711829]])

    • group1_std
      (chain, draw)
      float64
      18.85 19.7 19.42 ... 19.63 20.5

      array([[18.84835455, 19.69704983, 19.41918455, ..., 20.01921524,
              18.65785299, 19.15483946],
             [18.86504218, 19.0738071 , 19.02604428, ..., 19.97325902,
              19.36606096, 18.18683247],
             [18.86093128, 19.76545932, 19.01293661, ..., 18.69863395,
              18.64782274, 18.07374153],
             [19.85109875, 18.33233699, 19.37432035, ..., 18.65387814,
              19.63316695, 20.49899775]])

    • group0_std
      (chain, draw)
      float64
      21.67 24.27 22.53 ... 23.27 22.86

      array([[21.67239295, 24.2666873 , 22.5347102 , ..., 22.51954619,
              22.63627358, 23.8867639 ],
             [22.95116938, 21.96071642, 24.08381751, ..., 21.93756423,
              22.2248329 , 20.95865748],
             [24.35118461, 22.52929911, 22.24160556, ..., 25.47075698,
              27.4969653 , 27.52663764],
             [21.3096265 , 23.86557297, 22.06471522, ..., 24.19115855,
              23.27417464, 22.85612764]])

    • difference of means
      (chain, draw)
      float64
      10.84 10.33 9.393 ... 14.83 12.56

      array([[10.83839535, 10.32983512,  9.39274349, ...,  9.92605621,
              14.75748675, 12.10984436],
             [ 8.11696868, 15.14306922,  8.91263704, ..., 13.40023428,
              11.85972777, 10.78021044],
             [18.1114801 , 13.02401741,  9.98786502, ..., 10.41001761,
               6.0772016 ,  5.63927349],
             [ 8.28926756, 15.00790266, 10.70519444, ..., 11.10108846,
              14.83172814, 12.55935608]])

    • difference of stds
      (chain, draw)
      float64
      -2.824 -4.57 ... -3.641 -2.357

      array([[-2.8240384 , -4.56963747, -3.11552565, ..., -2.50033095,
              -3.97842059, -4.73192445],
             [-4.0861272 , -2.88690932, -5.05777323, ..., -1.96430521,
              -2.85877194, -2.77182501],
             [-5.49025333, -2.76383979, -3.22866895, ..., -6.77212303,
              -8.84914257, -9.45289611],
             [-1.45852775, -5.53323598, -2.69039487, ..., -5.53728041,
              -3.6410077 , -2.35712989]])

    • effect size
      (chain, draw)
      float64
      0.5337 0.4674 ... 0.6889 0.5785

      array([[0.53366087, 0.46740716, 0.44653548, ..., 0.46587883, 0.71145565,
              0.55933412],
             [0.38638085, 0.73624503, 0.41066814, ..., 0.63876295, 0.56896159,
              0.54940105],
             [0.83157433, 0.61455834, 0.48273071, ..., 0.46592341, 0.25868345,
              0.24218555],
             [0.40252295, 0.70527285, 0.51558649, ..., 0.51392295, 0.68886182,
              0.57851685]])

  • Indexes: (2)
    • chain
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3], dtype='int64', name='chain'))

    • draw
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             990, 991, 992, 993, 994, 995, 996, 997, 998, 999],
            dtype='int64', name='draw', length=1000))

  • Attributes: (2)

    created_at :

    2024-01-26T17:59:51.448577

    arviz_version :

    0.17.0

  
  
• sample_stats

  <xarray.Dataset>
  Dimensions:          (chain: 4, draw: 1000)
  Coordinates:
    * chain            (chain) int64 0 1 2 3
    * draw             (draw) int64 0 1 2 3 4 5 6 ... 993 994 995 996 997 998 999
  Data variables:
      acceptance_rate  (chain, draw) float64 0.9811 0.9027 ... 0.9055 0.9627
      step_size        (chain, draw) float64 0.6869 0.6869 ... 0.7466 0.7466
      diverging        (chain, draw) bool False False False ... False False False
      energy           (chain, draw) float64 1.922e+03 1.922e+03 ... 1.923e+03
      n_steps          (chain, draw) int64 7 3 7 3 3 7 7 7 7 ... 7 3 3 7 3 7 3 7 7
      tree_depth       (chain, draw) int64 3 2 3 2 2 3 3 3 3 ... 3 2 2 3 2 3 2 3 3
      lp               (chain, draw) float64 1.92e+03 1.921e+03 ... 1.922e+03
  Attributes:
      created_at:     2024-01-26T17:59:51.451453
      arviz_version:  0.17.0

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 1000
  • Coordinates: (2)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 5 ... 995 996 997 998 999

      array([  0,   1,   2, ..., 997, 998, 999])

  • Data variables: (7)
    • acceptance_rate
      (chain, draw)
      float64
      0.9811 0.9027 ... 0.9055 0.9627

      array([[0.98111554, 0.90272499, 0.98622356, ..., 0.82520578, 0.99590741,
              0.9481697 ],
             [0.92138161, 0.97176364, 0.97409675, ..., 0.9897707 , 1.        ,
              0.86580464],
             [0.57201637, 1.        , 0.71971671, ..., 0.85838273, 0.72092925,
              0.97441452],
             [0.84401806, 1.        , 0.98174208, ..., 0.87430664, 0.90554357,
              0.96272327]])

    • step_size
      (chain, draw)
      float64
      0.6869 0.6869 ... 0.7466 0.7466

      array([[0.68688308, 0.68688308, 0.68688308, ..., 0.68688308, 0.68688308,
              0.68688308],
             [0.69021349, 0.69021349, 0.69021349, ..., 0.69021349, 0.69021349,
              0.69021349],
             [0.80395877, 0.80395877, 0.80395877, ..., 0.80395877, 0.80395877,
              0.80395877],
             [0.74660382, 0.74660382, 0.74660382, ..., 0.74660382, 0.74660382,
              0.74660382]])

    • diverging
      (chain, draw)
      bool
      False False False ... False False

      array([[False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False]])

    • energy
      (chain, draw)
      float64
      1.922e+03 1.922e+03 ... 1.923e+03

      array([[1921.84780935, 1921.72580419, 1922.5103851 , ..., 1923.52603467,
              1924.55900393, 1923.96039287],
             [1922.99894355, 1922.88777587, 1923.66669621, ..., 1923.29335416,
              1921.1762929 , 1921.999636  ],
             [1928.42010591, 1924.22147444, 1924.13967984, ..., 1922.46848121,
              1926.98890382, 1926.52684919],
             [1922.8927189 , 1922.21488175, 1922.5760905 , ..., 1925.26086374,
              1922.46218974, 1922.98590566]])

    • n_steps
      (chain, draw)
      int64
      7 3 7 3 3 7 7 7 ... 3 3 7 3 7 3 7 7

      array([[7, 3, 7, ..., 7, 7, 7],
             [7, 7, 7, ..., 7, 3, 3],
             [3, 3, 3, ..., 7, 3, 7],
             [3, 7, 7, ..., 3, 7, 7]])

    • tree_depth
      (chain, draw)
      int64
      3 2 3 2 2 3 3 3 ... 2 2 3 2 3 2 3 3

      array([[3, 2, 3, ..., 3, 3, 3],
             [3, 3, 3, ..., 3, 2, 2],
             [2, 2, 2, ..., 3, 2, 3],
             [2, 3, 3, ..., 2, 3, 3]])

    • lp
      (chain, draw)
      float64
      1.92e+03 1.921e+03 ... 1.922e+03

      array([[1920.26351241, 1921.36292055, 1921.0149726 , ..., 1923.05744828,
              1922.42268354, 1921.65845682],
             [1921.70360617, 1921.78087669, 1921.9712219 , ..., 1921.27580506,
              1920.10955314, 1921.4619374 ],
             [1924.04958648, 1920.75763356, 1920.43466424, ..., 1921.92021646,
              1925.21430706, 1926.27865133],
             [1921.87110579, 1921.71870831, 1920.63643283, ..., 1920.64690082,
              1921.52329917, 1921.86660801]])

  • Indexes: (2)
    • chain
      PandasIndex

      PandasIndex(Index([0, 1, 2, 3], dtype='int64', name='chain'))

    • draw
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             990, 991, 992, 993, 994, 995, 996, 997, 998, 999],
            dtype='int64', name='draw', length=1000))

  • Attributes: (2)

    created_at :

    2024-01-26T17:59:51.451453

    arviz_version :

    0.17.0

  
  
• observed_data

  <xarray.Dataset>
  Dimensions:         (mom_hs_1_dim_0: 341, mom_hs_0_dim_0: 93)
  Coordinates:
    * mom_hs_1_dim_0  (mom_hs_1_dim_0) int64 0 1 2 3 4 5 ... 336 337 338 339 340
    * mom_hs_0_dim_0  (mom_hs_0_dim_0) int64 0 1 2 3 4 5 6 ... 87 88 89 90 91 92
  Data variables:
      mom_hs_1        (mom_hs_1_dim_0) float64 65.0 98.0 85.0 ... 76.0 88.0 70.0
      mom_hs_0        (mom_hs_0_dim_0) float64 98.0 102.0 101.0 ... 93.0 94.0 50.0
  Attributes:
      created_at:                 2024-01-26T17:59:51.452606
      arviz_version:              0.17.0
      inference_library:          numpyro
      inference_library_version:  0.13.2
      sampling_time:              1.844539

  xarray.Dataset

  • Dimensions:
    • mom_hs_1_dim_0: 341
    • mom_hs_0_dim_0: 93
  • Coordinates: (2)
    • mom_hs_1_dim_0
      (mom_hs_1_dim_0)
      int64
      0 1 2 3 4 5 ... 336 337 338 339 340

      array([  0,   1,   2, ..., 338, 339, 340])

    • mom_hs_0_dim_0
      (mom_hs_0_dim_0)
      int64
      0 1 2 3 4 5 6 ... 87 88 89 90 91 92

      array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
             18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
             36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53,
             54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71,
             72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89,
             90, 91, 92])

  • Data variables: (2)
    • mom_hs_1
      (mom_hs_1_dim_0)
      float64
      65.0 98.0 85.0 ... 76.0 88.0 70.0

      array([ 65.,  98.,  85.,  83., 115.,  69., 106., 102.,  95.,  91.,  58.,
              84.,  78., 110., 102.,  99., 105., 102., 115., 100.,  87.,  96.,
              72.,  78.,  77.,  98.,  69., 130., 109.,  92., 100., 107.,  86.,
              90., 110., 107., 113.,  65., 102., 111.,  42., 100.,  67.,  92.,
             100., 110.,  56., 107.,  97.,  95.,  86.,  79.,  81.,  79.,  56.,
              52.,  63.,  80.,  87.,  88.,  92., 117., 102., 107.,  99.,  56.,
             110., 109.,  86.,  92.,  91., 123., 102., 105., 114.,  96.,  66.,
             104., 108.,  84.,  83.,  93., 114., 106.,  87.,  65.,  95.,  61.,
              73., 112., 113., 105., 122.,  96.,  97.,  94., 117., 136., 116.,
             106.,  99.,  94.,  89., 119., 112.,  92.,  86.,  69.,  45.,  57.,
              94., 104., 144.,  52., 102., 106.,  98.,  97.,  94., 111., 100.,
             105.,  90.,  98., 121., 106., 121., 102.,  64.,  99.,  81.,  69.,
              84., 104., 104., 107.,  88.,  94., 109.,  94.,  98., 102., 104.,
              87., 102.,  77., 109.,  94.,  86.,  97.,  97.,  88., 103.,  87.,
              90.,  65., 109.,  87.,  58.,  87., 113.,  75.,  99., 108.,  95.,
             100.,  85., 108.,  90., 100.,  94.,  95., 119.,  98., 112., 122.,
             126., 116.,  98.,  94.,  93.,  90.,  70., 110., 104.,  83.,  99.,
              81., 104., 109., 113.,  74.,  81.,  89.,  93., 102.,  95.,  85.,
              97.,  92., 104., 120., 105.,  68., 104.,  80., 120.,  94.,  81.,
             101.,  61.,  68., 110., 113.,  86.,  83., 106., 104.,  78.,  99.,
              91.,  69.,  84.,  58.,  42.,  72.,  52., 101.,  63.,  73.,  68.,
              60.,  69.,  56.,  49.,  44.,  74.,  46.,  58.,  89.,  94.,  88.,
              79.,  87.,  46.,  95.,  92.,  62.,  52., 101.,  97.,  85.,  98.,
              94.,  90.,  72.,  92.,  83.,  64., 101.,  82.,  77., 101.,  50.,
              90., 103.,  96.,  50.,  47.,  73.,  77.,  64.,  52.,  61.,  86.,
              83.,  64.,  83., 116.,  74.,  76.,  92.,  98.,  96.,  67.,  84.,
             111.,  68., 107.,  89.,  83.,  73.,  74.,  94.,  58.,  76.,  61.,
              38., 100.,  84.,  99.,  86.,  94.,  90., 112.,  58.,  87.,  76.,
              68., 110.,  56.,  79.,  82.,  60., 102.,  87., 119.,  84.,  86.,
             113.,  72., 104.,  94., 104.,  64.,  95.,  56.,  91., 106.,  76.,
              90., 108.,  85.,  87.,  89.,  42., 102., 104.,  76.,  88.,  70.])

    • mom_hs_0
      (mom_hs_0_dim_0)
      float64
      98.0 102.0 101.0 ... 93.0 94.0 50.0

      array([ 98., 102., 101.,  99., 106., 103.,  56.,  78.,  76.,  79., 100.,
              94.,  73.,  78.,  94.,  83.,  92., 109.,  95.,  49.,  85., 104.,
              89.,  67., 103., 114.,  93.,  87., 111.,  64.,  78.,  97.,  95.,
              91.,  97.,  82., 100., 136.,  95.,  78.,  83.,  89.,  98.,  50.,
              57., 106.,  40.,  42.,  80.,  58.,  73.,  75.,  20.,  71.,  46.,
              54.,  54.,  58.,  76.,  43.,  60.,  43.,  42.,  75.,  62.,  41.,
             100.,  42.,  41.,  82.,  50.,  88.,  87.,  54.,  49.,  80.,  42.,
              78.,  80.,  67.,  96.,  65.,  75.,  86., 104.,  41., 106.,  76.,
             100.,  59.,  93.,  94.,  50.])

  • Indexes: (2)
    • mom_hs_1_dim_0
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             331, 332, 333, 334, 335, 336, 337, 338, 339, 340],
            dtype='int64', name='mom_hs_1_dim_0', length=341))

    • mom_hs_0_dim_0
      PandasIndex

      PandasIndex(Index([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
             18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
             36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53,
             54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71,
             72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89,
             90, 91, 92],
            dtype='int64', name='mom_hs_0_dim_0'))

  • Attributes: (5)

    created_at :

    2024-01-26T17:59:51.452606

    arviz_version :

    0.17.0

    inference_library :

    numpyro

    inference_library_version :

    0.13.2

    sampling_time :

    1.844539

  
  

Recuperiamo i valori a posteriori della media di ciascuno dei due gruppi.

mu1_post = idata.posterior["group1_mean"]
mu0_post = idata.posterior["group0_mean"]

Per ciacun gruppo, l’oggetto ottenuto è un array di dimensioni \(4 \times 1000\).

mu1_post.shape

(4, 1000)

Una volta ottenute le stime a posteriori dei parametri di interesse per ciascun gruppo, possiamo procedere con il test di ipotesi bayesiano. Ad esempio, ci possiamo chiedere quale sia la probabilità che un bambino la cui madre non ha completato la scuola superiore abbia un QI maggiore di un bambino la cui madre ha completato la scuola superiore.

print(np.mean(mu1_post < mu0_post))

<xarray.DataArray ()>
array(0.)

Una tale probabilità è effettivamente uguale a zero il che conferma il risultato precedente, ovvero l’iportanza del livello di istruzione della madre per il QI del figlio.

Commenti e considerazioni finali

In questo capitolo abbiamo esaminato la procedura bayesiana per calcolare la distribuzione a posteriori della differenza tra le medie di due gruppi indipendenti. Inoltre, abbiamo esplorato il calcolo della dimensione dell’effetto in termini bayesiani. Nell’esempio trattato, abbiamo considerato il caso in cui la verosimiglianza è descritta da una distribuzione Gaussiana. Tuttavia, va sottolineato che la scelta di una distribuzione specifica per la verosimiglianza non è vincolante nella statistica bayesiana. È possibile utilizzare qualsiasi distribuzione di probabilità, purché sia adeguata ai dati del campione.

Nel caso del confronto tra le medie di due gruppi indipendenti, una distribuzione molto utilizzata è la distribuzione \(t\) di Student. Questa distribuzione è particolarmente vantaggiosa quando si desidera condurre un’analisi statistica “robusta”, ovvero un’analisi che non sia influenzata da osservazioni anomale o outlier presenti nei dati. Per questo motivo, la distribuzione \(t\) di Student è spesso preferita quando si lavora con dati che potrebbero contenere valori anomali.

Per eseguire un’analisi bayesiana utilizzando la distribuzione \(t\) di Student per il confronto tra le medie di due gruppi indipendenti, si può fare riferimento al sito di PyMC, dove è disponibile una descrizione dettagliata del modello bayesiano da utilizzare.

%run ../wm.py

Watermark:
----------
Last updated: 2024-01-26T18:59:52.627569+01:00

Python implementation: CPython
Python version       : 3.11.7
IPython version      : 8.19.0

Compiler    : Clang 16.0.6 
OS          : Darwin
Release     : 23.3.0
Machine     : x86_64
Processor   : i386
CPU cores   : 8
Architecture: 64bit