Confronto tra le medie di due gruppi#
Nel capitolo Confronto tra due gruppi, abbiamo esaminato come effettuare l’inferenza sulla differenza tra le medie di due campioni indipendenti usando un modello bayesiano. In questo approccio, trattiamo i due gruppi come entità separate e calcoliamo direttamente la differenza tra le loro medie. Per fare ciò, utilizziamo distribuzioni a priori per i parametri del modello, come le deviazioni standard e le medie, e poi aggiorniamo queste distribuzioni a posteriori utilizzando i dati osservati.
Tuttavia, possiamo ottenere lo stesso risultato utilizzando un modello di regressione. Invece di calcolare direttamente la differenza tra le medie dei due gruppi, includiamo una variabile “dummy” nel modello di regressione. Questa variabile “dummy” ci permette di rappresentare in modo binario l’appartenenza dei dati ai gruppi: assegniamo il valore 0 a un gruppo di riferimento e il valore 1 al gruppo di confronto. Successivamente, il modello di regressione stima un parametro di regressione associato alla variabile “dummy”, che rappresenta proprio la differenza tra le medie dei due gruppi. Quindi, la variabile “dummy” funge da indicatore per il gruppo e ci consente di ottenere la stima della differenza tra le medie dei due gruppi in modo efficiente.
Entrambi gli approcci sono validi per l’inferenza sulla differenza tra le medie dei due gruppi indipendenti, ma il modello di regressione offre una maggiore flessibilità e possibilità di espansione. Infatti, possiamo includere ulteriori variabili esplicative nel modello di regressione per comprendere meglio i fattori che influenzano il risultato di interesse. Questa flessibilità è particolarmente utile quando vogliamo esplorare come altre variabili possano contribuire alla differenza tra le medie dei gruppi o quando desideriamo considerare più fattori simultaneamente.
Preparazione del Notebook#
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import scipy as sc
import statistics as st
import arviz as az
import bambi as bmb
import pymc as pm
import pymc.sampling_jax
from pymc import HalfNormal, Model, Normal, sample
import warnings
warnings.filterwarnings("ignore", category=UserWarning)
warnings.filterwarnings("ignore", category=FutureWarning)
warnings.filterwarnings("ignore", category=Warning)
/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")
sns.set_theme(palette="colorblind")
Regressione bayesiana per due gruppi indipendenti#
Nel contesto bayesiano, il modello di regressione può essere formulato nel modo seguente:
In questa rappresentazione:
\( \alpha \) agisce come intercetta,
\( \beta \) è il coefficiente angolare o la pendenza,
\( \sigma \) è l’errore standard associato alle osservazioni.
Nel caso specifico, la variabile \( x \) è una variabile indicatrice che assume i valori 0 o 1. Per il gruppo identificato da \( x = 0 \), il modello si riduce a:
Questo implica che \( \alpha \) rappresenta la media del gruppo codificato come \( x = 0 \).
Per il gruppo contrassegnato da \( x = 1 \), il modello diventa:
In termini dei parametri del modello, la media per il gruppo codificato con \( x = 1 \) è rappresentata da \( \alpha + \beta \). In questa configurazione, \( \beta \) indica la differenza tra la media del gruppo con \( x = 1 \) e quella del gruppo con \( x = 0 \). Di conseguenza, l’analisi della differenza tra le medie dei due gruppi può essere effettuata attraverso l’inferenza sul parametro \( \beta \). In sintesi, per confrontare le medie dei due gruppi indipendenti, si può esaminare la distribuzione a posteriori di \( \beta \).
Un esempio illustrativo#
Esaminiamo nuovamente i dati relativi al quoziente di intelligenza dei bambini le cui madri hanno completato oppure no la scuola superiore. Ci poniamo il problema di replicare i risultati ottenuti in precedenza usando l’analisi di regressione.
Leggiamo i dati:
kidiq = pd.read_stata("../data/kidiq.dta")
kidiq.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 |
kidiq.groupby(["mom_hs"]).size()
mom_hs
0.0 93
1.0 341
dtype: int64
Ci sono 93 bambini la cui madre non ha completato le superiori e 341 bambini la cui madre ha ottenuto il diploma di scuola superiore.
summary_stats = [st.mean, st.stdev]
kidiq.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 |
az.plot_violin(
{
"mom_hs=0": kidiq.loc[kidiq.mom_hs == 0, "kid_score"],
"mom_hs=1": kidiq.loc[kidiq.mom_hs == 1, "kid_score"],
}
);
Iniziamo l’inferenza statistica sulla differenza tra le medie dei due gruppi utilizzando bambi. Questo pacchetto offre una sintassi semplice per formulare il modello bayesiano di interesse. Un altro vantaggio è che bambi selezionerà automaticamente le distribuzioni a priori appropriate per i parametri del modello, rendendo il processo più intuitivo.
Il modello di regressione sopra descritto si scrive nel modo seguente.
mod = bmb.Model("kid_score ~ mom_hs", kidiq)
Effettuiamo il campionamento.
results = mod.fit(method="nuts_numpyro", idata_kwargs={"log_likelihood": True})
Show code cell output
Compiling...
Compilation time = 0:00:02.737080
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 0: 0%| | 0/2000 [00:02<?, ?it/s]
Running chain 3: 0%| | 0/2000 [00:02<?, ?it/s]
Running chain 2: 0%| | 0/2000 [00:02<?, ?it/s]
Running chain 1: 0%| | 0/2000 [00:02<?, ?it/s]
Running chain 0: 100%|█████████████████████████████████████████████████████| 2000/2000 [00:02<00:00, 870.55it/s]
Running chain 1: 100%|█████████████████████████████████████████████████████| 2000/2000 [00:02<00:00, 871.33it/s]
Running chain 2: 100%|█████████████████████████████████████████████████████| 2000/2000 [00:02<00:00, 872.08it/s]
Running chain 3: 100%|█████████████████████████████████████████████████████| 2000/2000 [00:02<00:00, 872.85it/s]
Sampling time = 0:00:02.766127
Transforming variables...
Transformation time = 0:00:00.235790
Computing Log Likelihood...
Log Likelihood time = 0:00:00.200011
Possiamo ispezionare le proprietà del modello nel modo seguente.
mod
Formula: kid_score ~ mom_hs
Family: gaussian
Link: mu = identity
Observations: 434
Priors:
target = mu
Common-level effects
Intercept ~ Normal(mu: 86.7972, sigma: 110.1032)
mom_hs ~ Normal(mu: 0.0, sigma: 124.2132)
Auxiliary parameters
sigma ~ HalfStudentT(nu: 4.0, sigma: 20.3872)
------
* To see a plot of the priors call the .plot_priors() method.
* To see a summary or plot of the posterior pass the object returned by .fit() to az.summary() or az.plot_trace()
Le distribuzioni a priori utilizzate di default dal modello possono essere visualizzate nel modo seguente.
mod.plot_priors();
Sampling: [Intercept, kid_score_sigma, mom_hs]
Per ispezionare il nostro posteriore e il processo di campionamento possiamo utilizzare az.plot_trace(). L’opzione kind='rank_vlines' ci fornisce una variante del grafico di rango che utilizza linee e punti e ci aiuta a ispezionare la stazionarietà delle catene. Poiché non c’è un modello chiaro o deviazioni serie dalle linee orizzontali, possiamo concludere che le catene sono stazionarie.
az.plot_trace(results, kind="rank_vlines")
plt.tight_layout()
az.summary(results, round_to=2)
| mean | sd | hdi_3% | hdi_97% | mcse_mean | mcse_sd | ess_bulk | ess_tail | r_hat | |
|---|---|---|---|---|---|---|---|---|---|
| Intercept | 77.52 | 2.08 | 73.91 | 81.67 | 0.03 | 0.02 | 4953.15 | 3053.09 | 1.0 |
| mom_hs | 11.81 | 2.33 | 7.09 | 15.82 | 0.03 | 0.02 | 4792.95 | 3193.39 | 1.0 |
| kid_score_sigma | 19.89 | 0.68 | 18.64 | 21.16 | 0.01 | 0.01 | 4761.95 | 2926.08 | 1.0 |
Il parametro “Intercept” rappresenta la stima a posteriori del punteggio del QI per il gruppo codificato con “mom_hs” uguale a 0. La media a posteriori di questo gruppo è di 77.6, che è praticamente identica al valore campionario corrispondente.
Il parametro “mom_hs” corrisponde alla stima a posteriori della differenza nei punteggi del QI tra il gruppo codificato con “mom_hs” uguale a 1 e il gruppo codificato con “mom_hs” uguale a 0. Anche in questo caso, la differenza a posteriori di 11.8 tra le medie dei due gruppi è molto simile alla differenza campionaria tra le medie dei due gruppi. La parte importante della tabella riguarda l’intervallo di credibilità al 94%, che è [7.5, 16.2], e che non include lo 0. Ciò significa che, con un livello di certezza soggettiva del 94%, possiamo essere sicuri che il QI dei bambini le cui madri hanno il diploma superiore sarà maggiore (in media) di almeno 7.5 punti, e tale differenza può arrivare fino a 16.2 punti, rispetto al QI dei bambini le cui madri non hanno completato la scuola superiore.
Se confrontiamo questi risultati con quelli ottenuti nel capitolo Confronto tra due gruppi, notiamo che sono quasi identici. Le piccole differenze che si osservano possono essere attribuite sia all’approssimazione numerica sia al fatto che nel modello precedente abbiamo consentito deviazioni standard diverse per i due gruppi, mentre nel caso attuale abbiamo assumo la stessa variabilità per entrambi i gruppi.
Il test bayesiano di ipotesi può essere svolto, per esempio, calcolando la probabilità che \(\beta_{mean\_diff} > 0\). Questa probabilità è 1, per cui concludiamo che la media del gruppo codificato con “mom_hs = 1” è maggiore della media del gruppo codificato con “mom_hs = 0”.
az.plot_posterior(results, var_names="mom_hs", ref_val=0, figsize=(6, 3));
Un valore numerico si ottiene nel modo seguente.
results.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 ... 993 994 995 996 997 998 999
Data variables:
Intercept (chain, draw) float64 77.67 76.74 79.07 ... 76.6 80.76
mom_hs (chain, draw) float64 11.13 12.84 11.07 ... 12.5 12.77 8.98
kid_score_sigma (chain, draw) float64 20.59 21.43 19.26 ... 20.92 18.91
Attributes:
created_at: 2024-01-26T21:45:38.422075
arviz_version: 0.17.0
modeling_interface: bambi
modeling_interface_version: 0.13.0# Probabiliy that posterior is > 0
(results.posterior["mom_hs"] > 0).mean().item()
1.0
Parametrizzazione alternativa#
Consideriamo adesso il caso in cui, per distinguere i gruppi, anziché una variabile dicotomica, con valori 0 e 1, usiamo una variabile qualitativa con i nomi dei due gruppi. Introduciamo questa nuova variabile nel data frame.
# Add a new column 'hs' with the categories based on 'mom_hs'
kidiq["hs"] = kidiq["mom_hs"].map({0: "not_completed", 1: "completed"})
kidiq.tail()
| kid_score | mom_hs | mom_iq | mom_work | mom_age | hs | |
|---|---|---|---|---|---|---|
| 429 | 94 | 0.0 | 84.877412 | 4 | 21 | not_completed |
| 430 | 76 | 1.0 | 92.990392 | 4 | 23 | completed |
| 431 | 50 | 0.0 | 94.859708 | 2 | 24 | not_completed |
| 432 | 88 | 1.0 | 96.856624 | 2 | 21 | completed |
| 433 | 70 | 1.0 | 91.253336 | 2 | 25 | completed |
Adattiamo il modello ai dati, usando questa nuova variabile e forziamo a zero l’intercetta che Bambi aggiunge di default al modello.
mod_2 = bmb.Model("kid_score ~ 0 + hs", kidiq)
results_2 = mod_2.fit(method="nuts_numpyro", idata_kwargs={"log_likelihood": True})
Show code cell output
Compiling...
Compilation time = 0:00:00.844433
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:02<?, ?it/s]
Running chain 0: 0%| | 0/2000 [00:02<?, ?it/s]
Running chain 1: 0%| | 0/2000 [00:02<?, ?it/s]
Running chain 2: 0%| | 0/2000 [00:02<?, ?it/s]
Running chain 0: 100%|█████████████████████████████████████████████████████| 2000/2000 [00:02<00:00, 943.46it/s]
Running chain 1: 100%|█████████████████████████████████████████████████████| 2000/2000 [00:02<00:00, 943.97it/s]
Running chain 2: 100%|█████████████████████████████████████████████████████| 2000/2000 [00:02<00:00, 944.66it/s]
Running chain 3: 100%|█████████████████████████████████████████████████████| 2000/2000 [00:02<00:00, 945.28it/s]
Sampling time = 0:00:02.261408
Transforming variables...
Transformation time = 0:00:00.092107
Computing Log Likelihood...
Log Likelihood time = 0:00:00.150272
Ispezionare il modello e le distribuzioni a priori.
mod_2
Formula: kid_score ~ 0 + hs
Family: gaussian
Link: mu = identity
Observations: 434
Priors:
target = mu
Common-level effects
hs ~ Normal(mu: [0. 0.], sigma: [124.2132 124.2132])
Auxiliary parameters
sigma ~ HalfStudentT(nu: 4.0, sigma: 20.3872)
------
* To see a plot of the priors call the .plot_priors() method.
* To see a summary or plot of the posterior pass the object returned by .fit() to az.summary() or az.plot_trace()
mod_2.plot_priors();
Sampling: [hs, kid_score_sigma]
Esaminiamo le distribuzioni a posteriori dei parametri del modello.
az.plot_trace(results_2, kind="rank_vlines");
az.summary(results_2)
| mean | sd | hdi_3% | hdi_97% | mcse_mean | mcse_sd | ess_bulk | ess_tail | r_hat | |
|---|---|---|---|---|---|---|---|---|---|
| hs[completed] | 89.285 | 1.095 | 87.247 | 91.353 | 0.017 | 0.012 | 4335.0 | 2655.0 | 1.0 |
| hs[not_completed] | 77.540 | 2.086 | 73.689 | 81.467 | 0.033 | 0.024 | 3918.0 | 2543.0 | 1.0 |
| kid_score_sigma | 19.877 | 0.691 | 18.639 | 21.224 | 0.011 | 0.008 | 3707.0 | 2699.0 | 1.0 |
In questo caso, notiamo che abbiamo ottenuto le distribuzioni a posteriori per i parametri hs[completed] e hs[not_completed] che corrispondono alle medie dei due gruppi. Tali distribuzioni a posteriori illustrano direttamente l’incertezza sulla media dei due gruppi, alla luce della variabilità campionaria e delle nostre credenze a priori.
Possiamo svolgere il test bayesiano di ipotesi sulla differenza tra le due medie a posteriori nel modo seguente.
post_group = results_2.posterior["hs"]
diff = post_group.sel(hs_dim="completed") - post_group.sel(hs_dim="not_completed")
az.plot_posterior(diff, ref_val=0, figsize=(6, 3));
# Probabiliy that posterior is > 0
(post_group > 0).mean().item()
1.0
I risultati sono ovviamente identici a quelli trovati in precedenza.
Analisi con PyMC#
Per completezza, ripetiamo l’analsi specificando il modello bayesiano con PyMC.
μ_m = kidiq.kid_score.mean()
μ_s = kidiq.kid_score.std() * 3
with Model() as model:
# Priors
alpha = Normal("alpha", mu=μ_m, sigma=μ_s)
beta = Normal("beta", mu=0, sigma=10)
sigma = pm.HalfNormal("sigma", sigma=30)
# Expected value of outcome
mu = alpha + beta * kidiq["mom_hs"]
# Likelihood of observations
Y_obs = pm.Normal("Y_obs", mu=mu, sigma=sigma, observed=kidiq["kid_score"])
pm.model_to_graphviz(model)
Esaminiamo la distribuzione predittia a priori.
with model:
prior_samples = pm.sample_prior_predictive(1000)
Sampling: [Y_obs, alpha, beta, sigma]
az.plot_dist(prior_samples.prior_predictive["Y_obs"], figsize=(6, 3));
Eseguiamo il campionamento.
with model:
trace = pm.sampling_jax.sample_numpyro_nuts()
Show code cell output
Compiling...
Compilation time = 0:00:00.690465
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 2: 0%| | 0/2000 [00:02<?, ?it/s]
Running chain 1: 0%| | 0/2000 [00:02<?, ?it/s]
Running chain 0: 0%| | 0/2000 [00:02<?, ?it/s]
Running chain 3: 0%| | 0/2000 [00:02<?, ?it/s]
Running chain 3: 95%|████████████████████████████████████████████████▍ | 1900/2000 [00:02<00:00, 18977.21it/s]
Running chain 0: 100%|█████████████████████████████████████████████████████| 2000/2000 [00:02<00:00, 829.85it/s]
Running chain 1: 100%|█████████████████████████████████████████████████████| 2000/2000 [00:02<00:00, 830.22it/s]
Running chain 2: 100%|█████████████████████████████████████████████████████| 2000/2000 [00:02<00:00, 830.78it/s]
Running chain 3: 100%|█████████████████████████████████████████████████████| 2000/2000 [00:02<00:00, 831.27it/s]
Sampling time = 0:00:02.517849
Transforming variables...
Transformation time = 0:00:00.067944
Esaminiamo i risultati.
az.summary(trace, round_to=2)
| mean | sd | hdi_3% | hdi_97% | mcse_mean | mcse_sd | ess_bulk | ess_tail | r_hat | |
|---|---|---|---|---|---|---|---|---|---|
| alpha | 78.07 | 1.94 | 74.42 | 81.69 | 0.05 | 0.04 | 1481.65 | 1780.99 | 1.01 |
| beta | 11.12 | 2.19 | 6.95 | 15.15 | 0.06 | 0.04 | 1514.23 | 1714.81 | 1.01 |
| sigma | 19.90 | 0.69 | 18.65 | 21.23 | 0.02 | 0.01 | 2019.05 | 1843.36 | 1.00 |
I risultati replicano quelli ottenuti in precedenza.
Aggiungiamo al modello bayesiano la stima della grandezza dell’effetto.
with Model() as model:
# Priors
alpha = Normal("alpha", mu=μ_m, sigma=μ_s)
beta = Normal("beta", mu=0, sigma=10)
sigma = pm.HalfNormal("sigma", sigma=30)
# Expected value of outcome
mu = alpha + beta * kidiq["mom_hs"]
# Likelihood of observations
Y_obs = pm.Normal("Y_obs", mu=mu, sigma=sigma, observed=kidiq["kid_score"])
effect_size = pm.Deterministic("effect size", beta / sigma)
trace = pm.sampling_jax.sample_numpyro_nuts()
Show code cell output
Compiling...
Compilation time = 0:00:00.664939
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 2: 0%| | 0/2000 [00:02<?, ?it/s]
Running chain 0: 0%| | 0/2000 [00:02<?, ?it/s]
Running chain 1: 0%| | 0/2000 [00:02<?, ?it/s]
Running chain 3: 0%| | 0/2000 [00:02<?, ?it/s]
Running chain 0: 100%|█████████████████████████████████████████████████████| 2000/2000 [00:02<00:00, 910.29it/s]
Running chain 1: 100%|█████████████████████████████████████████████████████| 2000/2000 [00:02<00:00, 910.70it/s]
Running chain 2: 100%|█████████████████████████████████████████████████████| 2000/2000 [00:02<00:00, 911.31it/s]
Running chain 3: 100%|█████████████████████████████████████████████████████| 2000/2000 [00:02<00:00, 911.91it/s]
Sampling time = 0:00:02.299994
Transforming variables...
Transformation time = 0:00:00.073710
az.plot_posterior(trace, ref_val=0, var_names=["effect size"], figsize=(6, 3));
az.summary(trace, round_to=2)
| mean | sd | hdi_3% | hdi_97% | mcse_mean | mcse_sd | ess_bulk | ess_tail | r_hat | |
|---|---|---|---|---|---|---|---|---|---|
| alpha | 78.12 | 1.99 | 74.56 | 82.01 | 0.05 | 0.04 | 1428.76 | 1882.23 | 1.0 |
| beta | 11.09 | 2.24 | 6.74 | 15.16 | 0.06 | 0.04 | 1541.00 | 1972.10 | 1.0 |
| sigma | 19.88 | 0.69 | 18.63 | 21.17 | 0.02 | 0.01 | 2031.33 | 1786.49 | 1.0 |
| effect size | 0.56 | 0.11 | 0.35 | 0.78 | 0.00 | 0.00 | 1545.12 | 1835.72 | 1.0 |
Anche in questo caso i risultati replicano quelli ottenuti in precedenza.
%load_ext watermark
%watermark -n -u -v -iv -w
Last updated: Fri Jan 26 2024
Python implementation: CPython
Python version : 3.11.7
IPython version : 8.19.0
numpy : 1.26.2
seaborn : 0.13.0
pymc : 5.10.3
arviz : 0.17.0
pandas : 2.1.4
matplotlib: 3.8.2
scipy : 1.11.4
bambi : 0.13.0
Watermark: 2.4.3