69  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).

69.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)

69.2 Introduzione

Nel Capitolo 67 abbiamo visto come rappresentare fattori a due livelli mediante variabili dummy, per includere effetti categoriali all’interno di un modello lineare. In questa sezione estendiamo tale approccio al caso in cui un singolo fattore abbia più di due categorie.

Supponiamo, ad esempio, di voler analizzare un fattore con tre gruppi. In questo caso possiamo utilizzare due variabili dummy per rappresentare le prime due categorie, e assumere implicitamente la terza come categoria di riferimento. Il modello lineare diventa:

\[ Y_i = \alpha + \gamma_1 D_{i1} + \gamma_2 D_{i2} + \varepsilon_i \tag{69.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 seguente tabella mostra come codificare le dummy per ciascun gruppo:

\[ \begin{array}{c|cc} \text{Gruppo} & D_{1} & D_{2} \\ \hline 1 & 1 & 0 \\ 2 & 0 & 1 \\ 3 & 0 & 0 \end{array} \tag{69.2}\]

69.2.1 Significato dei parametri

L’obiettivo è stimare le medie di popolazione per ciascun gruppo, indicate con \(\mu_j\) per il gruppo \(j\). Poiché l’errore \(\varepsilon\) ha media zero, possiamo calcolare l’aspettativa del modello per ciascun gruppo:

\[ \begin{aligned} \text{Gruppo 1: } \mu_1 &= \mathbb{E}(Y \mid D_1 = 1, D_2 = 0) = \alpha + \gamma_1 \\ \text{Gruppo 2: } \mu_2 &= \mathbb{E}(Y \mid D_1 = 0, D_2 = 1) = \alpha + \gamma_2 \\ \text{Gruppo 3: } \mu_3 &= \mathbb{E}(Y \mid D_1 = 0, D_2 = 0) = \alpha \end{aligned} \tag{69.3}\]

Queste espressioni mostrano come i coefficienti del modello siano relativi alla categoria di riferimento, cioè il Gruppo 3. Possiamo dunque esprimere i parametri del modello in funzione delle medie di gruppo:

\[ \alpha = \mu_3, \quad \gamma_1 = \mu_1 - \mu_3, \quad \gamma_2 = \mu_2 - \mu_3 . \tag{69.4}\]

In sintesi:

  • \(\alpha\) rappresenta la media del gruppo di riferimento (Gruppo 3),
  • \(\gamma_1\) è la differenza tra la media del Gruppo 1 e quella del Gruppo 3,
  • \(\gamma_2\) è la differenza tra la media del Gruppo 2 e quella del Gruppo 3.

Questa struttura rende il modello facilmente interpretabile: i coefficienti delle dummy indicano quanto ciascun gruppo si discosta dalla categoria di riferimento, mentre l’intercetta fornisce direttamente la media di quest’ultima. Questo schema può essere esteso a fattori con un numero arbitrario di categorie.

69.3 Simulazione

Per illustrare il funzionamento di un modello con un fattore a tre livelli, simuliamo un dataset con tre condizioni sperimentali: controllo, psicoterapia1 e psicoterapia2. Supponiamo che ogni gruppo abbia una diversa media, ma la stessa deviazione standard:

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

69.3.1 Esplorazione dei dati

Visualizziamo le distribuzioni con un violin plot arricchito da un boxplot:

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

69.4 Codifica 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 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.

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

69.5 Contrasti personalizzati

L’uso di contrasti personalizzati ci permette di testare ipotesi più specifiche. Ad esempio:

  • contrasto 1: controllare se la media del gruppo controllo differisce dalla media combinata delle due psicoterapie;
  • contrasto 2: confrontare direttamente psicoterapia1 con psicoterapia2.

69.5.1 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
mod_custom <- lm(punteggio ~ condizione, data = df)
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

69.6 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

69.7 Uso del pacchetto emmeans

Il pacchetto emmeans consente di ottenere gli stessi risultati in modo più diretto e modulare.

69.7.1 Stima con brms

Usiamo ora emmeans con brms:

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

69.7.2 Calcolo delle medie marginali

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

69.7.3 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

69.7.4 Contrasti personalizzati con emmeans

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)

69.8 Riflessioni Conclusive

Un’ANOVA a una via può essere vista come un caso speciale di regressione lineare, applicato a un fattore con più di due livelli. L’elemento più informativo dell’analisi non è tanto la verifica globale dell’effetto del fattore, ma l’esplorazione mirata di contrasti specifici tra le medie.

La possibilità di definire contrasti personalizzati, come quelli che confrontano un singolo gruppo con la media degli altri, oppure due gruppi tra loro, permette di rispondere a domande teoriche mirate. Questo approccio si integra perfettamente con il framework bayesiano, in cui i contrasti possono essere interpretati in termini di probabilità a posteriori, anziché solo di significatività.

Il pacchetto emmeans rende questo processo ancora più accessibile, offrendo strumenti per:

  • calcolare medie marginali stimate;
  • specificare contrasti personalizzati;
  • ottenere intervalli di confidenza (o credibilità) e p-value corretti.

In definitiva, l’uso combinato di regressione, contrasti e strumenti come emmeans o brms fornisce un approccio potente, flessibile e interpretabile per analizzare l’effetto di fattori categoriali in modelli lineari.

Informazioni sull’Ambiente di Sviluppo

sessionInfo()
#> R version 4.5.0 (2025-04-11)
#> Platform: aarch64-apple-darwin20
#> Running under: macOS Sequoia 15.4.1
#> 
#> 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.15.3
#>  [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.0    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.2.0        cli_3.6.5           
#> [37] mvtnorm_1.3-3        rmarkdown_2.29       generics_0.1.3      
#> [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

Bibliografia

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