Regressão Polinomial Global: o exemplo eólico
Este material apresenta, de forma progressiva, duas grandes famílias de modelos de regressão para descrever como \(E[Y]\) varia com uma covariável contínua \(x\):
Começamos pelos modelos globais polinomiais, do grau zero (a média) até o grau 5, comparando-os formalmente por meio de anova(). Na sequencia serão tratados modelos por partes sobre o mesmo conjunto de dados.
Vamos considerar um conjunto de dados com a variável explicativa \(X\) sendo a velocidade do vento (m/s) e a resposta Y$Y, a energia gerada.
eolica <- read.csv2("http://www.leg.ufpr.br/~paulojus/CE092/dados/eolica.csv")
str(eolica)'data.frame': 25 obs. of 2 variables:
$ vento : num 5 6 3.4 2.7 10 9.7 9.55 3.05 8.15 6.2 ...
$ energia: num 1.58 1.82 1.06 0.5 2.24 ...with(eolica, plot(energia ~ vento, pch = 19, col = "gray30",
xlab = "velocidade do vento (m/s)", ylab = "energia gerada"))
São \(n=25\) pares (velocidade do vento, energia gerada). A relação é crescente, mas a forma exata dessa curva, se uma reta, alguma curvatura suave, ou algo que satura em ventos altos é justamente o que os modelos a seguir tentam capturar e modelar. A seguir mostramos como ajustar diferentes modelos, comparar os ajustes e visualizar predições.
Em R vamos criar um objeto para armazenar as predições dos diferentes modelos.
preds <- data.frame(vento = seq(2, 11, length = 201))Um modelo de regressão polinomial de grau \(p\) especifica
\[E[Y_i] = \beta_0 + \beta_1 x_i + \beta_2 x_i^2 + \cdots + \beta_p x_i^p, \qquad i = 1, \ldots, n,\]
com \(\beta_0, \beta_1, \ldots, \beta_p\) estimados por mínimos quadrados. O grau \(p\) controla a flexibilidade do modelo: quanto maior \(p\), mais a curva pode se dobrar para acompanhar os dados. O caso \(p=0\) é o modelo mais simples de todos. Define um modelo sem nenhuma covariável, apenas estimando a média:
\[E[Y_i] = \beta_0.\]
Os seis modelos considerados aqui são:
| Grau \(p\) | \(E[Y]\) | Interpretação |
|---|---|---|
| 0 | \(\beta_0\) | uma única média |
| 1 | \(\beta_0+\beta_1 x\) | reta |
| 2 | \(\beta_0+\beta_1 x+\beta_2 x^2\) | parábola |
| 3 | \(\beta_0+\beta_1 x+\beta_2 x^2+\beta_3 x^3\) | cúbica |
| 4 | \(\beta_0+\cdots+\beta_4 x^4\) | quártica |
| 5 | \(\beta_0+\cdots+\beta_5 x^5\) | quíntica |
Cada modelo desta tabela é um caso particular do seguinte: o modelo de grau \(p\) é obtido do modelo de grau \(p+1\) impondo \(\beta_{p+1}=0\). Essa relação de aninhamento é o que torna válida a comparação sequencial via anova() feita adiante.
graus <- 0:5
modpol <- vector("list", length(graus))
names(modpol) <- paste0("grau", graus)
modpol$grau0 <- lm(energia ~ 1, data = eolica)
for (p in 1:5) {
## o grau "p" é inserido como valor literal na fórmula (e não como o símbolo "p")
## para que predict(newdata = ...) reavalie poly() corretamente mais adiante
form_p <- as.formula(sprintf("energia ~ poly(vento, %d, raw = TRUE)", p))
modpol[[paste0("grau", p)]] <- lm(form_p, data = eolica)
}
sapply(modpol, function(m) round(coef(m), 4))$grau0
(Intercept)
1.6096
$grau1
(Intercept) poly(vento, 1, raw = TRUE)
0.1309 0.2411
$grau2
(Intercept) poly(vento, 2, raw = TRUE)1
-1.1559 0.7229
poly(vento, 2, raw = TRUE)2
-0.0381
$grau3
(Intercept) poly(vento, 3, raw = TRUE)1
-2.3586 1.4269
poly(vento, 3, raw = TRUE)2 poly(vento, 3, raw = TRUE)3
-0.1604 0.0065
$grau4
(Intercept) poly(vento, 4, raw = TRUE)1
-4.7152 3.2773
poly(vento, 4, raw = TRUE)2 poly(vento, 4, raw = TRUE)3
-0.6621 0.0628
poly(vento, 4, raw = TRUE)4
-0.0022
$grau5
(Intercept) poly(vento, 5, raw = TRUE)1
-6.4744 5.0280
poly(vento, 5, raw = TRUE)2 poly(vento, 5, raw = TRUE)3
-1.3154 0.1775
poly(vento, 5, raw = TRUE)4 poly(vento, 5, raw = TRUE)5
-0.0118 0.0003 Por exemplo, os modelos de grau 1 a 3, com coeficientes estimados arredondados são:
\[\begin{align*} \text{(grau 1:)} \qquad \hat E[Y] &= 0.131 + 0.241\,x \\ \text{(grau 2:)} \qquad \hat E[Y] &= -1.156 + 0.723\,x -0.038\,x^2 \\ \text{(grau 3:)} \qquad \hat E[Y] &= -2.359 + 1.427\,x -0.16\,x^2 + 0.0065\,x^3 \end{align*}\]
Aqui os polinômios foram parametrizados na forma crua (raw = TRUE), isto é, \(X=[1,\,x,\,x^2,\,\ldots,\,x^p]\), exatamente como na tabela de fórmulas acima — o que torna os coeficientes diretamente legíveis na expressão de \(E[Y]\). O preço dessa legibilidade é que as colunas de \(X\) ficam fortemente correlacionadas entre si para \(p\) grande (multicolinearidade), o que não afeta o ajuste (valores preditos) mas torna os coeficientes individuais instáveis e de interpretação isolada pouco útil. A alternativa numérica — poly(x, p) sem raw = TRUE, que usa uma base ortogonal. Essa base produz exatamente os mesmos valores ajustados e o mesmo resultado de anova(), apenas com coeficientes reparametrizados.
cores <- c("gray50", "steelblue", "darkgreen", "tomato", "purple", "orange")
with(eolica, plot(energia ~ vento, pch = 19, col = "gray30",
xlab = "velocidade do vento (m/s)", ylab = "energia gerada"))
for (i in seq_along(modpol)) {
preds[[names(modpol)[i]]] <- predict(modpol[[i]], newdata = preds)
lines(preds$vento, preds[[names(modpol)[i]]], col = cores[i], lwd = 2)
}
legend("topleft", legend = paste("grau", graus), col = cores, lwd = 2, bty = "n", ncol = 2)
Uma forma útil de expor o comportamento de cada grau fora do intervalo observado dos dados (extrapolação) é ampliar a janela de \(x\):
preds_ext <- data.frame(vento = seq(0, 13, length = 301))
with(eolica, plot(energia ~ vento, pch = 19, col = "gray30",
xlim = c(0, 13), ylim = c(-2, 5),
xlab = "velocidade do vento (m/s)", ylab = "energia gerada"))
for (i in seq_along(modpol)) {
lines(preds_ext$vento, predict(modpol[[i]], newdata = preds_ext), col = cores[i], lwd = 2)
}
abline(v = range(eolica$vento), lty = 3, col = "gray60")
legend("topleft", legend = paste("grau", graus), col = cores, lwd = 2, bty = "n", ncol = 2)
Dentro do intervalo observado (linhas verticais pontilhadas) os polinômios de graus mais altos praticamente coincidem; enquanto que fora dele, divergem rapidamente. Isto ilustra o comportamento espúrio de polinômios globais de grau alto longe dos dados observados (ou mesmo entre eles, quando há poucos pontos). O comportamento errático dos polinômios globais de alto grau também pode ocorrer em situações nas quais há “gaps” entre valores de \(x\). Isso é um argumento importante a favor dos modelos por partes que virão a seguir: a flexibilidade local não exige extrapolação instável.
anova()Como os seis modelos ajustados até aqui são aninhados (Seção 3.1), anova() aplicado à sequência produz uma tabela de testes F sequenciais: cada linha testa se acrescentar o termo de grau \(p\) ao modelo de grau \(p-1\) reduz significativamente a soma de quadrados dos resíduos (RSS).
with(modpol, anova(grau0, grau1, grau2, grau3, grau4, grau5))Analysis of Variance Table
Model 1: energia ~ 1
Model 2: energia ~ poly(vento, 1, raw = TRUE)
Model 3: energia ~ poly(vento, 2, raw = TRUE)
Model 4: energia ~ poly(vento, 3, raw = TRUE)
Model 5: energia ~ poly(vento, 4, raw = TRUE)
Model 6: energia ~ poly(vento, 5, raw = TRUE)
Res.Df RSS Df Sum of Sq F Pr(>F)
1 24 10.2112
2 23 1.2816 1 8.9296 928.2858 < 2.2e-16 ***
3 22 0.3311 1 0.9505 98.8097 5.795e-09 ***
4 21 0.2354 1 0.0956 9.9432 0.005234 **
5 20 0.1862 1 0.0492 5.1196 0.035566 *
6 19 0.1828 1 0.0034 0.3545 0.558592
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Complementarmente, geramos uma tabela de \(R^2\), \(R^2\) ajustado, AIC e BIC que permite comparar todos os modelos entre si, não apenas pares consecutivos:
data.frame(
grau = graus,
gl = sapply(modpol, function(m) length(coef(m))),
R2 = sapply(modpol, function(m) summary(m)$r.squared),
R2ajustado = sapply(modpol, function(m) summary(m)$adj.r.squared),
AIC = sapply(modpol, AIC),
BIC = sapply(modpol, BIC)
) grau gl R2 R2ajustado AIC BIC
grau0 0 1 0.0000000 0.0000000 52.56213 54.999882
grau1 1 2 0.8744932 0.8690364 2.67724 6.333867
grau2 2 3 0.9675771 0.9646295 -29.16010 -24.284596
grau3 3 4 0.9769441 0.9736504 -35.68372 -29.589340
grau4 4 5 0.9817670 0.9781205 -39.55098 -32.237727
grau5 5 6 0.9821010 0.9773908 -38.01316 -29.481024Leitura dos resultados. O salto de grau 0 para grau 1 é enorme (RSS cai de 10.21 para 1.28, \(R^2\) passa de 0 para 0.874). Isto indica que a maior parte da variação de energia já é explicada por uma reta. Os graus 2, 3 e 4 continuam significativos a 5% na tabela sequencial (cada termo adicional ainda reduz o RSS mais do que se esperaria por acaso), mas com ganhos de \(R^2\) cada vez menores; o termo de grau 5 deixa de ser significativo (Pr(>F) bem acima de 0.05) e AIC/BIC pioram de grau 4 para grau 5. A interpretação é de que o grau 5 está ajustando ruído, não sinal. Tanto AIC quanto BIC preferem o grau 4 entre os seis candidatos.
Vale notar a diferença de propósito entre as duas tabelas. A anova() sequencial responde à pergunta: “vale a pena adicionar este termo ao modelo imediatamente anterior?” Cada linha da tabela é então um teste de hipótese, com os problemas usuais de testes múltiplos em sequência (a taxa de erro tipo I se acumula ao longo dos 5 testes). Já os critérios AIC/BIC comparam os modelos diretamente por um critério de ajuste penalizado pela complexidade, sem depender de uma ordem específica de comparação. Os dois caminhos aqui apontam para conclusões parecidas (grau 4 como um limite razoável, grau 5 desnecessário), mas isso não é garantido em geral — é um bom motivo para reportar mais de um critério. Métodos de regularização (não explorados aqui) também podem ser avaliados para comparar os modelos.
O resultado obtido aqui selecionando um grau relativamente alto (4) pelos critérios de significância/AIC num conjunto de apenas 25 pontos deve se visto com cautela. Com poucos dados, é fácil que um polinômio de grau alto capture particularidades da amostra em vez de um padrão real. Tópicos seguintes sobre modelos por partes, oferecem uma via alternativa de ganhar flexibilidade sem depender de potências altas de \(x\) em todo o domínio.