Ci sono dei casi nei quali dobbiamo riuscire a coniugare la razionalità statistica teorica con le evidenze imposte dai dati empirici, vediamo un esempio.
L'emoglobina, la proteina che trasporta l'ossigeno dai polmoni ai tessuti, è contenuta esclusivamente nei globuli rossi (o eritrociti) del sangue e ci si attende quindi una relazione di proporzionalità lineare tra la concentrazione dell'emoglobina e la concentrazione degli eritrociti. Inoltre nel caso limite di una concentrazione degli eritrociti ridotta a zero (0) anche l'emoglobina deve avere una concentrazione pari a zero e quindi la relazione di proporzionalità lineare tra le due grandezze dovrebbe essere rappresentata da una retta che passa per l'origine degli assi (x,y = 0,0).
Nel set di dati ais contenuto nel pacchetto DAAG [1] abbiamo 202 valori di concentrazione degli eritrociti, espressa in migliaia di miliardi per litro di sangue (10¹²/L) – o milioni per microlitro di sangue (10⁶/µL), le due espressioni sono numericamente equivalenti – e i corrispondenti valori di emoglobina, la cui concentrazione è espressa in grammi per decilitro di sangue (g/dL).
Per analizzare i dati impieghiamo inizialmente la regressione lineare ordinaria.
# REGRESSIONE LINEARE
#
library (DAAG) # carica il pacchetto che include il set di dati ais
#
x <- ais$rcc # riporta gli eritrociti in ascisse
y <- ais$hg # riporta l'emoglobina in ordinate
reglin <- lm(y ~ x) # calcola la regressione lineare ordinaria
#
summary(reglin) # risultati della regressione
#
Dopo avere caricato il pacchetto che include il set di dati ais per chiarezza semantica gli eritrociti sono ridenominati x e l'emoglobia è ridenominata y.
Quindi viene calcolata con la funzione lm() la retta di regressione che esprime la y in funzione della x (y ~ x), e con la funzione summary() sono riportate le statistiche della regressione.
> summary(reglin) # risultati della regressione
Call:
lm(formula = y ~ x, data = ais)
Residuals:
Min 1Q Median 3Q Max
-1.66939 -0.44231 0.03338 0.44605 1.50594
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.08973 0.45703 4.572 8.44e-06 ***
x 2.64412 0.09641 27.427 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.626 on 200 degrees of freedom
Multiple R-squared: 0.79, Adjusted R-squared: 0.7889
F-statistic: 752.2 on 1 and 200 DF, p-value: < 2.2e-16
Le cose che ci interessano principalmente in questa analisi numerica sono due:
→ la prima è che la retta di regressione y=2.08973+2.64412·x quando poniamo x=0 intercetta l'asse delle y al valore = 2.08973 e quindi non passa per l'origine degli assi, un fatto che cozza contro quanto atteso in base all'evidenza empirica;
→ la seconda è che il valore di R² – per definizione compreso tra 0 e 1 – che rappresenta la quota di varianza spiegata dalla retta di regressione, in questo caso è uguale a 0.79 (ci servirà tra poco come termine di paragone).
A questo punto aggiungiamo un grafico con i dati e con la retta di regressione.
#
plot(x, y, xlim=c(0,8), ylim=c(0,20), xaxs="i", yaxs="i", xlab="rcc", ylab="hg", main="Regressione lineare ordinaria e passante per l'origine degli assi", cex.main=0.9) # grafico dei dati
abline(reglin, col="blue") # retta di regressione
#
Con la funzione plot() viene realizzato il grafico impiegando, accanto agli usuali argomenti, gli argomenti xaxs="i" e yaxs="i" che allineano inizio e fine delle scale al riquadro grafico (gli argomenti di default sono xaxs="r" e yaxs="r", e lasciano il 6% di spazio aggiuntivo tra inizio e fine delle scale e riquadro grafico, si può fare una prova con questi valori per visualizzare la differenza).
Con la funzione abline() viene poi aggiunta al grafico la retta di regressione, che conferma il risultato nient'affatto soddisfacente.
A fronte della nostra certezza che in assenza di eritrociti deve essere assente anche l'emoglobina, ripetiamo i calcoli impiegando un modello di regressione che impone il passaggio della retta per l'origine degli assi.
#
reglin0 <- lm(y ~ 0 + x) # regressione lineare passante per x=0 e y=0
abline(reglin0, col="red") # retta di regressione passante per x=0 e y=0
summary(reglin0) # risultati della regressione
#
Per farlo impieghiamo nella funzione lm() l'argomento y ~ 0 + x quindi con la funzione abline() riportiamo sul grafico la corrispondente retta di regressione.
Infine validiamo dal punto di vista statistico i risultati della regressione impiegando la funzione summary().
> summary(reglin0) # risultati della regressione
Call:
lm(formula = y ~ 0 + x, data = ais)
Residuals:
Min 1Q Median 3Q Max
-2.18149 -0.45742 0.05747 0.46888 1.60372
Coefficients:
Estimate Std. Error t value Pr(>|t|)
x 3.08288 0.00974 316.5 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.6562 on 201 degrees of freedom
Multiple R-squared: 0.998, Adjusted R-squared: 0.998
F-statistic: 1.002e+05 on 1 and 201 DF, p-value: < 2.2e-16
La conclusione?
Impiegando un modello alternativo che, sulla base di un dato biologico certo, impone il passaggio della retta di regressione per l'origine degli assi:
→ le differenze residue tra i singoli punti e la retta di regressione (Residuals) diventano un poco superiori a quelle della regressione lineare ordinaria, ma si tratta di un fatto atteso in quanto questa, impiegando il metodo dei minimi quadrati, è per definizione quella che li minimizza;
→ il valore di R² – che con la regressione ordinaria era, lo ricordiamo, 0.79 – migliora aumentando a 0.998;
→ la retta di regressione passante per l'origine degli assi spiega il 99.8% della varianza osservata, un risultato che sfiora il massimo teorico e che conferma la validità del razionale e del modello alternativo impiegato.
----------
[1] Vedere il post Il set di dati ais nel quale trovate anche come caricare i dati senza impiegare il pacchetto DAAG.
Nessun commento:
Posta un commento