here::here("code", "_common.R") |>
source()
# Load packages
if (!requireNamespace("pacman")) install.packages("pacman")
pacman::p_load(cmdstanr, posterior, bayestestR, brms, emmeans)
68 ANOVA ad una via
68.1 Preparazione del Notebook
68.2 Introduzione
Nel Capitolo 65 ci siamo concentrati sul confronto tra due gruppi utilizzando una regressione lineare con variabili dummy. Questo approccio ci ha permesso di modellare in modo semplice l’effetto di un fattore binario e di stimare con incertezza l’ampiezza della differenza. Ora estendiamo quella logica al caso in cui il fattore abbia più di due livelli.
Questo passaggio ci introduce al cuore dell’ANOVA a una via, che non è altro che un modello lineare con un fattore categoriale a \(k\) livelli. In questo contesto, ci interessa capire quanta variabilità nei dati può essere attribuita alle differenze tra gruppi, e quanto invece rimane all’interno dei gruppi stessi. Come sempre in questo manuale, manterremo una lettura orientata all’incertezza e alla variabilità intra- e inter-individuale, trattando l’inferenza come uno strumento per quantificare la credibilità delle ipotesi.
68.3 Codifica del modello con variabili dummy
Supponiamo un esperimento con tre gruppi. Per rappresentare questo fattore all’interno di un modello lineare, usiamo due variabili dummy e consideriamo il terzo gruppo come riferimento implicito. Il modello assume la forma:
\[ Y_i = \alpha + \gamma_1 D_{i1} + \gamma_2 D_{i2} + \varepsilon_i \tag{68.1}\]
dove:
- \(\alpha\) è l’intercetta del modello,
- \(\gamma_1\) e \(\gamma_2\) sono i coefficienti associati alle variabili dummy,
- \(D_{i1}\) e \(D_{i2}\) indicano l’appartenenza dell’osservazione \(i\) ai gruppi 1 e 2, rispettivamente,
- \(\varepsilon_i\) è l’errore aleatorio.
La codifica delle dummy è la seguente:
\[ \begin{array}{c|cc} \text{Gruppo} & D_{1} & D_{2} \\ \hline 1 & 1 & 0 \\ 2 & 0 & 1 \\ 3 & 0 & 0 \end{array} \tag{68.2}\]
68.3.1 Interpretazione dei parametri
Con questa codifica, possiamo esprimere le medie di ciascun gruppo come:
\[ \begin{aligned} \mu_1 &= \alpha + \gamma_1 \\ \mu_2 &= \alpha + \gamma_2 \\ \mu_3 &= \alpha \end{aligned} \]
Da cui otteniamo:
\[ \alpha = \mu_3, \quad \gamma_1 = \mu_1 - \mu_3, \quad \gamma_2 = \mu_2 - \mu_3. \]
Quindi:
- \(\alpha\): media del gruppo 3 (riferimento),
- \(\gamma_1\): quanto il gruppo 1 si discosta da \(\mu_3\),
- \(\gamma_2\): quanto il gruppo 2 si discosta da \(\mu_3\).
In un’ottica bayesiana, questi coefficienti possono essere pensati come distribuzioni: esprimono quanto crediamo che ciascuna differenza sia plausibile, date le osservazioni. Passiamo ora a una simulazione.
68.4 Simulazione
Simuliamo un esperimento con tre condizioni: controllo
, psicoterapia1
e psicoterapia2
. Ogni gruppo ha una media diversa ma la stessa deviazione standard. Ci interessa modellare la variabilità tra le condizioni e interpretare le differenze in modo probabilistico.
set.seed(123)
n <- 30 # numero di osservazioni per gruppo
# Medie di ciascun gruppo
mean_control <- 30
mean_psico1 <- 25
mean_psico2 <- 20
# Deviazione standard comune
sd_value <- 5
# Generazione dei dati
controllo <- rnorm(n, mean_control, sd_value)
psicoterapia1 <- rnorm(n, mean_psico1, sd_value)
psicoterapia2 <- rnorm(n, mean_psico2, sd_value)
# Creazione del data frame
df <- data.frame(
condizione = rep(c("controllo", "psicoterapia1", "psicoterapia2"), each = n),
punteggio = c(controllo, psicoterapia1, psicoterapia2)
)
df |> head()
#> condizione punteggio
#> 1 controllo 27.20
#> 2 controllo 28.85
#> 3 controllo 37.79
#> 4 controllo 30.35
#> 5 controllo 30.65
#> 6 controllo 38.58
68.4.1 Esplorazione iniziale
Visualizziamo le distribuzioni dei punteggi:
ggplot(df, aes(x = condizione, y = punteggio, fill = condizione)) +
geom_violin(trim = FALSE) +
geom_boxplot(width = 0.2, outlier.shape = NA) +
labs(
title = "Distribuzione dei punteggi di depressione per gruppo",
x = "Condizione sperimentale",
y = "Punteggio di depressione"
) +
theme(legend.position = "none")
Calcoliamo media e deviazione standard per ogni gruppo:
68.5 Modello lineare con variabili dummy
Convertiamo condizione
in fattore e definiamo controllo
come categoria di riferimento:
Il modello di regressione con le variabili dummy sarà:
\[ Y_i = \beta_0 + \beta_1 \cdot \text{psicoterapia1}_i + \beta_2 \cdot \text{psicoterapia2}_i + \varepsilon_i, \]
dove:
- \(\beta_0\) è la media del gruppo di controllo;
- \(\beta_1\) e \(\beta_2\) sono le differenze tra le rispettive psicoterapie e il gruppo di controllo.
68.5.1 Stima del modello
Eseguiamo una prima analisi usando il metodo di massima verosimiglianza:
fm1 <- lm(punteggio ~ condizione, data = df)
summary(fm1)
#>
#> Call:
#> lm(formula = punteggio ~ condizione, data = df)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -11.668 -2.620 -0.183 2.681 10.128
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 29.764 0.819 36.33 < 2e-16
#> condizionepsicoterapia1 -3.873 1.159 -3.34 0.0012
#> condizionepsicoterapia2 -9.642 1.159 -8.32 1.1e-12
#>
#> Residual standard error: 4.49 on 87 degrees of freedom
#> Multiple R-squared: 0.446, Adjusted R-squared: 0.434
#> F-statistic: 35.1 on 2 and 87 DF, p-value: 6.75e-12
Verifica delle medie e differenze tra i gruppi:
out <- tapply(df$punteggio, df$condizione, mean)
out[2] - out[1] # psicoterapia1 - controllo
#> psicoterapia1
#> -3.873
out[3] - out[1] # psicoterapia2 - controllo
#> psicoterapia2
#> -9.642
68.6 Contrasti personalizzati
I contrasti ci permettono di andare oltre il test globale e formulare ipotesi teoriche mirate. Ad esempio:
- la media del gruppo controllo è diversa dalla media delle due psicoterapie?
- le due psicoterapie differiscono tra loro?
A questo fine, specifichiamo la seguente matrice dei contrasti:
my_contrasts <- matrix(c(
0.6667, 0, # controllo
-0.3333, 0.5, # psicoterapia1
-0.3333, -0.5 # psicoterapia2
), ncol = 2, byrow = TRUE)
colnames(my_contrasts) <- c("Ctrl_vs_PsicoMean", "P1_vs_P2")
rownames(my_contrasts) <- c("controllo", "psicoterapia1", "psicoterapia2")
contrasts(df$condizione) <- my_contrasts
Adattiamo il modello:
mod_custom <- lm(punteggio ~ condizione, data = df)
Esaminiamo i coefficienti:
summary(mod_custom)
#>
#> Call:
#> lm(formula = punteggio ~ condizione, data = df)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -11.668 -2.620 -0.183 2.681 10.128
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 25.259 0.473 53.40 < 2e-16
#> condizioneCtrl_vs_PsicoMean 6.758 1.003 6.73 1.7e-09
#> condizioneP1_vs_P2 5.770 1.159 4.98 3.2e-06
#>
#> Residual standard error: 4.49 on 87 degrees of freedom
#> Multiple R-squared: 0.446, Adjusted R-squared: 0.434
#> F-statistic: 35.1 on 2 and 87 DF, p-value: 6.75e-12
Interpretazione dei coefficienti:
- Intercetta: non rappresenta più una singola media, ma una combinazione lineare dei gruppi.
-
Ctrl_vs_PsicoMean: confronta la media di
controllo
con la media combinata delle due psicoterapie. - P1_vs_P2: differenza tra le due psicoterapie.
Verifica manuale:
# Controllo - media delle psicoterapie
out[1] - (out[2] + out[3]) / 2
#> controllo
#> 6.758
# Psicoterapia1 - Psicoterapia2
out[2] - out[3]
#> psicoterapia1
#> 5.77
68.7 Estensione bayesiana con brms
e emmeans
Usiamo ora il modello bayesiano:
mod <- brm(punteggio ~ condizione, data = df, backend = "cmdstanr")
summary(mod)
#> Family: gaussian
#> Links: mu = identity; sigma = identity
#> Formula: punteggio ~ condizione
#> Data: df (Number of observations: 90)
#> 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
#> Intercept 25.26 0.48 24.33 26.15 1.00
#> condizioneCtrl_vs_PsicoMean 6.78 1.04 4.73 8.85 1.00
#> condizioneP1_vs_P2 5.76 1.16 3.49 8.08 1.00
#> Bulk_ESS Tail_ESS
#> Intercept 4321 2937
#> condizioneCtrl_vs_PsicoMean 4260 2964
#> condizioneP1_vs_P2 4598 2785
#>
#> Further Distributional Parameters:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma 4.54 0.34 3.93 5.26 1.00 4287 3279
#>
#> 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).
Le medie marginali e i confronti possono essere ottenuti con il pacchetto emmeans
:
em <- emmeans(mod, specs = "condizione")
em
#> condizione emmean lower.HPD upper.HPD
#> controllo 29.8 28.1 31.4
#> psicoterapia1 25.9 24.3 27.5
#> psicoterapia2 20.1 18.4 21.7
#>
#> Point estimate displayed: median
#> HPD interval probability: 0.95
Confronti tra gruppi:
pairs(em) # confronti a coppie
#> contrast estimate lower.HPD upper.HPD
#> controllo - psicoterapia1 3.90 1.70 6.21
#> controllo - psicoterapia2 9.65 7.31 12.03
#> psicoterapia1 - psicoterapia2 5.75 3.57 8.14
#>
#> Point estimate displayed: median
#> HPD interval probability: 0.95
Contrasti personalizzati:
contrast(em, method = my_list)
#> contrast estimate lower.HPD upper.HPD
#> Ctrl_vs_PsicoMean 6.77 4.77 8.88
#> P1_vs_P2 5.75 3.57 8.14
#>
#> Point estimate displayed: median
#> HPD interval probability: 0.95
# Visualizzazione
plot(em)
68.8 Riflessioni Conclusive
L’ANOVA a una via è un esempio fondamentale di come un modello lineare possa rappresentare la variabilità tra gruppi. Tuttavia, il suo valore non sta nel test globale, ma nella possibilità di analizzare differenze mirate tra medie.
Attraverso contrasti personalizzati, possiamo porre domande teoriche precise e ottenere risposte in termini di effetti stimati con incertezza. Questo approccio si integra naturalmente con la prospettiva bayesiana, che ci permette di esprimere la probabilità che una certa differenza superi una soglia di interesse pratico.
Il pacchetto emmeans
(insieme a brms
) consente di navigare questa complessità in modo modulare e trasparente, producendo stime interpretabili e inference compatibili con i nostri modelli teorici. L’obiettivo non è semplicemente sapere se c’è una differenza, ma capire quanto, con quale incertezza e tra chi.
Informazioni sull’Ambiente di Sviluppo
sessionInfo()
#> R version 4.5.0 (2025-04-11)
#> Platform: aarch64-apple-darwin20
#> Running under: macOS Sequoia 15.5
#>
#> Matrix products: default
#> BLAS: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
#>
#> locale:
#> [1] C/UTF-8/C/C/C/C
#>
#> time zone: Europe/Rome
#> tzcode source: internal
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] emmeans_1.11.1 brms_2.22.0 Rcpp_1.0.14 bayestestR_0.16.0
#> [5] posterior_1.6.1 cmdstanr_0.9.0 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 estimability_1.5.1
#> [10] lifecycle_1.0.4 StanHeaders_2.32.10 processx_3.8.6
#> [13] magrittr_2.0.3 compiler_4.5.0 rlang_1.1.6
#> [16] tools_4.5.0 utf8_1.2.5 yaml_2.3.10
#> [19] data.table_1.17.2 labeling_0.4.3 bridgesampling_1.1-2
#> [22] htmlwidgets_1.6.4 curl_6.2.2 pkgbuild_1.4.7
#> [25] mnormt_2.1.1 plyr_1.8.9 RColorBrewer_1.1-3
#> [28] abind_1.4-8 withr_3.0.2 grid_4.5.0
#> [31] stats4_4.5.0 colorspace_2.1-1 xtable_1.8-4
#> [34] inline_0.3.21 insight_1.3.0 cli_3.6.5
#> [37] mvtnorm_1.3-3 rmarkdown_2.29 generics_0.1.4
#> [40] RcppParallel_5.1.10 rstudioapi_0.17.1 reshape2_1.4.4
#> [43] tzdb_0.5.0 rstan_2.32.7 parallel_4.5.0
#> [46] matrixStats_1.5.0 vctrs_0.6.5 V8_6.0.3
#> [49] Matrix_1.7-3 jsonlite_2.0.0 hms_1.1.3
#> [52] glue_1.8.0 codetools_0.2-20 ps_1.9.1
#> [55] distributional_0.5.0 stringi_1.8.7 gtable_0.3.6
#> [58] QuickJSR_1.7.0 pillar_1.10.2 htmltools_0.5.8.1
#> [61] Brobdingnag_1.2-9 R6_2.6.1 rprojroot_2.0.4
#> [64] evaluate_1.0.3 lattice_0.22-7 backports_1.5.0
#> [67] rstantools_2.4.0 coda_0.19-4.1 nlme_3.1-168
#> [70] checkmate_2.3.2 xfun_0.52 pkgconfig_2.0.3