here::here("code", "_common.R") |>
source()
# Load packages
if (!requireNamespace("pacman")) install.packages("pacman")
pacman::p_load(cmdstanr, posterior, bayestestR, brms, emmeans,
loo, bridgesampling)69 ANOVA ad due vie
In questo capitolo imparerai a
- fare inferenza sulla media di un campione;
- trovare le distribuzioni a posteriori usando
brms; - verificare il modello usando i pp-check plots.
Prerequisiti
- Leggere il capitolo Extending the Normal Regression Model del testo di Johnson et al. (2022).
Preparazione del Notebook
69.1 Introduzione
L’ANOVA a due o più vie estende il caso dell’ANOVA ad una via alla presenza di molteplici criteri di classificazione. Qui ci concentreremo sull’ANOVA a due vie.
69.1.1 Pattern delle medie nella classificazione a due vie
Immaginiamo di disporre delle medie di popolazione. La notazione per la classificazione a due vie è illustrata nella seguente tabella:
| C1 | C2 | … | Cc | Totale colonna | |
|---|---|---|---|---|---|
| R1 | µ11 | µ12 | … | µ1c | µ1: |
| R2 | µ21 | µ22 | … | µ2c | µ2: |
| … | … | … | … | … | … |
| Rr | µr1 | µr2 | … | µrc | µr: |
| ———— | ——- | ——- | —— | ——- | —————- |
| Totale riga | µ:1 | µ:2 | … | µ:c | µ:: |
Qui, i fattori \(R\) e \(C\) (così denominati in riferimento alle righe e alle colonne della tabella delle medie) presentano rispettivamente \(r\) e \(c\) categorie. Indichiamo le categorie dei fattori come \(R_j\) e \(C_k\).
All’interno di ogni cella del disegno sperimentale—cioè per ciascuna combinazione di categorie \(\{R_j, C_k\}\) dei due fattori—si trova una media di popolazione \(\mu_{jk}\) relativa alla variabile di risposta.
Per descrivere le medie su righe, colonne e quella complessiva, utilizziamo la notazione a “punti”:
\[ \mu_{j:} \;=\; \frac{\sum_{k=1}^{c} \mu_{jk}}{c} \quad\text{(media marginale sulla riga $j$)}, \]
\[ \mu_{:k} \;=\; \frac{\sum_{j=1}^{r} \mu_{jk}}{r} \quad\text{(media marginale sulla colonna $k$)}, \]
\[ \mu_{::} \;=\; \frac{\sum_{j=1}^{r}\sum_{k=1}^{c} \mu_{jk}}{r\,c} \;=\; \frac{\sum_{j=1}^{r} \mu_{j:}}{r} \;=\; \frac{\sum_{k=1}^{c} \mu_{:k}}{c} \quad\text{(media generale, o grand mean)}. \]
69.1.2 Effetti principali e interazione
Se i fattori \(R\) e \(C\) non interagiscono nel determinare la variabile di risposta, allora l’effetto di uno di essi, a parità di categoria dell’altro, rimane costante. In termini pratici, la differenza fra le medie di cella \(\mu_{jk} - \mu_{j'k}\)—quando confrontiamo due categorie di \(R\), ad esempio \(R_j\) e \(R_{j'}\)—è la stessa per tutte le categorie di \(C\) (cioè per \(k = 1, 2, \dots, c\)).
Di conseguenza, la differenza fra le medie nelle righe è uguale alla differenza fra le corrispondenti medie marginali di riga:
\[ \mu_{jk} - \mu_{j'k} \;=\; \mu_{jk'} - \mu_{j'k'} \;=\; \mu_{j:} - \mu_{j':} \quad\text{per ogni } j, j' \text{ e } k, k'. \]
Un altro modo di vedere questa assenza di interazione è attraverso i profili di medie di cella, che in questo caso risultano “paralleli”. Se i profili sono paralleli, allora la differenza fra \(\mu_{j1}\) e \(\mu_{j2}\) (categorie \(C_1\) e \(C_2\)) resta costante fra le diverse righe \(j = 1, 2\), e coincide con la differenza fra le medie marginali di colonna, \(\mu_{:1} - \mu_{:2}\).
L’interazione è un concetto simmetrico: se il fattore \(R\) interagisce con il fattore \(C\), vale anche il contrario. Quando non si verifica interazione, l’effetto principale (o main effect) di ogni fattore corrisponde semplicemente alla differenza fra le medie marginali di popolazione relative a quel fattore.
La figura seguente presentata diversi scenari possibili.
69.2 Simulazione
Per fare un esempio concreto, simuliamo dei dati seguendo lo schema utilizzato con l’ANOVA ad una via. In questo caso, aggiungiamo un secondo fattore: la gravità dei sintomi.
# Imposta un seme per riproducibilità (opzionale)
set.seed(123)
# Definiamo il numero di osservazioni per ogni cella
n <- 30
# Definiamo le categorie dei fattori
gravita <- c("molto_gravi", "poco_gravi")
condizione <- c("controllo", "psicoterapia1", "psicoterapia2")
# Definiamo la deviazione standard (ipotizziamo la stessa per tutte le celle)
sd_value <- 5
# Definiamo le medie attese per ciascuna combinazione (gravita x condizione)
# Esempio:
# - Pazienti molto gravi: (controllo=30, psico1=25, psico2=20)
# - Pazienti poco gravi: (controllo=25, psico1=20, psico2=15)
mean_table <- matrix(
c(30, 25, 20, # molto_gravi
25, 20, 15), # poco_gravi
nrow = 2, # 2 righe per "molto_gravi" e "poco_gravi"
ncol = 3, # 3 colonne per "controllo", "psicoterapia1", "psicoterapia2"
byrow = TRUE
)
# Crea un data frame vuoto che riempiremo con i dati simulati
df <- data.frame()
for (i in seq_along(gravita)) {
for (j in seq_along(condizione)) {
# Estraiamo la media corrispondente alla combinazione (i, j)
current_mean <- mean_table[i, j]
# Generiamo 'n' osservazioni normali con media = current_mean e sd = sd_value
simulated_data <- rnorm(n, mean = current_mean, sd = sd_value)
# Creiamo un data frame temporaneo per questa combinazione
temp_df <- data.frame(
gravita = gravita[i],
condizione = condizione[j],
punteggio = simulated_data
)
# Append al data frame finale
df <- rbind(df, temp_df)
}
}Esaminiamo le medie:
# Esempio di sintesi statistica
aggregate(punteggio ~ gravita + condizione, data = df, mean)
#> gravita condizione punteggio
#> 1 molto_gravi controllo 29.76
#> 2 poco_gravi controllo 24.53
#> 3 molto_gravi psicoterapia1 25.89
#> 4 poco_gravi psicoterapia1 19.08
#> 5 molto_gravi psicoterapia2 20.12
#> 6 poco_gravi psicoterapia2 15.77Rappresentiamo graficamente la distribuzione dei dati nelle varie condizioni:
ggplot(df, aes(x = condizione, y = punteggio, fill = gravita)) +
geom_boxplot(position = position_dodge()) +
labs(title = "Simulazione dati: Effetto gravità e condizione",
x = "Condizione sperimentale",
y = "Punteggio")
Addattiamo ai dati un modello bayesiano:
mod <- brm(
punteggio ~ gravita * condizione,
data = df,
backend = "cmdstanr"
)Esaminiamo le stime del modello:
conditional_effects(mod, "condizione:gravita")
Un sommario delle stime a posteriori si ottiene nel modo seguente:
summary(mod)
#> Family: gaussian
#> Links: mu = identity; sigma = identity
#> Formula: punteggio ~ gravita * condizione
#> Data: df (Number of observations: 180)
#> Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
#> total post-warmup draws = 4000
#>
#> Regression Coefficients:
#> Estimate Est.Error l-95% CI
#> Intercept 29.75 0.86 28.03
#> gravitapoco_gravi -5.20 1.21 -7.61
#> condizionepsicoterapia1 -3.86 1.23 -6.24
#> condizionepsicoterapia2 -9.61 1.23 -12.05
#> gravitapoco_gravi:condizionepsicoterapia1 -1.61 1.72 -4.99
#> gravitapoco_gravi:condizionepsicoterapia2 0.83 1.72 -2.53
#> u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept 31.45 1.00 2133 2829
#> gravitapoco_gravi -2.87 1.00 1866 2352
#> condizionepsicoterapia1 -1.36 1.00 2240 2404
#> condizionepsicoterapia2 -7.16 1.00 2274 2624
#> gravitapoco_gravi:condizionepsicoterapia1 1.85 1.00 1924 2338
#> gravitapoco_gravi:condizionepsicoterapia2 4.23 1.00 2071 2416
#>
#> Further Distributional Parameters:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma 4.79 0.25 4.33 5.32 1.00 3481 2699
#>
#> Draws were sampled using sample(hmc). 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).L’esame dell’output precedente mostra come sia difficile rispondere alle domende di interesse mediante l’esame dei singoli coefficienti. Per valutare gli effetti principali e la presenza di interazione si usa invece un metodo basato sul confronto tra modelli.
69.3 Confronto tra Modelli
Il test degli effetti principali e dell’interazione viene eseguito, nell’approccio frequentista, calcolando R^2 per i vari modelli per poi fare una differenza tra gli R^2 pesati per i gradi di libertà. Questa differenza si distribuisce come F e quindi disponiamo di una distribuzione campionaria nota per questa statistica. Il test si fa nel solito modo, ovvero ci si chiede se la differenza in R^2 osservata è sufficientemente piccola da potere essere attribuita al caso, sotto l’ipotesi che i due modelli sono equivalenti, oppure no.
I test bayesiani si basano su una logica diversa, anche se ha alcuni punti in comune. Il confronto tra modelli si basa su una quantità chiamata LOO (Leave-One-Out cross-validation).
69.4 Cos’è LOO (Leave-One-Out cross-validation)?
LOO è un metodo per stimare quanto bene un modello predice nuovi dati. In particolare, misura quanto il modello è in grado di generalizzare oltre i dati usati per costruirlo.
Come funziona:
Si rimuove iterativamente una osservazione dal dataset, si stima il modello sui dati rimanenti e si valuta quanto bene il modello predice l’osservazione esclusa. Questo processo viene ripetuto per tutte le osservazioni.Nel contesto Bayesiano:
Poiché ricalcolare il modello per ogni possibile sottogruppo di dati sarebbe computazionalmente oneroso, si usa una stima approssimativa di LOO basata sul concetto di PSIS-LOO (Pareto Smoothed Importance Sampling). Questa stima è disponibile direttamente inbrmstramite il comandoloo().
69.5 Indicatori principali nei risultati di LOO
elpd(expected log predictive density):
È la somma logaritmica delle probabilità predittive del modello per ogni osservazione, calcolata dopo aver escluso quella stessa osservazione. Valori più alti indicano un modello che predice meglio.elpd_diff:
Differenza dielpdtra due modelli. Se \(\text{elpd}_{\text{model1}} > \text{elpd}_{\text{model2}}\), il modello 1 è preferito.se_diff(standard error of difference):
L’incertezza associata aelpd_diff. Se|elpd_diff|è molto piccolo rispetto ase_diff, la differenza tra i modelli non è considerata robusta.
Nel caso presente, confrontiamo due modelli:
-
mod1: Modello completo (con interazione).
- mod: Modello senza interazione (solo effetti principali).
mod1 <- brm(
punteggio ~ gravita + condizione,
data = df,
backend = "cmdstanr"
)69.5.1 Confronto LOO
# Confronto via LOO
loo_full <- loo(mod)
loo_noint <- loo(mod1)
loo_compare(loo_full, loo_noint)
#> elpd_diff se_diff
#> mod1 0.0 0.0
#> mod -0.7 1.5elpd_diff = -0.9: Il modello con interazione (mod) ha un elpd leggermente più basso rispetto al modello senza interazione (mod1). Questo significa che il modello senza interazione è leggermente preferibile in termini di capacità predittiva sui nuovi dati.se_diff = 1.5: L’incertezza associata alla differenza di elpd è molto più grande della differenza stessa (elpd_diff/se_diff). Questo implica che non c’è evidenza sufficiente per concludere che uno dei due modelli sia chiaramente migliore.
In conclusione,
Il modello senza interazione (
mod1) è leggermente preferibile in termini di LOO-IC, ma la differenza è trascurabile e non robusta, data l’elevata incertezza (se_diff= 1.5).-
Interpretazione pratica:
- Non ci sono prove sufficienti per preferire il modello con interazione (mod) rispetto al modello senza interazione (mod1).
- In assenza di altre considerazioni teoriche che giustifichino l’inclusione dell’interazione, il modello più semplice (mod1, senza interazione) è probabilmente una scelta migliore.
Si procede in un modo simile per i test degli effetti principali. In quel caso potremmo confrontrare un modello con gli effetti additivi dei due fattori con un modello con un unico fattore.
Informazioni sull’Ambiente di Sviluppo
sessionInfo()
#> R version 4.4.2 (2024-10-31)
#> Platform: aarch64-apple-darwin20
#> Running under: macOS Sequoia 15.3.1
#>
#> Matrix products: default
#> BLAS: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
#>
#> 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] rstan_2.32.6 StanHeaders_2.32.10 bridgesampling_1.1-2
#> [4] loo_2.8.0 emmeans_1.10.7 brms_2.22.0
#> [7] Rcpp_1.0.14 bayestestR_0.15.2 posterior_1.6.1
#> [10] cmdstanr_0.8.1 thematic_0.1.6 MetBrewer_0.2.0
#> [13] ggokabeito_0.1.0 see_0.10.0 gridExtra_2.3
#> [16] patchwork_1.3.0 bayesplot_1.11.1 psych_2.4.12
#> [19] scales_1.3.0 markdown_1.13 knitr_1.49
#> [22] lubridate_1.9.4 forcats_1.0.0 stringr_1.5.1
#> [25] dplyr_1.1.4 purrr_1.0.4 readr_2.1.5
#> [28] tidyr_1.3.1 tibble_3.2.1 ggplot2_3.5.1
#> [31] tidyverse_2.0.0 rio_1.2.3 here_1.0.1
#>
#> loaded via a namespace (and not attached):
#> [1] tidyselect_1.2.1 farver_2.1.2 fastmap_1.2.0
#> [4] TH.data_1.1-3 tensorA_0.36.2.1 pacman_0.5.1
#> [7] digest_0.6.37 timechange_0.3.0 estimability_1.5.1
#> [10] lifecycle_1.0.4 survival_3.8-3 processx_3.8.6
#> [13] magrittr_2.0.3 compiler_4.4.2 rlang_1.1.5
#> [16] tools_4.4.2 yaml_2.3.10 data.table_1.17.0
#> [19] labeling_0.4.3 htmlwidgets_1.6.4 curl_6.2.1
#> [22] pkgbuild_1.4.6 mnormt_2.1.1 plyr_1.8.9
#> [25] abind_1.4-8 multcomp_1.4-28 withr_3.0.2
#> [28] stats4_4.4.2 grid_4.4.2 inline_0.3.21
#> [31] xtable_1.8-4 colorspace_2.1-1 MASS_7.3-65
#> [34] insight_1.1.0 cli_3.6.4 mvtnorm_1.3-3
#> [37] rmarkdown_2.29 generics_0.1.3 RcppParallel_5.1.10
#> [40] rstudioapi_0.17.1 reshape2_1.4.4 tzdb_0.4.0
#> [43] splines_4.4.2 parallel_4.4.2 matrixStats_1.5.0
#> [46] vctrs_0.6.5 V8_6.0.1 Matrix_1.7-2
#> [49] sandwich_3.1-1 jsonlite_1.9.1 hms_1.1.3
#> [52] glue_1.8.0 codetools_0.2-20 ps_1.9.0
#> [55] distributional_0.5.0 stringi_1.8.4 gtable_0.3.6
#> [58] QuickJSR_1.6.0 munsell_0.5.1 pillar_1.10.1
#> [61] htmltools_0.5.8.1 Brobdingnag_1.2-9 R6_2.6.1
#> [64] rprojroot_2.0.4 evaluate_1.0.3 lattice_0.22-6
#> [67] backports_1.5.0 rstantools_2.4.0 coda_0.19-4.1
#> [70] nlme_3.1-167 checkmate_2.3.2 xfun_0.51
#> [73] zoo_1.8-13 pkgconfig_2.0.3Bibliografia
Johnson, A. A., Ott, M., & Dogucu, M. (2022). Bayes Rules! An Introduction to Bayesian Modeling with R. CRC Press.