66  ANOVA ad una via

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 Geocentric models di Statistical rethinking (McElreath, 2020).

66.1 Preparazione del Notebook

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

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

66.2 Introduzione

Nel Capitolo 63 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.

66.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{66.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{66.2}\]

66.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.

66.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

66.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:

df |> 
  group_by(condizione) |> 
  summarize(
    media = mean(punteggio),
    sd = sd(punteggio)
  )
#> # A tibble: 3 × 3
#>   condizione    media    sd
#>   <chr>         <dbl> <dbl>
#> 1 controllo      29.8  4.91
#> 2 psicoterapia1  25.9  4.18
#> 3 psicoterapia2  20.1  4.35

66.5 Modello lineare con variabili dummy

Convertiamo condizione in fattore e definiamo controllo come categoria di riferimento:

df$condizione <- factor(df$condizione)
df$condizione <- relevel(df$condizione, ref = "controllo")
contrasts(df$condizione)
#>               psicoterapia1 psicoterapia2
#> controllo                 0             0
#> psicoterapia1             1             0
#> psicoterapia2             0             1

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.

66.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

66.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

66.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:

my_list <- list(
  "Ctrl_vs_PsicoMean" = c(
    "controllo" = 1, "psicoterapia1" = -0.5, "psicoterapia2" = -0.5
  ),
  "P1_vs_P2" = c(
    "controllo" = 0, "psicoterapia1" = 1, "psicoterapia2" = -1
  )
)
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)

66.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.7    MetBrewer_0.2.0  
#>  [9] ggokabeito_0.1.0  see_0.11.0        gridExtra_2.3     patchwork_1.3.0  
#> [13] bayesplot_1.13.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.3.0      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        TH.data_1.1-3        tensorA_0.36.2.1    
#>  [7] pacman_0.5.1         digest_0.6.37        timechange_0.3.0    
#> [10] estimability_1.5.1   lifecycle_1.0.4      StanHeaders_2.32.10 
#> [13] survival_3.8-3       processx_3.8.6       magrittr_2.0.3      
#> [16] compiler_4.5.0       rlang_1.1.6          tools_4.5.0         
#> [19] utf8_1.2.6           yaml_2.3.10          data.table_1.17.6   
#> [22] labeling_0.4.3       bridgesampling_1.1-2 htmlwidgets_1.6.4   
#> [25] curl_6.3.0           pkgbuild_1.4.8       mnormt_2.1.1        
#> [28] plyr_1.8.9           RColorBrewer_1.1-3   abind_1.4-8         
#> [31] multcomp_1.4-28      withr_3.0.2          stats4_4.5.0        
#> [34] grid_4.5.0           colorspace_2.1-1     inline_0.3.21       
#> [37] xtable_1.8-4         MASS_7.3-65          insight_1.3.0       
#> [40] cli_3.6.5            mvtnorm_1.3-3        rmarkdown_2.29      
#> [43] generics_0.1.4       RcppParallel_5.1.10  rstudioapi_0.17.1   
#> [46] reshape2_1.4.4       tzdb_0.5.0           rstan_2.32.7        
#> [49] splines_4.5.0        parallel_4.5.0       matrixStats_1.5.0   
#> [52] vctrs_0.6.5          V8_6.0.4             Matrix_1.7-3        
#> [55] sandwich_3.1-1       jsonlite_2.0.0       hms_1.1.3           
#> [58] glue_1.8.0           codetools_0.2-20     ps_1.9.1            
#> [61] distributional_0.5.0 stringi_1.8.7        gtable_0.3.6        
#> [64] QuickJSR_1.8.0       pillar_1.10.2        htmltools_0.5.8.1   
#> [67] Brobdingnag_1.2-9    R6_2.6.1             rprojroot_2.0.4     
#> [70] evaluate_1.0.4       lattice_0.22-7       backports_1.5.0     
#> [73] rstantools_2.4.0     coda_0.19-4.1        nlme_3.1-168        
#> [76] checkmate_2.3.2      xfun_0.52            zoo_1.8-14          
#> [79] pkgconfig_2.0.3

Bibliografia

McElreath, R. (2020). Statistical rethinking: A Bayesian course with examples in R and Stan (2nd Edition). CRC Press.