Regressió lineal i diagnòstic: mínims quadrats, predicció i límits del model

La regressió lineal és el model estadístic més utilitzat del món i també el més mal utilitzat. Matemàticament és una projecció ortogonal i prou; la dificultat real està a saber què autoritza a concloure i què no. Aquesta pàgina fa les dues coses: el càlcul complet i, després, la llista de les hipòtesis que el sostenen i què passa quan fallen.

El model i les equacions normals

El model lineal s'escriu $y = X\beta + \varepsilon$, amb la matriu de disseny $X$ de rang columna complet. Les hipòtesis bàsiques són que l'error té esperança condicional nul·la i, per a les fórmules habituals d'inferència, que la seva variància condicional és constant i igual a $\sigma^{2}$ (homoscedasticitat) i que els errors no estan correlacionats. La normalitat dels errors és una hipòtesi addicional, que no cal per a l'estimació però sí per a la inferència exacta amb mostres petites.

Mínims quadrats minimitza la suma de quadrats dels residus. Derivant respecte de $\beta$ i igualant a zero s'obtenen les equacions normals i, formalment,

$$X^{T}X\,\widehat\beta = X^{T}y, \qquad \widehat\beta = (X^{T}X)^{-1}X^{T}y.$$

Geomètricament, el vector ajustat és la projecció ortogonal de $y$ sobre l'espai generat per les columnes de $X$, i el vector de residus és el que queda perpendicular. D'aquí surten totes les propietats que veuràs.

Regressió simple

Amb un sol predictor i terme independent,

$$\widehat\beta_1 = \frac{\sum_i (x_i-\bar x)(y_i-\bar y)}{\sum_i (x_i-\bar x)^{2}} = \frac{S_{xy}}{S_{xx}}, \qquad \widehat\beta_0 = \bar y - \widehat\beta_1\,\bar x.$$

Cal que els valors del predictor no siguin tots iguals, perquè si no el denominador s'anul·la. De la segona fórmula es dedueix immediatament que la recta passa pel centre de gravetat dels punts. I de l'ortogonalitat amb la columna de uns surt que els residus sumen zero.

Incertesa dels coeficients i predicció

Sota el model, l'estimador és no esbiaixat i la seva matriu de variàncies condicional és $\sigma^{2}(X^{T}X)^{-1}$. La variància de l'error s'estima amb

$$s^{2} = \frac{\text{SSE}}{n-p},$$

on $p$ és el nombre de paràmetres. Amb errors normals, cada coeficient estandarditzat pel seu error estàndard segueix una $t$ amb $n-p$ graus de llibertat, i d'aquí surten els contrastos i els intervals dels coeficients.

Per a la predicció hi ha dos intervals diferents i confondre'ls és un error clàssic. En regressió simple, per a una abscissa $x_0$, l'interval per a la mitjana de la resposta té semiamplada

$$t_{n-2,\,1-\alpha/2}\; s\,\sqrt{\frac{1}{n} + \frac{(x_0-\bar x)^{2}}{S_{xx}}},$$

mentre que per a una nova observació individual cal sumar un $1$ dins de l'arrel. La diferència és conceptual: el primer només recull la incertesa de la recta ajustada, i el segon hi afegeix l'error aleatori propi de l'observació nova, que no desapareix per molta mostra que tinguis.

Fixa't també en el terme $(x_0-\bar x)^{2}$: tots dos intervals s'eixamplen a mesura que t'allunyes del centre de les dades. La recta està molt millor determinada al mig del núvol de punts que als extrems.

Ajust, residus i límits

Amb terme independent i resposta no constant, el coeficient de determinació és

$$R^{2} = 1 - \frac{\text{SSE}}{\text{SST}}.$$

Mesura la proporció de la variabilitat mostral explicada per l'ajust. Tres avisos imprescindibles: no demostra causalitat, no garanteix bona predicció fora de la mostra, i mai no disminueix en afegir predictors, encara que siguin soroll pur. Per comparar models amb nombres diferents de variables cal una mesura penalitzada o una validació sobre dades no utilitzades per ajustar.

El diagnòstic no és opcional. Cal mirar els residus contra els valors ajustats (per detectar curvatura no modelada i heteroscedasticitat), la dependència si les dades tenen ordre temporal, les observacions influents (un sol punt allunyat pot dominar l'ajust) i la multicol·linearitat entre predictors, que infla enormement els errors estàndard sense afectar la predicció.

Sobre l'heteroscedasticitat convé una precisió: invalida els errors estàndard homoscedàstics i per tant els contrastos i els intervals, però no trenca necessàriament la manca de biaix dels coeficients, que només depèn de l'esperança condicional nul·la. Per això la solució habitual és corregir els errors estàndard, no llençar el model.

I finalment, extrapolar fora del rang observat de les dades depèn completament de suposar que la forma lineal continua valent allà on no has mirat. És l'operació més arriscada de tota la regressió.

Exercici resolt: per què els residus sumen zero

Per què una recta de regressió amb terme independent té residus de suma zero?

Les equacions normals es poden escriure com $X^{T}(y - X\widehat\beta) = X^{T}e = 0$, és a dir, el vector de residus és ortogonal a totes les columnes de la matriu de disseny.

Quan el model inclou terme independent, una d'aquestes columnes és el vector de uns. L'ortogonalitat amb ell diu exactament que la suma dels residus és zero.

Conseqüència pràctica: si el teu model no inclou terme independent, aquesta propietat desapareix, els residus poden tenir mitjana no nul·la i la descomposició habitual de sumes de quadrats deixa de ser vàlida. Per això el coeficient de determinació d'un model sense intercepció no es pot interpretar de la manera habitual, i alguns programes fins i tot el calculen amb una fórmula diferent.

Exercici resolt: els dos intervals de predicció

Amb $n = 50$, desviació residual $s = 2$, mitjana del predictor $\bar x = 10$ i $S_{xx} = 200$, compara els dos intervals a $x_0 = 12$.

L'arrel per a la mitjana de la resposta val

$$\sqrt{\frac{1}{50} + \frac{(12-10)^{2}}{200}} = \sqrt{0{,}02 + 0{,}02} = \sqrt{0{,}04} = 0{,}2,$$

i amb $s = 2$ i un quantil aproximat de $2{,}01$ la semiamplada és uns $0{,}80$.

Per a una nova observació, l'arrel és $\sqrt{1{,}04} \approx 1{,}0198$ i la semiamplada puja a uns $4{,}10$.

Cinc vegades més ample. El motiu és que, amb cinquanta observacions, la recta ja està molt ben determinada (el terme $0{,}04$ és petit), però l'error individual d'una observació nova continua tenint desviació $2$ i això no millora mai. Si t'interessa on cau la mitjana, el primer interval; si t'interessa on caurà el pròxim punt, el segon. Donar el primer quan et demanaven el segon és subestimar enormement la incertesa.