10 - Regressione lineare semplice e multipla

Hai misurato il tempo di risposta di un servizio al variare del carico e vuoi prevedere quanto impiegherà con un carico che non hai provato. Tracci una retta sui punti, ci aggiungi una fascia di incertezza al 95%, e quella fascia sbaglia: se usi l'intervallo di confidenza che il tuo software chiama «95%», una misura futura ci cade dentro solo il 44,5% delle volte, come mostra la simulazione qui sotto. Esistono due intervalli diversi, uno per la media della risposta e uno per una singola misura futura, e confonderli è l'errore più frequente nella regressione.

Questo articolo, il decimo della serie Probabilità e Statistica, costruisce la retta dei minimi quadrati partendo da cinque punti, ti dice che cosa garantiscono le sue assunzioni e dove si rompono. Servono la teoria degli stimatori e il t-test degli articoli precedenti; dal chi-quadro ti porti l'idea di confrontare ciò che osservi con ciò che il modello prevede, che qui si applica ai residui.

La retta dei minimi quadrati

Il modello della regressione lineare semplice dice che la risposta y_i è una retta più un errore: y_i = β₀ + β₁x_i + ε_i, con errori ε_i indipendenti, di media zero e varianza σ² costante. I parametri β₀ (intercetta) e β₁ (pendenza: di quanto cambia in media y per un'unità di x) sono incogniti.

Il metodo dei minimi quadrati sceglie la retta che minimizza la somma dei quadrati dei residui, SSE = Σ(y_i − b₀ − b₁x_i)², dove il residuo è la distanza verticale di ogni punto dalla retta. Annullando le derivate si ottiene b₁ = S_xy/S_xx con S_xy = Σ(x_i − x̄)(y_i − ȳ) e S_xx = Σ(x_i − x̄)², e b₀ = ȳ − b₁x̄: la retta passa sempre per il punto delle medie.

L'esempio è volutamente piccolo, perché si possa rifare a mano: x = 1, 2, 3, 4, 5 e y = 2,1; 3,9; 6,2; 7,8; 10,1. Le medie sono 3 e 6,02, S_xx = 10 e S_xy = 19,90, quindi b₁ = 1,99 e b₀ = 6,02 − 1,99·3 = 0,05. La stima della varianza dell'errore è s² = SSE/(n − 2): si divide per n − 2 e non per n perché due parametri sono stati stimati dai dati, come per la varianza campionaria dell'articolo sugli stimatori.

L'errore standard della pendenza è s/√S_xx. Il programma ripete i conti, scrive lo stesso problema in forma matriciale (equazioni normali X′Xb = X′y, risolte anche con la decomposizione QR e con np.linalg.lstsq) e confronta con scipy.stats.linregress e con statsmodels.

# Regressione lineare semplice: esercizio svolto con 5 punti, formule a mano, forma matriciale e controllo con scipy e statsmodels
import numpy as np
from scipy import stats
import statsmodels.api as sm

x = np.array([1, 2, 3, 4, 5.0])
y = np.array([2.1, 3.9, 6.2, 7.8, 10.1])
n = len(x)

Sxx = ((x - x.mean()) ** 2).sum()
Sxy = ((x - x.mean()) * (y - y.mean())).sum()
b1 = Sxy / Sxx                               # pendenza
b0 = y.mean() - b1 * x.mean()                # intercetta: la retta passa per (x medio, y medio)
e = y - (b0 + b1 * x)                        # residui
sse = (e ** 2).sum()
sst = ((y - y.mean()) ** 2).sum()
s2 = sse / (n - 2)                           # stima di sigma^2: n - 2 gradi di liberta' (due parametri stimati)
se1 = np.sqrt(s2 / Sxx)
t1 = b1 / se1
p1 = 2 * stats.t.sf(abs(t1), n - 2)
q = stats.t.ppf(0.975, n - 2)
print(f"Sxx = {Sxx:.1f}, Sxy = {Sxy:.2f}, pendenza = {b1:.3f}, intercetta = {b0:.3f}")
print(f"residui = {np.round(e, 3).tolist()}, somma dei residui nulla: {abs(e.sum()) < 1e-12}, SSE = {sse:.4f}, SST = {sst:.4f}")
print(f"R^2 = {1 - sse / sst:.5f}, sigma^2 stimata = {s2:.5f}")
print(f"errore standard della pendenza = {se1:.4f}, t = {t1:.2f} con {n - 2} g.d.l., p = {p1:.2e}")
print(f"IC 95% della pendenza = [{b1 - q * se1:.3f}, {b1 + q * se1:.3f}]")

# Forma matriciale: equazioni normali X'X b = X'y, decomposizione QR, lstsq
X = np.column_stack([np.ones(n), x])
b_normali = np.linalg.solve(X.T @ X, X.T @ y)
Q, R = np.linalg.qr(X)
b_qr = np.linalg.solve(R, Q.T @ y)
b_lstsq = np.linalg.lstsq(X, y, rcond=None)[0]
print(f"coefficienti: normali {np.round(b_normali, 6).tolist()}, QR {np.round(b_qr, 6).tolist()}, lstsq {np.round(b_lstsq, 6).tolist()}")

r = stats.linregress(x, y)
m = sm.OLS(y, X).fit()
print(f"scipy.linregress: pendenza {r.slope:.3f}, errore standard {r.stderr:.4f}, p = {r.pvalue:.2e}, R^2 = {r.rvalue ** 2:.5f}")
print(f"statsmodels:      pendenza {m.params[1]:.3f}, errore standard {m.bse[1]:.4f}, p = {m.pvalues[1]:.2e}, R^2 = {m.rsquared:.5f}")

# Output:
# Sxx = 10.0, Sxy = 19.90, pendenza = 1.990, intercetta = 0.050
# residui = [0.06, -0.13, 0.18, -0.21, 0.1], somma dei residui nulla: True, SSE = 0.1070, SST = 39.7080
# R^2 = 0.99731, sigma^2 stimata = 0.03567
# errore standard della pendenza = 0.0597, t = 33.32 con 3 g.d.l., p = 5.94e-05
# IC 95% della pendenza = [1.800, 2.180]
# coefficienti: normali [0.05, 1.99], QR [0.05, 1.99], lstsq [0.05, 1.99]
# scipy.linregress: pendenza 1.990, errore standard 0.0597, p = 5.94e-05, R^2 = 0.99731
# statsmodels:      pendenza 1.990, errore standard 0.0597, p = 5.94e-05, R^2 = 0.99731

I tre metodi matriciali danno gli stessi coefficienti (0,05 e 1,99), e anche SciPy e statsmodels coincidono con i conti a mano: errore standard 0,0597, t = 33,32 con 3 gradi di libertà, p = 5,94 · 10⁻⁵, intervallo al 95% [1,800; 2,180]. I residui sommano a zero, come devono quando c'è l'intercetta. L'indice R² = 1 − SSE/SST = 0,99731 dice che la retta spiega il 99,7% della variabilità totale delle y.

Per le equazioni normali conviene usare la decomposizione QR e non invertire X′X: il numero di condizione di X′X è il quadrato di quello di X, e con variabili quasi collineari si perdono cifre. Le formule del caso generale con p coefficienti, e la proprietà di proiezione ortogonale, sono nella lezione sulla regressione del corso MIT 18.650.

Che cosa garantiscono le assunzioni

Con errori di media zero, varianza costante e incorrelati, gli stimatori dei minimi quadrati sono i migliori fra quelli lineari e non distorti (teorema di Gauss-Markov). Per fare inferenza esatta serve di più: errori Normali.

Allora le stime sono Normali, (n − p)s²/σ² segue una chi-quadro con n − p gradi di libertà, ed è indipendente dalle stime; da qui il test t sul singolo coefficiente con n − p gradi di libertà, enunciato nella stessa lezione MIT. Per la retta semplice p = 2. Per testare tutti i coefficienti insieme si usa il test F (che per una sola variabile esplicativa è il quadrato del t).

Il valore pratico sta negli intervalli. Per un punto x₀ si possono costruire tre intervalli al 95%, tutti con il quantile t a n − 2 gradi di libertà. L'intervallo per la pendenza è b₁ ± t·s/√S_xx. L'intervallo per la media della risposta in x₀ è ŷ₀ ± t·s√h₀, con h₀ = 1/n + (x₀ − x̄)²/S_xx, la leva del punto: cresce allontanandosi dalla media delle x. L'intervallo di previsione per una nuova misura è ŷ₀ ± t·s√(1 + h₀): l'«1» in più è l'errore ε della misura nuova, che nessun numero di dati elimina.

Il programma simula 100.000 esperimenti con 30 punti fissi fra 0 e 10, pendenza vera 2, intercetta 1 ed errori Normali con deviazione standard 1,5, e per ciascuno costruisce i tre intervalli in x₀ = 9 con una nuova misura futura. Conta quante volte ciascun intervallo contiene il suo bersaglio. Aggiunge una quarta riga, quella dell'errore descritto in apertura: l'intervallo per la media usato come se fosse un intervallo di previsione. In coda confronta con get_prediction di statsmodels.

# Intervallo di confidenza per la media, intervallo di previsione e intervallo per la pendenza: copertura reale
import numpy as np
from scipy import stats
import statsmodels.api as sm

n, R = 30, 100_000
beta0, beta1, sigma = 1.0, 2.0, 1.5
rng = np.random.default_rng(11)
x = np.linspace(0, 10, n)                             # disegno fisso
x0 = 9.0                                              # punto in cui si vuole prevedere
Sxx = ((x - x.mean()) ** 2).sum()
q = stats.t.ppf(0.975, n - 2)

Y = beta0 + beta1 * x + sigma * rng.standard_normal((R, n))
xm = x.mean()
b1 = ((x - xm) * (Y - Y.mean(axis=1, keepdims=True))).sum(axis=1) / Sxx
b0 = Y.mean(axis=1) - b1 * xm
e = Y - b0[:, None] - b1[:, None] * x
s = np.sqrt((e ** 2).sum(axis=1) / (n - 2))
yhat0 = b0 + b1 * x0
h0 = 1 / n + (x0 - xm) ** 2 / Sxx                     # leva del punto x0
se_media = s * np.sqrt(h0)
se_prev = s * np.sqrt(1 + h0)
media_vera = beta0 + beta1 * x0
y_nuovo = media_vera + sigma * rng.standard_normal(R)  # un'osservazione futura per ciascun esperimento

print(f"n = {n}, x0 = {x0}, esperimenti = {R}; errore standard di una copertura: {np.sqrt(0.95 * 0.05 / R):.4f}")
print(f"IC della pendenza copre il valore vero:        {(abs(b1 - beta1) <= q * s / np.sqrt(Sxx)).mean():.4f}")
print(f"IC della media in x0 copre la media vera:      {(abs(yhat0 - media_vera) <= q * se_media).mean():.4f}")
print(f"IC della media usato come previsione:          {(abs(yhat0 - y_nuovo) <= q * se_media).mean():.4f}")
print(f"intervallo di previsione copre il valore nuovo: {(abs(yhat0 - y_nuovo) <= q * se_prev).mean():.4f}")
print(f"rapporto fra le semiampiezze, previsione/media: esatto {np.sqrt((1 + h0) / h0):.3f}, simulato {(se_prev / se_media).mean():.3f}")

# Controllo con statsmodels su un campione
m = sm.OLS(Y[0], sm.add_constant(x)).fit()
pr = m.get_prediction(np.array([[1.0, x0]])).summary_frame(alpha=0.05)
print(f"controllo statsmodels (primo esperimento): media in x0 {pr['mean'].iloc[0]:.4f} contro {yhat0[0]:.4f}; "
      f"semiampiezza IC {(pr['mean_ci_upper'].iloc[0] - pr['mean'].iloc[0]):.4f} contro {q * se_media[0]:.4f}; "
      f"semiampiezza previsione {(pr['obs_ci_upper'].iloc[0] - pr['mean'].iloc[0]):.4f} contro {q * se_prev[0]:.4f}")

# Output:
# n = 30, x0 = 9.0, esperimenti = 100000; errore standard di una copertura: 0.0007
# IC della pendenza copre il valore vero:        0.9504
# IC della media in x0 copre la media vera:      0.9506
# IC della media usato come previsione:          0.4453
# intervallo di previsione copre il valore nuovo: 0.9497
# rapporto fra le semiampiezze, previsione/media: esatto 3.425, simulato 3.425
# controllo statsmodels (primo esperimento): media in x0 17.9370 contro 17.9370; semiampiezza IC 0.7687 contro 0.7687; semiampiezza previsione 2.6326 contro 2.6326

I tre intervalli corretti coprono come promesso: 95,04% per la pendenza, 95,06% per la media, 94,97% per la previsione (errore standard 0,0007). L'intervallo per la media, usato come previsione, copre solo il 44,53% delle misure future: è troppo stretto, perché ignora la variabilità della singola misura. Le semiampiezze stanno nel rapporto √((1 + h₀)/h₀) = 3,425, sia con la formula sia nella simulazione: l'intervallo di previsione è più di tre volte più largo. Il confronto con statsmodels in coda (semiampiezze 0,7687 e 2,6326) dà gli stessi valori, quindi la formula implementata è quella della libreria.

Più variabili: tre cose che cambiano

Nella regressione multipla la risposta dipende da p − 1 variabili: y = β₀ + β₁x₁ + … + ε, e ogni coefficiente β_j è la variazione media di y per un'unità di x_j a parità delle altre variabili. Stima, test e intervalli sono gli stessi, con n − p gradi di libertà. Tre fenomeni nascono solo con più variabili, e il programma li misura.

Primo: aggiungere variabili alza R² anche se sono pura casualità. Con n = 30 e 10 variabili indipendenti da y, il valore atteso di R² è p/(n − 1) = 10/29, un risultato esatto (la legge di R² è una Beta). Il rimedio è l'R² corretto, 1 − (1 − R²)(n − 1)/(n − p − 1), che compensa i gradi di libertà persi.

Secondo: l'omissione di una variabile correlata con quelle incluse sposta il coefficiente. Se il modello vero ha x₁ e x₂ (coefficienti 2 e 3, correlazione 0,6 fra le due variabili standardizzate) e tu stimi solo x₁, la pendenza attesa non è 2 ma 2 + 3 · 0,6 = 3,8: assorbe l'effetto della variabile omessa. Terzo: la collinearità. Quando due variabili esplicative sono molto correlate (correlazione ρ), la deviazione standard di ciascun coefficiente cresce di un fattore √VIF, con VIF = 1/(1 − ρ²) (variance inflation factor): il modello continua a prevedere bene, ma non sa attribuire l'effetto a una variabile o all'altra.

# Regressione multipla: R^2 con variabili di solo rumore, distorsione da variabile omessa, collinearita'
import numpy as np
from scipy import stats

rng = np.random.default_rng(0)

def ols(X, y):
    """Coefficienti per QR; restituisce anche gli errori standard classici."""
    Q, R = np.linalg.qr(X)
    b = np.linalg.solve(R, Q.T @ y)
    e = y - X @ b
    s2 = e @ e / (len(y) - X.shape[1])
    Rinv = np.linalg.inv(R)
    return b, np.sqrt(s2 * np.diag(Rinv @ Rinv.T)), e

# 1) Aggiungere p variabili di puro rumore fa salire R^2 anche se non c'e' nessun legame
n, p, Rep = 30, 10, 20_000
r2 = np.empty(Rep); r2_agg = np.empty(Rep)
for i in range(Rep):
    X = np.column_stack([np.ones(n), rng.standard_normal((n, p))])
    y = rng.standard_normal(n)                              # y non dipende da nessuna X
    _, _, e = ols(X, y)
    r2[i] = 1 - e @ e / ((y - y.mean()) @ (y - y.mean()))
    r2_agg[i] = 1 - (1 - r2[i]) * (n - 1) / (n - p - 1)       # R^2 corretto per i gradi di liberta'
print(f"1) n = {n}, p = {p} variabili di rumore: R^2 medio simulato {r2.mean():.4f}, esatto p/(n-1) = {p / (n - 1):.4f}; "
      f"R^2 corretto medio {r2_agg.mean():.4f}")
print(f"   R^2 supera 0.5 nel {(r2 > 0.5).mean():.4f} dei casi (esatto, legge Beta: {stats.beta.sf(0.5, p / 2, (n - p - 1) / 2):.4f})")

# 2) Variabile omessa: y = 1 + 2 x1 + 3 x2 + e, con corr(x1, x2) = 0.6; si stima y su x1 da sola
n, rho, Rep = 200, 0.6, 20_000
b_corto = np.empty(Rep); b_pieno = np.empty(Rep)
for i in range(Rep):
    z = rng.standard_normal((n, 2))
    x1 = z[:, 0]; x2 = rho * z[:, 0] + np.sqrt(1 - rho ** 2) * z[:, 1]
    y = 1 + 2 * x1 + 3 * x2 + rng.standard_normal(n)
    b_corto[i] = ols(np.column_stack([np.ones(n), x1]), y)[0][1]
    b_pieno[i] = ols(np.column_stack([np.ones(n), x1, x2]), y)[0][1]
print(f"2) pendenza di x1 (vero 2): senza x2 media {b_corto.mean():.4f} (attesa 2 + 3*0.6 = {2 + 3 * rho:.4f}); con x2 media {b_pieno.mean():.4f}")

# 3) Collinearita': deviazione standard della pendenza di x1 al crescere della correlazione con x2
print("3) corr(x1, x2)  VIF = 1/(1-r^2)  dev. std. simulata  attesa 1/sqrt(n(1-r^2))")
for rho in (0.0, 0.9, 0.99):
    b = np.empty(5_000)
    for i in range(b.size):
        z = rng.standard_normal((n, 2))
        x1 = z[:, 0]; x2 = rho * z[:, 0] + np.sqrt(1 - rho ** 2) * z[:, 1]
        y = 1 + 2 * x1 + 3 * x2 + rng.standard_normal(n)
        b[i] = ols(np.column_stack([np.ones(n), x1, x2]), y)[0][1]
    print(f"   {rho:13.2f}  {1 / (1 - rho ** 2):12.2f}  {b.std():16.4f}  {1 / np.sqrt(n * (1 - rho ** 2)):.4f}")

# Output:
# 1) n = 30, p = 10 variabili di rumore: R^2 medio simulato 0.3442, esatto p/(n-1) = 0.3448; R^2 corretto medio -0.0010
#    R^2 supera 0.5 nel 0.1080 dei casi (esatto, legge Beta: 0.1097)
# 2) pendenza di x1 (vero 2): senza x2 media 3.8013 (attesa 2 + 3*0.6 = 3.8000); con x2 media 1.9992
# 3) corr(x1, x2)  VIF = 1/(1-r^2)  dev. std. simulata  attesa 1/sqrt(n(1-r^2))
#             0.00          1.00            0.0717  0.0707
#             0.90          5.26            0.1656  0.1622
#             0.99         50.25            0.5109  0.5013

Con 10 variabili di rumore, R² vale in media 0,3442 (esatto 0,3448) e supera 0,5 nel 10,8% dei casi (esatto 0,1097): un modello senza alcun legame con la risposta sembra spiegare un terzo della varianza. L'R² corretto ha media −0,0010, cioè zero. Omettendo x₂, la pendenza di x₁ ha media 3,8013, contro 3,8 attesi; includendola, 1,9992, vicino al vero 2.

Con correlazione 0,90 il VIF è 5,26 e la deviazione standard del coefficiente passa da 0,0717 a 0,1656; con 0,99 il VIF è 50,25 e la deviazione standard è 0,5109, 7 volte maggiore. Un coefficiente non significativo in presenza di forte collinearità non significa che la variabile non conti: significa che i dati non bastano a separarla dall'altra.

Quando la varianza dell'errore non è costante

La seconda assunzione, varianza costante (omoschedasticità), è quella che si rompe più spesso: l'errore di una misura di tempo, di prezzo o di reddito tende a crescere con il livello. I coefficienti stimati restano non distorti, ma l'errore standard classico s√(X′X)⁻¹_jj sbaglia, e con esso intervalli e p-value. Gli errori standard robusti all'eteroschedasticità di White, nella variante HC3, stimano la varianza con (X′X)⁻¹ [Σ e_i²/(1 − h_i)² x_ix_i′] (X′X)⁻¹, dove h_i è la leva dell'osservazione i; statsmodels li offre con cov_type="HC3", descritti nella documentazione di RegressionResults.

Il programma simula 20.000 esperimenti con 100 punti in [0, 1], pendenza vera 2, e una deviazione standard dell'errore che cresce da 0,2 a 3,2 (σ_i = 0,2 + 3x_i²). Con un disegno fisso la vera deviazione standard della pendenza si calcola esattamente con la formula a «sandwich» con le varianze vere; il programma confronta la formula con la simulazione, poi misura quanto copre l'intervallo al 95% costruito con l'errore standard classico e con HC3. Verifica la propria implementazione di HC3 contro statsmodels.

# Quando la varianza dell'errore cresce con x: errori standard classici contro errori robusti HC3
import numpy as np
import statsmodels.api as sm
from scipy import stats

n, R = 100, 20_000
rng = np.random.default_rng(99)
x = np.linspace(0, 1, n)                              # disegno fisso
X = np.column_stack([np.ones(n), x])
sigma_i = 0.2 + 3.0 * x ** 2                          # deviazione standard dell'errore: da 0,2 a 3,2
XtXi = np.linalg.inv(X.T @ X)
H = X @ XtXi @ X.T
h = np.diag(H)

# Deviazione standard vera dello stimatore della pendenza (formula esatta, fissato il disegno)
cov_vera = XtXi @ (X.T * sigma_i ** 2) @ X @ XtXi
sd_vera = np.sqrt(cov_vera[1, 1])

Y = 1.0 + 2.0 * x + sigma_i * rng.standard_normal((R, n))
B = Y @ X @ XtXi                                      # coefficienti di tutti gli esperimenti, forma (R, 2)
E = Y - B @ X.T
s2 = (E ** 2).sum(axis=1) / (n - 2)
se_classico = np.sqrt(s2 * XtXi[1, 1])
# HC3: (X'X)^-1 [sum e_i^2 / (1-h_i)^2 x_i x_i'] (X'X)^-1
w = E ** 2 / (1 - h) ** 2
meat = np.einsum('ri,ij,ik->rjk', w, X, X)
se_hc3 = np.sqrt(np.einsum('a,rab,b->r', XtXi[1], meat, XtXi[:, 1]))

q = stats.t.ppf(0.975, n - 2)
cop = lambda se: (abs(B[:, 1] - 2.0) <= q * se).mean()
print(f"dev. std. vera della pendenza (formula) = {sd_vera:.4f}; simulata = {B[:, 1].std():.4f}")
print(f"errore standard classico: media {se_classico.mean():.4f}  -> IC al 95% copre {cop(se_classico):.4f}")
print(f"errore standard HC3:      media {se_hc3.mean():.4f}  -> IC al 95% copre {cop(se_hc3):.4f}")

# Controllo di HC3 con statsmodels sul primo esperimento
m = sm.OLS(Y[0], X).fit(cov_type="HC3")
print(f"controllo statsmodels (primo esperimento): HC3 {m.bse[1]:.6f} contro {se_hc3[0]:.6f}")

# Output:
# dev. std. vera della pendenza (formula) = 0.6326; simulata = 0.6378
# errore standard classico: media 0.5116  -> IC al 95% copre 0.8852
# errore standard HC3:      media 0.6316  -> IC al 95% copre 0.9432
# controllo statsmodels (primo esperimento): HC3 0.469747 contro 0.469747

La deviazione standard vera della pendenza è 0,6326 (simulata 0,6378). L'errore standard classico vale in media 0,5116, cioè sottostima del 19%, e l'intervallo che ne risulta copre l'88,52% delle volte invece del 95%. HC3 vale in media 0,6316, vicino al valore vero, e l'intervallo copre il 94,32%. Il controllo con statsmodels coincide a sei cifre (0,469747) sul primo esperimento. Con una retta e 100 punti il costo di HC3 è nullo: se hai dubbi sull'omoschedasticità, usare l'errore robusto non sbaglia quando la varianza è costante e protegge quando non lo è, come in questo caso.

Quattro errori che si commettono davvero

• Presentare l'intervallo di confidenza della media come fascia di previsione: copre le misure future il 44,5% delle volte.
• Aggiungere variabili finché R² sale: con 10 variabili di rumore su 30 punti sale a 0,34 in media.
• Leggere un coefficiente senza chiedersi quali variabili correlate mancano: la pendenza di 2 diventa 3,8.
• Fidarsi degli errori standard classici quando l'errore cresce con x: l'intervallo «al 95%» copre l'88,5%.

I quattro errori sono errori di lettura: il modello fa quello che deve, e il guaio è che il numero stampato (95%, R², coefficiente) è letto come se valesse in condizioni che non valgono. Per ciascuno hai un numero misurato che dice di quanto: se il tuo problema ha una variabile omessa plausibile, una varianza che cresce o due variabili molto correlate, hai già un'idea dell'ordine di grandezza dell'errore.

Adesso sai costruire a mano e con QR la retta dei minimi quadrati, calcolare errori standard, test t e R², distinguere intervallo per la media e intervallo di previsione e verificarne la copertura con una simulazione, riconoscere gli effetti di variabili inutili, omesse e collineari, e passare a errori robusti HC3 quando la varianza non è costante. Il modello lineare è anche il punto di arrivo dell'ANOVA, che è una regressione con variabili di gruppo, e nell'articolo sulla stima bayesiana vedrai come una distribuzione a priori cambia il modo di stimare un parametro.

Nella stessa serie trovi anche Analisi dati esplorativa (EDA) con Python.