Tecniche di regressione polinomiale

Nel software statistico R esistono numerosi pacchetti e numerose funzioni per stimare un modello di regressione polinomiale lineare. Questo è però uno degli argomenti più oscuri, a mio avviso, nelle guide e nei tutorial di R, i quali trattano solo di sfuggita questo capitolo. Un altro problema che ho dovuto affrontare, con non poche difficoltà, è la rappresentazione grafica dei modelli trovati. Mi duole ammettere che R (il software statistico che personalmente preferisco) ha non poche limitazioni per quanto riguarda questo specifico argomento.
Vi chiedo scusa per questo sfogo personale, ma ho perso quasi un pomeriggio per giungere alle conclusioni che ora vi presenterò :-)

Supponiamo di voler creare un polinomio che possa approssimare al meglio il seguente dataset, relativo alla popolazione di una certa città italiana nell'arco di 10 anni. In tabella sono riassunti i dati:



Per prima cosa importiamo i dati in R. Se questa è una tabella di Excel, possiamo utilizzare la funzione read.delim("clipboard"), come discusso in un post precedente. Otterremo così un dataframe contenente tutti i valori (questa è la scelta che io consiglio, cioè quella di creare un dataframe). In alternativa possiamo inserire i dati nel modo consueto:

Anno <- c(1959, 1960, 1961, 1962, 1963, 1964, 1965, 1966, 1967, 1968, 1969)
Popolazione <- c(4835, 4970, 5085, 5160, 5310, 5260, 5235, 5255, 5235, 5210, 5175)

A questo punto creiamo il dataframe di nome data:

data <- data.frame(Anno, Popolazione)
data
   Anno Popolazione
1  1959        4835
2  1960        4970
3  1961        5085
4  1962        5160
5  1963        5310
6  1964        5260
7  1965        5235
8  1966        5255
9  1967        5235
10 1968        5210
11 1969        5175

Io consiglio di procedere lavorando con i dataframe, ma questa scelta non è obbligatoria: potremmo anche lavorare con i vettori.

A questo punto può essere utile mettere in grafico questi valori, per osservare l'andamento e farci un'idea della forma del grafico finale. Per comodità modifichiamo la colonna Anno che contiene numeri molto elevati, creando un intorno dello zero, in questo modo:

data$Anno <- data$Anno - 1964
data
   Anno Popolazione
1    -5        4835
2    -4        4970
3    -3        5085
4    -2        5160
5    -1        5310
6     0        5260
7     1        5235
8     2        5255
9     3        5235
10    4        5210
11    5        5175

Richiamiamo il grafico così:

plot(data$Anno, data$Popolazione, type="b")



A questo punto possiamo cominciare con la ricerca di un modello polinomiale che approssimi adeguatamente i nostri dati. Innanzitutto specifichiamo che quello che vogliamo è un polinomio in funzione di X, ossia un raw polynomial in inglese, che è differente dall'orthogonal polynomial. Questa è una precisazione importante perchè i comandi e i risultati cambieranno nei due casi in R. Quindi noi cerchiamo una funzione di X del tipo:



A che grado del polinomio fermarci? Dipende dal grado di precisione che cerchiamo. Maggiore è il grado del polinomio, maggiore sarà la precisione del modello, ma maggiori saranno le difficoltà di calcolo; inoltre bisogna verificare la significatività dei coefficienti che vengono trovati. Ma veniamo subito al dunque.

In R per fittare un modello di regressione polinomiale (non ortogonale) esistono due metodi, tra loro identici. Supponiamo di cercare i valori dei coefficienti beta per un polinomio di 1° grado, quindi di 2° grado, quindi di 3° grado:

fit1 <- lm(data$Popolazione ~ data$Anno)
fit2 <- lm(data$Popolazione ~ data$Anno + I(data$Anno^2))
fit3 <- lm(data$Popolazione ~ data$Anno + I(data$Anno^2) + I(data$Anno^3))

Oppure possiamo scrivere più rapidamente per i polinomi di grado 2 e 3:

fit2b <- lm(data$Popolazione ~ poly(data$Anno, 2, raw=TRUE))
fit3b <- lm(data$Popolazione ~ poly(data$Anno, 3, raw=TRUE))

La funzione poly che vediamo applicata è utile qualora si desideri ottenere un polinomio di grado elevato, perchè evita di scrivere esplicitamente la formula. Se specifichiamo raw=TRUE, i due metodi visti forniscono lo stesso output; se invece non specifichiamo raw=TRUE, la funzione poly ci fornirà i valori dei parametri beta di un polinomio ortogonale, che è differente dalla formula generale che ho scritto sopra, sebbene i modelli siano entrambi efficaci. Rimando a testi di statistica per la differenza tra le due metodiche.

Analizziamo ora l'output.

> summary(fit2)

Call:
lm(formula = data$Popolazione ~ data$Anno + I(data$Anno^2))

Residuals:
    Min      1Q  Median      3Q     Max
-46.888 -18.834  -3.159   2.040  86.748

Coefficients:
               Estimate Std. Error t value Pr(>|t|)
(Intercept)    5263.159     17.655 298.110  < 2e-16 ***
data$Anno        29.318      3.696   7.933 4.64e-05 ***
I(data$Anno^2)  -10.589      1.323  -8.002 4.36e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 38.76 on 8 degrees of freedom
Multiple R-squared: 0.9407, Adjusted R-squared: 0.9259
F-statistic: 63.48 on 2 and 8 DF, p-value: 1.235e-05

L'output di summary(fit2b) è assolutamente identico. Abbiamo ottenuto i valori di beta0 (5263.159), beta1 (29.318) e beta2 (-10.589), i quali risultano essere tutt'e 3 significativi. L'equazione del polinomio di 2° grado del nostro modello è quindi:



Se cercassimo un polinomio di 3° grado, avremmo:

> summary(fit3)

Call:
lm(formula = data$Popolazione ~ data$Anno + I(data$Anno^2) +
    I(data$Anno^3))

Residuals:
    Min      1Q  Median      3Q     Max
-32.774 -14.802  -1.253   3.199  72.634

Coefficients:
                Estimate Std. Error t value Pr(>|t|)
(Intercept)    5263.1585    15.0667 349.324 4.16e-16 ***
data$Anno        14.3638     8.1282   1.767   0.1205
I(data$Anno^2)  -10.5886     1.1293  -9.376 3.27e-05 ***
I(data$Anno^3)    0.8401     0.4209   1.996   0.0861 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 33.08 on 7 degrees of freedom
Multiple R-squared: 0.9622, Adjusted R-squared: 0.946
F-statistic: 59.44 on 3 and 7 DF, p-value: 2.403e-05

E l'equazione sarebbe:



In quest'ultimo caso però i coefficienti beta1 e beta3 non sono significativi, quindi il modello migliore è il polinomio di 2° grado.

Il problema più grande è adesso quello di rappresentare graficamente quanto ottenuto. In R difatti non esiste (per quanto ne so io) una funzione per plottare i polinomi trovati fittando i dati. Dobbiamo quindi procedere con artefatti grafici comunque validi, ma alquanto laboriosi.

Innanzitutto mettiamo in grafico i valori, con il comando visto prima. Questa volta visualizzeremo solamente le linee e non i punti, per comodità grafica:

plot(data$Anno, data$Popolazione, type="l", lwd=3)

Ora aggiungiamo a questo grafico l'andamento del polinomio di 2° grado in questo modo:

points(data$Anno, predict(fit2), type="l", col="red", lwd=2)

In pratica la funzione predict() non fa altro che sostituire tutti i valori contenuti nella colonna data$Anno nel polinomio di 2° grado che abbiamo calcolato. Quindi le varie coppie di valori vengono unite da segmenti. Quello che otteniamo (vedi figura sotto) è quindi l'unione dei punti virtuali, ossia del valori y = f(x). Non è disegnato il grafico nel continuo, ma nel discreto. Con pochi valori, questo metodo è fortemente debilitante.



Proviamo ad aggiungere anche il grafico del polinomio di 3° grado:

points(data$Anno, predict(fit3), type="l", col="blue", lwd=2)



Come vedete i due modelli hanno andamenti molto simili.

Se volessimo invece ottenere il grafico continuo delle funzioni ottenute, dobbiamo procedere in questo modo. Innanzitutto si crea la funzione polinomiale trovata:

andamento2 <- function(x) fit2$coefficient[3]*x^2 + fit2$coefficient[2]*x + fit2$coefficient[1]

La funzione function(x) precede la formula che è quella del polinomio di 2° grado. fit2$coefficient[3] è beta2, il coefficient[2] è beta1 e coefficient[1] è beta0.

A questo punto mettiamo in grafico solo i punti, e poi visualizziamo la funzione creata con la sintassi curve(x):

plot(data$Anno, data$Popolazione, type="p", lwd=3)
andamento2 <- function(x) fit2$coefficient[3]*x^2 + fit2$coefficient[2]*x + fit2$coefficient[1]
curve(andamento2, col="red", lwd=2)

I punti però scompaiono, ma possiamo riposizionarli con il comando points:

points(data$Anno, data$Popolazione, type="p", lwd=3)

Una nota: è necessario seguire l'ordine dei comandi come ve l'ho descritto, altrimenti la funzione curve genera un grafico errato. Quindi riassumendo i comandi per ottenere la funzione continua, e i punti sperimentali sullo stesso grafico sono i seguenti:

plot(data$Anno, data$Popolazione, type="p", lwd=3)
andamento2 <- function(x) fit2$coefficient[3]*x^2 + fit2$coefficient[2]*x + fit2$coefficient[1]
curve(andamento2, col="red", lwd=2)
points(data$Anno, data$Popolazione, type="p", lwd=3)

Il grafico che otteniamo è il seguente:



Adesso cerchiamo il grafico del polinomio di 3° grado:

plot(data$Anno, data$Popolazione, type="p", lwd=3)
andamento3 <- function(x) fit3$coefficient[4]*x^3 + fit3$coefficient[3]*x^2 + fit3$coefficient[2]*x + fit3$coefficient[1]
curve(andamento3, col="red", lwd=2)
points(data$Anno, data$Popolazione, type="p", lwd=3)



Ribadisco quindi che per ottenere un output corretto è necessario seguire la serie di istruzioni viste finora nell'ordine con cui vi è stato presentato, altrimenti i grafici risulteranno alterati.
Nonostante queste piccole difficoltà, il risultato è pur sempre notevole.

E adesso un piccolo esercizio per voi :-D
Con la serie di dati che qui vi presento, scrivete nei commenti il grado del polinomio e il polinomio stesso che ho trovato per ottenere la curva rossa (dovrete andare a tentativi, così potrete esercitarvi con le istruzioni viste finora, fino a generare un grafico simile al mio!).

x <- c(1:18)
y <- c(7, 6, 4, 3, 4, 5, 4, 3, 2, 3, 4, 5, 6, 4, 2, 3, 5, 7)