Come Gestire le Covariate in Statistica: Dalla Regressione ai Modelli a Effetti Misti in Python

Quante volte ci è capitato di guardare il coefficiente di un modello e chiederci:

“Ma questo è l’effetto reale, o c’è sotto qualcos’altro?”

Lavorare con i dati significa fare costantemente i conti con il rumore e con variabili che si intrecciano tra loro.

L’istinto, spesso, ci suggerisce di buttare nel calderone del modello tutte le variabili a disposizione per “controllare” il più possibile.

Eppure, aggiungere una variabile non è mai un’operazione neutra.

In statistica, quando inseriamo una variabile extra per pulire il nostro segnale o per ridurre la varianza, stiamo maneggiando una covariata. Non si tratta di un semplice trucco matematico per alzare l'[math]R^2[/math], ma di una lente che cambia letteralmente la domanda a cui il nostro modello sta rispondendo.

In questo articolo imparerai a:

• Distinguere tra covariata e confondente
• Stimare l’effetto di una covariata continua in una regressione lineare e riconoscere l’omitted variable bias
• Condurre un’ANCOVA in Python, verificare l’assunzione delle pendenze parallele e confrontare codifiche ordinali e categoriche
• Gestire covariate time-varying in studi longitudinali con Linear Mixed Models
• Separare effetti between-person e within-person e interpretare l’ICC
• Scegliere il modello adeguato in base alla struttura dei dati

---

1. Cos’è una covariata?

In statistica, una covariata è una variabile che viene inclusa in un modello insieme alla variabile esplicativa principale e alla variabile risposta.

Immaginiamo di voler studiare l’effetto di un trattamento [math]T[/math] su un outcome [math]Y[/math]. Un modello lineare semplice è:

[math]\displaystyle Y_i = \beta_0 + \beta_1 T_i + \varepsilon_i[/math]

Se esiste una seconda variabile [math]X[/math] — per esempio l’età del soggetto — possiamo scrivere:

[math]\displaystyle Y_i = \beta_0 + \beta_1 T_i + \beta_2 X_i + \varepsilon_i[/math]

Qui [math]X[/math] è una covariata. Il coefficiente [math]\beta_1[/math] rappresenta l’associazione tra trattamento e outcome a parità di covariata inclusa nel modello.

Covariata e confondente non sono sinonimi

Questa distinzione è fondamentale.

Una covariata è una variabile inserita nel modello.

Un confondente è una variabile che, in un problema causale specifico, è associata sia all’esposizione/trattamento sia all’outcome e può distorcere la relazione osservata.

Una variabile confondente può essere utilizzata come covariata, ma non tutte le covariate sono confondenti. Per esempio, in uno studio randomizzato l’età può essere una covariata prognostica anche quando la randomizzazione ha già eliminato il confondimento sistematico.

---

2. Esercizio 1 — Regressione lineare con covariata continua

Il problema

Un’azienda farmaceutica testa un farmaco per ridurre la pressione sanguigna. Vengono reclutati 200 pazienti e raccolte:

• pressione sistolica dopo 3 mesi
• tipo di trattamento (farmaco o placebo)
• età del paziente

La domanda: perché l’età dovrebbe essere considerata una covariata? E cosa succede se viene ignorata?

Simuliamo i dati con omitted variable bias

Per capire il problema costruiamo artificialmente un dataset in cui il trattamento è correlato all’età: i pazienti più anziani hanno maggiore probabilità di ricevere il farmaco. Questo scenario è realistico in studi osservazionali o in campioni randomizzati di piccole dimensioni dove si verificano sbilanciamenti casuali.

import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression

np.random.seed(42)
n = 200

eta = np.random.uniform(50, 80, n)

# Pazienti più anziani hanno maggiore probabilità di ricevere il farmaco
prob_tratt = 0.2 + 0.6 * (eta - 50) / 30
trattamento = (np.random.uniform(0, 1, n) < prob_tratt).astype(int)

# Modello generativo vero
pressione = (
    120
    + 10 * trattamento      # vero effetto del farmaco
    + 0.5 * eta             # effetto dell'età
    + np.random.normal(0, 5, n)
)

df = pd.DataFrame({
    "pressione": pressione,
    "trattamento": trattamento,
    "eta": eta
})

Il vantaggio della simulazione è che conosciamo la verità generativa: effetto del trattamento = 10, effetto dell’età = 0.5.

Modello 1: ignoriamo la covariata

X1 = df[["trattamento"]]
y = df["pressione"]

mod1 = LinearRegression().fit(X1, y)
print(f"Coefficiente trattamento: {mod1.coef_[0]:.2f}")
print(f"R²: {mod1.score(X1, y):.3f}")

Output:

Coefficiente trattamento: 13.41
R²: 0.523

Il valore stimato (13.41) è superiore al vero effetto di 10. Questo è l’omitted variable bias: poiché il trattamento è correlato positivamente con l’età, e l’età aumenta la pressione, il coefficiente del trattamento assorbe anche parte dell’effetto dell’età.

Modello 2: inseriamo la covariata

X2 = df[["trattamento", "eta"]]
mod2 = LinearRegression().fit(X2, y)

print(f"Coefficiente trattamento: {mod2.coef_[0]:.2f}")
print(f"Coefficiente età: {mod2.coef_[1]:.2f}")
print(f"R²: {mod2.score(X2, y):.3f}")

Output:

Coefficiente trattamento: 10.20
Coefficiente età: 0.50
R²: 0.718

Il modello con covariata recupera molto bene la struttura generativa. Il coefficiente del trattamento si avvicina al vero valore 10, e l'[math]R^2[/math] aumenta da 0.523 a 0.718.

Controllo delle assunzioni

Prima di fidarci dei coefficienti, verifichiamo visivamente che la relazione tra età e pressione sia approssimativamente lineare in entrambi i gruppi:

import seaborn as sns
import matplotlib.pyplot as plt

sns.lmplot(x='eta', y='pressione', hue='trattamento', data=df,
           palette={0: '#E74C3C', 1: '#3498DB'}, height=5, aspect=1.2)
plt.title('Regressione di pressione su età, per gruppo')
plt.show()

Se le rette nei due gruppi fossero chiaramente non parallele, dovremmo considerare un’interazione [math]\text{trattamento} \times \text{età}[/math].

Cosa abbiamo imparato?

La covariata ha due funzioni distinte:

• Spiegare variabilità — L’età contribuisce a spiegare la pressione, aumentando l'[math]R^2[/math].
• Isolare l’associazione di interesse — Il coefficiente del trattamento viene stimato tenendo conto dell’età, riducendo il bias.

Randomizzazione: le covariate sono inutili?

No. La randomizzazione protegge dal confondimento sistematico, ma in campioni finiti possono verificarsi sbilanciamenti casuali. Una covariata fortemente prognostica aumenta la precisione della stima dell’effetto del trattamento anche in studi randomizzati.

---

3. Esercizio 2 — ANCOVA con covariata ordinale

Il problema

Un ricercatore confronta tre metodi di insegnamento (A, B, C) per studiarne l’effetto sul punteggio d’esame. Raccoglie anche il livello di istruzione dei genitori: Basso, Medio, Alto.

La domanda: il metodo ha effetto sul punteggio dopo aver tenuto conto dell’istruzione genitoriale?

ANCOVA: un ponte tra ANOVA e regressione

L’ANCOVA (Analysis of Covariance) confronta i gruppi tenendo contemporaneamente conto di una o più covariate:

[math]\displaystyle \text{Punteggio}_i = \beta_0 + \beta_1 \text{Metodo}_i + \beta_2 \text{Istruzione}_i + \varepsilon_i[/math]

Simuliamo un dataset realistico

Generiamo 60 studenti (20 per gruppo):

import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.formula.api import ols

np.random.seed(456)

righe = []
for metodo in ['A', 'B', 'C']:
    for _ in range(20):
        istruzione = np.random.choice([1, 2, 3], p=[0.3, 0.4, 0.3])
        effetto_metodo = {'A': 0, 'B': 3, 'C': 1.5}[metodo]
        punteggio = 50 + effetto_metodo + 4 * istruzione + np.random.normal(0, 3)
        righe.append({'metodo': metodo, 'istruzione': istruzione, 'punteggio': punteggio})

df2 = pd.DataFrame(righe)

Codifica ordinale vs categorica

L’istruzione è ordinale. Possiamo codificarla come numerica (1, 2, 3), ma questa scelta assume che il passaggio da Basso a Medio sia equivalente a quello da Medio ad Alto. Verifichiamo se questa assunzione è sostenuta dai dati:

# Modello 1: istruzione come ordinale (lineare)
model_ord = ols("punteggio ~ C(metodo) + istruzione", data=df2).fit()

# Modello 2: istruzione come categorica
model_cat = ols("punteggio ~ C(metodo) + C(istruzione)", data=df2).fit()

print(f"AIC ordinale: {model_ord.aic:.1f}")
print(f"AIC categorico: {model_cat.aic:.1f}")

Output:

AIC ordinale: 298.0
AIC categorico: 299.9

L’AIC più basso del modello ordinale suggerisce che l’assunzione di linearità è ragionevole, ma la differenza è modesta. In un articolo scientifico converrebbe riportare entrambi i modelli o giustificare la scelta con ragionamento sostanziale.

Verifica dell’assunzione delle pendenze parallele

L’ANCOVA assume che l’effetto della covariata sia lo stesso in tutti i gruppi. Testiamo l’interazione:

model_int = ols("punteggio ~ C(metodo) * istruzione", data=df2).fit()
print(sm.stats.anova_lm(model_int, typ=2))

Output:

                          sum_sq    df         F    PR(>F)
C(metodo)             202.368     2.0  12.602  3.22e-05
istruzione            508.098     1.0  63.282  1.17e-10
C(metodo):istruzione    7.485     2.0   0.466     0.630
Residual              433.572    54.0     NaN       NaN

L’interazione non è significativa ([math]p = 0.630[/math]): le pendenze sono parallele e l’ANCOVA è appropriata.

Stima e confronti post-hoc

# Tabella ANOVA
anova_table = sm.stats.anova_lm(model_ord, typ=2)
print(anova_table)

Output:

               sum_sq    df          F        PR(>F)
C(metodo)    202.368   2.0  12.602119  3.224215e-05
istruzione   508.098   1.0  63.282020  1.166003e-10
Residual     433.572  54.0        NaN           NaN

Il metodo di insegnamento rimane significativo dopo aver aggiustato per l’istruzione.

Estimated Marginal Means (EMM)

Le medie grezze non sono il confronto più informativo. Calcoliamo le medie marginali stimate fissando la covariata al livello medio ([math]\text{istruzione} = 2[/math]):

pred = model_ord.predict(pd.DataFrame({
    "metodo": ["A", "B", "C"],
    "istruzione": [2, 2, 2]
}))

for metodo, valore in zip(["A", "B", "C"], pred):
    print(f"Metodo {metodo}: {valore:.2f}")

Output:

Metodo A: 57.63
Metodo B: 62.08
Metodo C: 59.19

Confronti post-hoc corretti

Non usare Tukey sulle medie grezze: i confronti devono essere coerenti con il modello.

from statsmodels.stats.multitest import multipletests
import numpy as np

contrasts = np.array([
    [0, 1, 0, 0],   # B - A
    [0, 0, 1, 0],   # C - A
    [0, 1, -1, 0]   # B - C
])

test = model_ord.t_test(contrasts)
p_values = np.asarray(test.pvalue).ravel()
reject, p_adjusted, _, _ = multipletests(p_values, method="holm")

print(test)
print("p-value corretti:", p_adjusted)

Output:

p-value corretti: [1.82e-05, 8.40e-02, 3.93e-03]

Dopo l’aggiustamento per l’istruzione e la correzione di Holm: B > A ([math]p_{\text{adj}} < 0.001[/math]), B > C ([math]p_{\text{adj}} = 0.004[/math]), mentre C vs A non è significativo ([math]p_{\text{adj}} = 0.084[/math]).

---

4. Esercizio 3 — Covariata time-varying in uno studio longitudinale

Il problema

Studio clinico di 12 mesi su 200 pazienti con diabete di tipo 2. Ogni paziente è osservato a 0, 3, 6, 9, 12 mesi. Misuriamo:

• HbA1c (outcome)
• BMI (covariata che cambia nel tempo)
• tempo
• trattamento (fisso nel tempo)

Perché la regressione lineare standard non basta

Le 1000 righe del dataset non sono indipendenti: 5 osservazioni appartengono allo stesso paziente. Questa dipendenza gerarchica viola un’assunzione fondamentale dell’OLS.

Simulazione del dataset longitudinale

import numpy as np
import pandas as pd

np.random.seed(123)

n_pazienti = 200
tempi = np.array([0, 3, 6, 9, 12])
righe = []

for paziente in range(n_pazienti):
    trattamento = np.random.binomial(1, 0.5)
    bmi_base = np.random.normal(29, 3)
    effetto_individuale = np.random.normal(0, 0.6)  # intercetta casuale
    slope_individuale = np.random.normal(0, 0.03)   # pendenza casuale

    for tempo in tempi:
        bmi = bmi_base - 0.10 * tempo + np.random.normal(0, 0.35)
        hba1c = (
            8.0
            + effetto_individuale
            + slope_individuale * tempo
            - 0.08 * tempo
            - 0.45 * trattamento
            + 0.12 * (bmi - 29)
            + np.random.normal(0, 0.25)
        )
        righe.append({
            "paziente": paziente, "tempo": tempo,
            "trattamento": trattamento, "bmi": bmi, "hba1c": hba1c
        })

df_long = pd.DataFrame(righe)

Linear Mixed Model

from statsmodels.regression.mixed_linear_model import MixedLM

model_full = MixedLM.from_formula(
    "hba1c ~ trattamento + tempo + bmi",
    groups="paziente",
    re_formula="~tempo",   # intercetta + pendenza casuali
    data=df_long
)
result_full = model_full.fit()
print(result_full.summary())

Output:

                     Coef.  Std.Err.    z    P>|z| [0.025 0.975]
----------------------------------------------------------------
Intercept             4.322    0.384  11.255 0.000  3.569  5.074
trattamento          -0.473    0.086  -5.513 0.000 -0.642 -0.305
tempo                -0.076    0.003 -25.111 0.000 -0.082 -0.070
bmi                   0.128    0.013   9.953 0.000  0.103  0.153
paziente Var          0.323    0.169
paziente x tempo Cov  0.001    0.007
tempo Var             0.001    0.001

Attenzione al warning di convergenza

statsmodels emette un ConvergenceWarning: The MLE may be on the boundary of the parameter space. Questo accade perché la varianza della pendenza temporale (tempo Var ≈ 0.001) è praticamente zero. In altre parole, le traiettorie individuali non differiscono molto nella loro velocità di cambiamento.

Verifichiamo se un modello più semplice (solo intercetta casuale) è sufficiente:

model_int_only = MixedLM.from_formula(
    "hba1c ~ trattamento + tempo + bmi",
    groups="paziente",
    re_formula="~1",   # solo intercetta casuale
    data=df_long
)
result_int_only = model_int_only.fit()
print(f"Log-Likelihood full: {result_full.llf:.1f}")
print(f"Log-Likelihood int-only: {result_int_only.llf:.1f}")

Output:

Log-Likelihood full: -448.4
Log-Likelihood int-only: -482.7

Il modello full ha una log-likelihood nettamente superiore, quindi la struttura a intercetta e pendenza casuali è giustificata nonostante la varianza sul boundary. In un’analisi reale, si potrebbe anche considerare un modello con solo intercetta casuale se il guadagno in log-likelihood non compensasse la complessità aggiuntiva.

ICC: quanto conta la struttura gerarchica?

L’Intraclass Correlation Coefficient (ICC) misura la quota di variabilità dovuta alle differenze tra pazienti:

var_paziente = result_int_only.cov_re.iloc[0, 0]
var_residuo = result_int_only.scale
icc = var_paziente / (var_paziente + var_residuo)
print(f"ICC = {icc:.3f}")

Output:

ICC = 0.819

Un ICC dell’82% significa che l’82% della variabilità totale è spiegata dalle differenze tra pazienti. Ignorare questa struttura gerarchica sarebbe un errore grave: una regressione OLS sottostimerebbe gli errori standard e produrrebbe inferenze distorte.

Interpretazione degli effetti fissi

• Trattamento: [math]\hat{\beta} = -0.473[/math], [math]p < 0.001[/math]. A parità di tempo e BMI, il trattamento riduce l’HbA1c di circa 0.47 punti.
• Tempo: [math]\hat{\beta} = -0.076[/math] al mese, [math]p < 0.001[/math]. Tendenza lineare alla riduzione.
• BMI: [math]\hat{\beta} = 0.128[/math], [math]p < 0.001[/math]. Vicinissimo al vero valore generativo (0.12).

Decomposizione within/between

Una covariata time-varying contiene due informazioni distinte:

• Between-person: i pazienti con BMI medio più alto hanno HbA1c più alta?
• Within-person: quando un paziente aumenta il proprio BMI rispetto alla sua media, aumenta anche la sua HbA1c?

Decomponiamo il BMI:

df_long["bmi_mean"] = df_long.groupby("paziente")["bmi"].transform("mean")
df_long["bmi_within"] = df_long["bmi"] - df_long["bmi_mean"]

model_decomp = MixedLM.from_formula(
    "hba1c ~ trattamento + tempo + bmi_mean + bmi_within",
    groups="paziente",
    re_formula="~1",
    data=df_long
)
result_decomp = model_decomp.fit()
print(result_decomp.summary())

Output:

             Coef.  Std.Err.    z    P>|z| [0.025 0.975]
--------------------------------------------------------
bmi_mean      0.124    0.015   8.285 0.000  0.095  0.154
bmi_within    0.136    0.029   4.651 0.000  0.079  0.193

Entrambi gli effetti sono significativi e coerenti con il vero meccanismo generativo. Questa decomposizione è molto più informativa di un singolo coefficiente del BMI.

---

5. Errori frequenti nell’uso delle covariate

Errore	Descrizione	Soluzione
Confondere covariata e confondente	Pensare che siano sinonimi	Ricorda: covariata = variabile nel modello; confondente = concetto causale
Codifica ordinale senza riflettere	Basso=1, Medio=2, Alto=3 assume linearità	Verifica con modello categorico e confronta AIC
Correlazione ≠ controllo causale	Inserire una variabile non elimina automaticamente il bias	Ragiona sulla struttura causale (DAG)
Ignorare le interazioni	ANCOVA assume pendenze parallele	Testa sempre metodo × covariata
Post-hoc sulle medie grezze	Tukey sui dati grezzi ≠ confronti aggiustati	Effettua i contrasti sul modello stimato
Trattare dati longitudinali come indipendenti	OLS su misure ripetute viola le assunzioni	Usa Mixed Models e calcola l’ICC
Boundary warning ignorato	Varianza casuale ≈ 0 senza verifica	Confronta modelli nested con log-likelihood o test del rapporto di verosimiglianza

---

6. Covariate e causal inference

Un errore frequente è pensare: “Più covariate inserisco, meglio è.”

In un modello predittivo questo può essere ragionevole. In un’analisi causale, controllare per una variabile che è conseguenza del trattamento (un mediatore) può introdurre distorsioni.

Esempio:

Trattamento → BMI → Outcome

Qui il BMI è un mediatore. Inserirlo come covariata “blocca” parte del percorso causale e sottostima l’effetto totale del trattamento.

Regola pratica: prima di inserire una covariata, chiediti:

• È misurata prima del trattamento?
• È associata sia al trattamento sia all’outcome?
• È un confondente o un mediatore?
• È una collider (causa comune di due variabili già incluse)?

Un DAG (Directed Acyclic Graph) aiuta a visualizzare queste relazioni.

---

7. Quale modello scegliere?

Struttura dati	Modello	Covariata	Effetto casuale
1 osservazione per soggetto	OLS / Regressione lineare	Continua, ordinale, categorica	No
Confronto tra gruppi + covariata	ANCOVA / Regressione	Continua o ordinale	No
Covariata ordinale	Regressione con codifica appropriata	Ordinale (testa linearità)	No
Covariata categoriale	Variabili dummy / C(variabile)	Categorica	No
Misure ripetute, covariata fissa	Mixed Model	Continua	Intercetta
Misure ripetute, covariata time-varying	Mixed Model	Continua (decomposta)	Intercetta ± Pendenza
Interesse per effetti within/between	Mixed Model + decomposizione	Media + deviazione	Intercetta
Effetto diverso tra gruppi	Modello con interazione	Qualsiasi	A seconda dei dati

La domanda non è “quale funzione Python uso?” ma “qual è la struttura dei miei dati e quale domanda voglio rispondere?”

---

---

8. Checklist per l’analisi

Prima di inserire una covariata nel modello, verifica:

• [ ] Qual è la mia variabile risposta?
• [ ] Qual è l’effetto che voglio stimare?
• [ ] La covariata è misurata prima del trattamento?
• [ ] È continua, nominale o ordinale?
• [ ] La relazione con l’outcome è plausibilmente lineare?
• [ ] L’effetto della covariata può cambiare tra i gruppi? (testa l’interazione)
• [ ] La covariata cambia nel tempo?
• [ ] Le osservazioni sono indipendenti o annidate?
• [ ] Ci sono dati mancanti? (MCAR, MAR, MNAR?)
• [ ] Sto facendo previsione o inferenza causale?

---

9. Conclusione

La covariata non è un “aggiustamento automatico” da inserire per ripulire i dati.

È una scelta statistica che modifica:

• l’interpretazione dei coefficienti
• la precisione delle stime
• le assunzioni del modello
• in alcuni casi, la domanda causale stessa

Con Python possiamo partire da una semplice regressione:

[math]\displaystyle Y_i = \beta_0 + \beta_1 X_i + \varepsilon_i[/math]

e arrivare a modelli sofisticati:

[math]\displaystyle Y_{it} = \beta X_{it} + b_{0i} + b_{1i} \text{Time}_{it} + \varepsilon_{it}[/math]

Ma il punto fondamentale non è memorizzare LinearRegression, ols o MixedLM.

È capire perché un certo modello è appropriato per una certa struttura dei dati.

Prima di inserire una covariata, fermati e chiediti:

che ruolo ha questa variabile nel processo che sto osservando?

Da questa domanda dipende molto più della semplice formula del modello.