Regressione polinomiale

La regressione polinomiale utilizza lo stesso metodo della regressione lineare, ma assume che la funzione che meglio descrive l'andamento dei dati non sia una retta, ma un polinomio. Quindi è adatta quando lo scatterplot di una relazione bivariata, ad esempio, mostra una forma diversa da quella della retta, ad esempio una curva, come nell'esempio che segue.

# per i grafici e il piping
library( tidyverse)

La relazione fra le due variabili non sembra ben rappresentata da una retta (anche se potrebbe dipendere dagli outliers, visto che i punti in coda sono pochi).

$$\hat Y = a + b_1X + b_2X^2$$

lm(mpg ~ hp + I(hp^2), data = mtcars) %>% 
summary()

equivale a:

lm(mpg ~ poly(hp, 2, raw=TRUE), data = mtcars) %>% 
summary()

## Call:
## lm(formula = mpg ~ hp + I(hp^2), data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.5512 -1.6027 -0.6977  1.5509  8.7213 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.041e+01  2.741e+00  14.744 5.23e-15 ***
## hp          -2.133e-01  3.488e-02  -6.115 1.16e-06 ***
## I(hp^2)      4.208e-04  9.844e-05   4.275 0.000189 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.077 on 29 degrees of freedom
## Multiple R-squared:  0.7561,	Adjusted R-squared:  0.7393 
## F-statistic: 44.95 on 2 and 29 DF,  p-value: 1.301e-09

$$\hat {\text{mpg}} = 4.041e{+01} -2.13e{-01} hp + 4.20e{-04} hp^2$$

Ovvero:

x = mtcars$hp
fitted.y = 4.041e+01 + (-2.133e-01 * x) + (4.208e-04 * x^2)

Per avere polinomi di grado superiore, è sufficiente scrivere:

• poly(hp, 3, raw=TRUE); oppure
• hp + I(hp^2) + I(hp^3)

La prima forma è quindi più sintetica e comoda da scrivere.

Se non specifichiamo raw = TRUE, con poly(x, n) otteniamo una regressione ortogonale, in cui cioè i coefficienti sono indipendenti l'uno dall'altro:

fit1 <- lm(mpg ~ poly(hp, 2), data = mtcars)
summary(fit1)

## Call:
## lm(formula = mpg ~ poly(hp, 2), data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.5512 -1.6027 -0.6977  1.5509  8.7213 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    20.091      0.544  36.931  < 2e-16 ***
## poly(hp, 2)1  -26.046      3.077  -8.464 2.51e-09 ***
## poly(hp, 2)2   13.155      3.077   4.275 0.000189 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.077 on 29 degrees of freedom
## Multiple R-squared:  0.7561,	Adjusted R-squared:  0.7393 
## F-statistic: 44.95 on 2 and 29 DF,  p-value: 1.301e-09

Vedi:

• Basi dei grafici con ggplot;
• il grafico alla pagina Regressione lineare bivariata

Per rappresentare la curva di interpolazione:

ggplot(mtcars, aes(hp, mpg)) +
  geom_point(alpha = 0.7) +
  geom_smooth(method = "lm", se = F) +                 # retta dei minimi quadrati
  geom_smooth(method = "lm", 
              formula = y ~ poly(x,2), col = "red")    # geom_smooth con la formula

E' possibile scaricare ed eseguire lo script dell'esempio:

Es-regr-poli.R

library(tidyverse)
 
# con I()
lm(mpg ~ hp + I(hp^2), data = mtcars) %>% 
summary()
 
# con poly()
lm(mpg ~ poly(hp, 2, raw=TRUE), data = mtcars) %>% 
summary()
 
# ortogonale, salvando il modello 'fit'
fit <- lm(mpg ~ poly(hp, 2), data = mtcars)
summary(fit)
 
 
# grafico della curva
ggplot(mtcars, aes(hp, mpg)) +
  geom_point(alpha = 0.7) +
  geom_smooth(method = "lm", se = F) +                 # retta dei minimi quadrati
  geom_smooth(method = "lm", 
              formula = y ~ poly(x,2), col = "red")    # geom_smooth con la formula