Revisão: Regressão Linear
Este curso trata de algumas extensões do modelo de regressão linear: splines e regressão suavizada, modelos aditivos generalizados (GAM), modelos não lineares, modelos heterocedásticos, modelos de efeitos aleatórios e regressão para dados no intervalo unitário, dentre outras. Todas essas extensões partem do modelo de regressão linear clássico e relaxam alguma(s) de suas suposições. Por exemplo:
Por isso, antes de avançar, vale revisar tópicos de regressão linear simples e múltipla: especificação, estimação, propriedades dos resíduos e critérios de escolha de modelo. Essa revisão fixa a notação e o vocabulário que serão reaproveitados (e estendidos) em todos os tópicos seguintes.
Para \(i = 1, \dots, n\) observações de uma resposta \(y_i\) e uma covariável \(x_i\), o modelo de regressão linear simples é
\[ y_i = \beta_0 + \beta_1 x_i + e_i \; , \]
com as suposições usuais sobre o erro \(e_i\):
Note que \(\beta_0\) e \(\beta_1\) são parâmetros fixos e desconhecidos. A suposição de linearidade está em \(\beta_0 + \beta_1 x_i\), ou seja, o modelo é linear nos parâmetros (\(\beta\)’s), não necessariamente em \(x\). Poderíamos ter \(x_i^2\) ou \(\log(x_i)\) como covariável e o modelo continuaria sendo linear.
Os estimadores \(\hat\beta_0, \hat\beta_1\) minimizam a soma de quadrados dos resíduos
\[ SQ(\beta_0, \beta_1) = \sum_{i=1}^n \left( y_i - \beta_0 - \beta_1 x_i \right)^2 . \]
Derivando em relação a \(\beta_0\) e \(\beta_1\) e igualando a zero (equações normais), obtemos
\[\begin{align*} \hat\beta_1 &= \frac{\sum_{i=1}^n (x_i - \bar x)(y_i - \bar y)}{\sum_{i=1}^n (x_i - \bar x)^2} = \frac{\sum_{i=1}^n x_i y_i - \frac{(\sum_{i=1}^n x_i)(\sum_{i=1}^n y_i)}{n}}{\sum_{i=1}^n x_i^2 - \frac{(\sum_{i=1}^n x_i)^2}{n}} = \frac{S_{xy}}{S_{xx}} \\ \hat\beta_0 &= \bar y - \hat\beta_1 \bar x \; . \end{align*}\]
Sob as suposições 1–3, o estimador de mínimos quadrados é o melhor estimador linear não viesado (teorema de Gauss-Markov), com
\[ \text{Var}(\hat\beta_1) = \frac{\sigma^2}{S_{xx}} \; . \]
Um estimador não viesado de \(\sigma^2\) é dado pelo quadrado médio residual:
\[ \hat\sigma^2 = \frac{1}{n-2}\sum_{i=1}^n (y_i - \hat y_i)^2 = \frac{SQRes}{n-2}\; , \]
com \(n-2\) graus de liberdade que correspondem aos dois parâmetros (\(\beta_0\) e \(\beta_1\)) da equação de regressão que foram estimados.
Vamos usar o conjunto de dados trees, disponível na R, com medidas de 31 cerejeiras: Girth (diâmetro do tronco, em polegadas), Height (altura, em pés) e Volume (volume de madeira, em pés cúbicos). Começamos com a regressão simples de Volume sobre Girth.
data(trees)
# Estimação "na mão", via fórmulas escalares
x <- trees$Girth
y <- trees$Volume
xbar <- mean(x)
ybar <- mean(y)
Sxy <- sum((x - xbar) * (y - ybar))
Sxx <- sum((x - xbar)^2)
beta1_hat <- Sxy / Sxx
beta0_hat <- ybar - beta1_hat * xbar
c(beta0_hat = beta0_hat, beta1_hat = beta1_hat) beta0_hat beta1_hat
-36.943459 5.065856 Comparando com lm():
rls <- lm(Volume ~ Girth, data = trees)
coef(rls)(Intercept) Girth
-36.943459 5.065856 Os valores coincidem, como esperado.
Interpretação de \(\hat\beta_1\): para cada polegada adicional de diâmetro do tronco, o volume esperado de madeira aumenta em aproximadamente 5.07 pés cúbicos.
plot(Volume ~ Girth, data = trees, pch = 19, col = "steelblue",
xlab = "Diâmetro do tronco", ylab = "Volume")
abline(rls, col = "firebrick", lwd = 2)
Da minimização de \(SQ(\beta_0,\beta_1)\) decorrem as equações normais
\[ \sum_{i=1}^n \hat e_i = 0 \qquad \text{e} \qquad \sum_{i=1}^n x_i \hat e_i = 0 \; , \]
que valem exatamente pela forma como \(\hat\beta_0\) e \(\hat\beta_1\) foram construídos. Usando essas duas igualdades é possível mostrar que a variabilidade total de \(y\) em torno de sua média se decompõe em uma parte explicada pelo modelo e uma parte residual:
\[ \underbrace{\sum_{i=1}^n (y_i - \bar y)^2}_{SQTotal} = \underbrace{\sum_{i=1}^n (\hat y_i - \bar y)^2}_{SQReg} + \underbrace{\sum_{i=1}^n (y_i - \hat y_i)^2}_{SQRes} \; . \]
Essa é a identidade fundamental da análise de variância (ANOVA) da regressão:
Esse último é o mesmo \(SQRes\) que aparece em \(\hat\sigma^2 = SQRes/(n-2)\).
Os graus de liberdade seguem a mesma lógica de contagem de parâmetros:
Organiza-se isso na tabela de análise de variância:
| Fonte | SQ | gl | QM | F |
|---|---|---|---|---|
| Regressão | \(SQReg\) | \(1\) | \(QMReg = SQReg/1\) | \(QMReg/QMRes\) |
| Resíduo | \(SQRes\) | \(n-2\) | \(QMRes = SQRes/(n-2)\) | |
| Total | \(SQTotal\) | \(n-1\) |
Sob \(H_0: \beta_1 = 0\), \(F = QMReg/QMRes \sim F_{1,\,n-2}\). Na regressão simples, essa estatística \(F\) é equivalente ao quadrado da estatística \(t\) usada para testar \(\beta_1\) (Seção sobre testes de coeficientes, adiante): \(F = t^2\), ou seja, os dois testes respondem exatamente à mesma pergunta.
# Decomposição manual da soma de quadrados
yhat <- fitted(rls)
SQTotal <- sum((y - mean(y))^2)
SQReg <- sum((yhat - mean(y))^2)
SQRes <- sum((y - yhat)^2)
c(SQTotal = SQTotal, SQReg = SQReg, SQRes = SQRes) SQTotal SQReg SQRes
8106.0839 7581.7813 524.3025 gl_reg <- 1
gl_res <- nrow(trees) - 2
QMReg <- SQReg/gl_reg
QMRes <- SQRes/gl_res
Fstat <- QMReg/QMRes
c(QMReg = QMReg, QMRes = QMRes, F = Fstat,
p_valor = pf(Fstat, gl_reg, gl_res, lower.tail = FALSE)) QMReg QMRes F p_valor
7.581781e+03 1.807940e+01 4.193603e+02 8.644334e-19 Conferindo com a tabela produzida pelo R usando a função anova(), e com o quadrado da estatística \(t\) de Girth em summary(rls):
anova(rls)Analysis of Variance Table
Response: Volume
Df Sum Sq Mean Sq F value Pr(>F)
Girth 1 7581.8 7581.8 419.36 < 2.2e-16 ***
Residuals 29 524.3 18.1
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tval <- summary(rls)$coefficients["Girth", "t value"]
c(t2 = tval^2, F_manual = Fstat) t2 F_manual
419.3603 419.3603 A notação escalar fica pouco prática assim que temos mais de uma covariável. A notação matricial generaliza a regressão simples e múltipla com a mesma formulação além de ser base base para algumas extensões vistas no curso (splines, por exemplo, nada mais são representados por regressões lineares em uma base de funções transformada de \(x\) e portanto escritas na mesma notação matricial).
Empilhando as \(n\) observações,
\[ \mathbf{y} = \mathbf{X}\boldsymbol\beta + \mathbf{e} \; , \]
em que:
Para a regressão simples, \(\mathbf{X}\) tem \(p=2\) colunas: uma de 1’s e uma com os valores de \(x_i\).
X <- model.matrix(~ Girth, data = trees)
head(X) (Intercept) Girth
1 1 8.3
2 1 8.6
3 1 8.8
4 1 10.5
5 1 10.7
6 1 10.8O vetor \(\hat{\boldsymbol\beta}\) minimiza
\[ SQ(\boldsymbol\beta) = (\mathbf{y} - \mathbf{X}\boldsymbol\beta)' (\mathbf{y} - \mathbf{X}\boldsymbol\beta) \; , \]
cuja solução (assumindo \(\mathbf{X}'\mathbf{X}\) não singular) é
\[ \hat{\boldsymbol\beta} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{y} \; . \]
y <- trees$Volume
beta_hat <- solve(crossprod(X), crossprod(X, y))
drop(beta_hat)(Intercept) Girth
-36.943459 5.065856 que reproduz o mesmo resultado de lm() acima.
Note que a fórmula escalar \(\hat\beta_1 = S_{xy}/S_{xx}\) é apenas o caso particular de \((\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{y}\) quando \(\mathbf{X}\) tem duas colunas.
Sob as suposições do modelo:
O estimador de \(\sigma^2\) é
\[ \hat\sigma^2 = \frac{SQRes}{n-p} = \frac{(\mathbf{y}-\mathbf{X}\hat{\boldsymbol\beta})'(\mathbf{y}-\mathbf{X}\hat{\boldsymbol\beta})}{n-p} \; , \]
com \(p\) o número de colunas de \(\mathbf{X}\) (parâmetros estimados, incluindo o intercepto) e \(SQRes\) a mesma soma de quadrados residual definida na análise de variância da seção anterior. A forma matricial apenas generaliza os graus de liberdade de \(n-2\) para \(n-p\).
O erro \(e_i\) é uma quantidade teórica, não observável. O resíduo \(\hat e_i = y_i - \hat y_i\) é sua contrapartida observável, calculada a partir do modelo ajustado. Em notação matricial,
\[ \hat{\mathbf{e}} = \mathbf{y} - \mathbf{X}\hat{\boldsymbol\beta} = \mathbf{y} - \mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{y} = \mathbf{y} - \mathbf{H}\mathbf{y} = (\mathbf{I} - \mathbf{H})\mathbf{y} \; , \]
em que \(\mathbf{H} = \mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\) é a matriz de projeção (hat matrix), que projeta \(\mathbf{y}\) no espaço gerado pelas colunas de \(\mathbf{X}\).
Quando o modelo inclui intercepto:
Os resíduos brutos têm variâncias diferentes entre si mesmo sob homocedasticidade dos erros, \(\text{Var}(\hat e_i) = \sigma^2(1-h_{ii})\), por isso é comum trabalhar com resíduos padronizados (dividindo por \(\hat\sigma\sqrt{1-h_{ii}}\)) ao inspecionar diagnósticos.
par(mfrow = c(2,2), mar = c(3,3,1.5,0.5), mgp = c(2,1,0))
plot(rls)
par(mfrow = c(1, 1))Os quatro painéis padrão do R mostram, respectivamente:
A regressão múltipla nada mais é que a mesma formulação matricial \(\mathbf{y} = \mathbf{X}\boldsymbol\beta + \mathbf{e}\) com mais colunas em \(\mathbf{X}\). Toda a álgebra de estimação ( \(\hat{\boldsymbol\beta} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{y}\)) e as propriedades da seção anterior permanecem idênticas.
Fazendo o cálculo passo a passo temos que, uma vez especificada a matriz \(X\), as demais expressões e/ou códigos não se alteram.
X <- model.matrix(~ Girth + Height, data = trees)
y <- trees$Volume
drop(solve(crossprod(X), crossprod(X, y)))(Intercept) Girth Height
-57.9876589 4.7081605 0.3392512 Usando a função lm() do R chega-se no mesmo resultado.
rlm <- lm(Volume ~ Girth + Height, data = trees)
summary(rlm)
Call:
lm(formula = Volume ~ Girth + Height, data = trees)
Residuals:
Min 1Q Median 3Q Max
-6.4065 -2.6493 -0.2876 2.2003 8.4847
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -57.9877 8.6382 -6.713 2.75e-07 ***
Girth 4.7082 0.2643 17.816 < 2e-16 ***
Height 0.3393 0.1302 2.607 0.0145 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.882 on 28 degrees of freedom
Multiple R-squared: 0.948, Adjusted R-squared: 0.9442
F-statistic: 255 on 2 and 28 DF, p-value: < 2.2e-16Na regressão múltipla, cada coeficiente \(\beta_j\) é interpretado como o efeito parcial de \(x_j\) sobre \(y\), mantendo as demais covariáveis constantes (ceteris paribus). Isto é diferente da regressão simples, em que \(\beta_1\) mistura o efeito de \(x\) com o de qualquer variável omitida correlacionada a ele. Para o modelo acima: mantendo Height fixo, cada polegada adicional de Girth aumenta o volume esperado em 4.71 pés cúbicos; mantendo Girth fixo, cada pé adicional de altura aumenta o volume esperado em 0.34 pés cúbicos.
Um cuidado importante: quando as covariáveis são correlacionadas entre si (multicolinearidade), os coeficientes individuais ficam menos precisos (variâncias maiores) e mais difíceis de interpretar isoladamente, mesmo que o modelo como um todo faça boas previsões. O fator de inflação de variância (VIF) é a ferramenta usual para diagnosticar isso.
A mesma decomposição \(SQTotal = SQReg + SQRes\) vale para a regressão múltipla, apenas generalizando os graus de liberdade: com \(p\) colunas em \(\mathbf{X}\) (incluindo o intercepto), \(SQReg\) tem \(p-1\) graus de liberdade (um para cada coeficiente de inclinação) e \(SQRes\) tem \(n-p\).
| Fonte | SQ | gl | QM | F |
|---|---|---|---|---|
| Regressão | \(SQReg\) | \(p-1\) | \(QMReg = SQReg/(p-1)\) | \(QMReg/QMRes\) |
| Resíduo | \(SQRes\) | \(n-p\) | \(QMRes = SQRes/(n-p)\) | |
| Total | \(SQTotal\) | \(n-1\) |
Sob \(H_0: \beta_1 = \dots = \beta_{p-1} = 0\) (todos os coeficientes de inclinação são nulos — o teste F global, retomado na próxima seção), \(F = QMReg/QMRes \sim F_{p-1,\,n-p}\).
Chamar anova() em um único modelo múltiplo, em vez de comparar dois modelos, tem uma leitura diferente da que vimos para modelos aninhados: o R decompõe \(SQReg\) sequencialmente, termo a termo, na ordem em que aparecem na fórmula:
anova(rlm)Analysis of Variance Table
Response: Volume
Df Sum Sq Mean Sq F value Pr(>F)
Girth 1 7581.8 7581.8 503.1503 < 2e-16 ***
Height 1 102.4 102.4 6.7943 0.01449 *
Residuals 28 421.9 15.1
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1A linha Girth mostra a soma de quadrados explicada por Girth sozinho (coincide com o \(SQReg\) do modelo simples), e a linha Height mostra a soma de quadrados adicional explicada por Height, dado que Girth já está no modelo. Quando as covariáveis são correlacionadas como no caso de Girth e Height, essa soma de quadrados sequencial depende da ordem em que os termos entram na fórmula, o mesmo fenômeno de multicolinearidade mencionado anteriormente.
Sob normalidade dos erros, para testar \(H_0: \beta_j = 0\) usa-se a estatística
\[ t_j = \frac{\hat\beta_j}{\text{ep}(\hat\beta_j)} \sim t_{n-p} \text{ sob } H_0 \; , \]
em que \(\text{ep}(\hat\beta_j)\) é a raiz do elemento \(j\) da diagonal de \(\hat\sigma^2(\mathbf{X}'\mathbf{X})^{-1}\). Essa é a coluna Pr(>|t|) na saída de summary().
O teste \(F\) global é exatamente a linha “Regressão” da tabela de análise de variância vista nas seções anteriores: compara o modelo ajustado a um modelo apenas com intercepto (nenhuma covariável explica \(y\)), testando \(H_0: \beta_1 = \dots = \beta_{p-1} = 0\). A partir das mesmas somas de quadrados \(SQReg\), \(SQRes\) e \(SQTotal\) já definidas, o coeficiente de determinação é
\[ R^2 = 1 - \frac{SQRes}{SQTotal} = \frac{SQReg}{SQTotal} \]
mede a proporção da variância de \(y\) explicada pelo modelo. Como \(R^2\) sempre aumenta (ou permanece igual) ao adicionar covariáveis, mesmo que irrelevantes, usa-se o \(R^2\) ajustado, que penaliza pelo número de parâmetros, o que permite comparar modelos com números diferentes de covariáveis:
\[ R^2_{aj} = 1 - \frac{SQRes/(n-p)}{SQTotal/(n-1)} \; . \]
Para os modelos ajustados até aqui temos:
c(R2.rls = summary(rls)$r.squared,
R2aj.rls = summary(rls)$adj.r.squared,
R2.rlm = summary(rlm)$r.squared,
R2aj.rlm = summary(rlm)$adj.r.squared) R2.rls R2aj.rls R2.rlm R2aj.rlm
0.9353199 0.9330895 0.9479500 0.9442322 Os testes \(t\) individuais e o teste \(F\) global são casos particulares de uma pergunta mais geral: será que uma combinação linear dos parâmetros é igual a um valor de interesse? Por exemplo, no modelo múltiplo, faria sentido perguntar se o efeito de Girth é igual ao efeito de Height (\(H_0: \beta_{Girth} - \beta_{Height} = 0\)), ou se a soma de dois efeitos atinge um valor específico.
Uma combinação linear. Seja \(\mathbf{c}\) um vetor \(p \times 1\) de constantes conhecidas. Para testar \(H_0: \mathbf{c}'\boldsymbol\beta = m\), usa-se
\[ t = \frac{\mathbf{c}'\hat{\boldsymbol\beta} - m} {\sqrt{\hat\sigma^2\, \mathbf{c}'(\mathbf{X}'\mathbf{X})^{-1}\mathbf{c}}} \sim t_{n-p} \text{ sob } H_0 \; , \]
que decorre diretamente de \(\text{Var}(\mathbf{c}'\hat{\boldsymbol\beta}) = \sigma^2\, \mathbf{c}'(\mathbf{X}'\mathbf{X})^{-1}\mathbf{c}\). O teste \(t\) de um coeficiente individual é o caso particular em que \(\mathbf{c}\) é um vetor com 1 na posição de \(\beta_j\) e 0 nas demais.
Várias combinações simultaneamente. De forma mais geral, para testar \(q\) combinações lineares ao mesmo tempo, \(H_0: \mathbf{C}\boldsymbol\beta = \mathbf{m}\), com \(\mathbf{C}\) uma matriz \(q \times p\) de posto completo, usa-se
\[ F = \frac{(\mathbf{C}\hat{\boldsymbol\beta}-\mathbf{m})' \left[\mathbf{C}(\mathbf{X}'\mathbf{X})^{-1}\mathbf{C}'\right]^{-1} (\mathbf{C}\hat{\boldsymbol\beta}-\mathbf{m}) / q}{\hat\sigma^2} \sim F_{q,\,n-p} \text{ sob } H_0 \; , \]
que se reduz à estatística \(t^2\) acima quando \(q=1\), e ao teste F global já visto quando \(\mathbf{C}\boldsymbol\beta\) é o vetor de todos os coeficientes de inclinação. Esse teste também é equivalente a comparar, via anova(), o modelo completo a um modelo reduzido em que a restrição \(\mathbf{C}\boldsymbol\beta=\mathbf{m}\) já está imposta, o que segue a mesma lógica de modelos aninhados vista na Seção de Escolha de Modelos. Pacotes como car (função linearHypothesis()) automatizam esse cálculo para hipóteses gerais; aqui fazemos o cálculo manualmente para explicitar de onde vêm os números.
Exemplo. Testar, no modelo múltiplo Volume ~ Girth + Height, se o efeito de Girth é igual ao de Height (\(\mathbf{c} = (0, 1, -1)'\), \(m=0\)):
Xm <- model.matrix(rlm)
XtX_inv <- solve(crossprod(Xm))
beta_hat_m <- coef(rlm)
sigma2_hat_m <- SQRes_m / gl_res_m
c_vec <- c(0, 1, -1) # H0: beta_Girth - beta_Height = 0
est <- sum(c_vec * beta_hat_m)
var_est <- sigma2_hat_m * drop(t(c_vec) %*% XtX_inv %*% c_vec)
t_stat <- est / sqrt(var_est)
p_valor <- 2 * pt(abs(t_stat), df = gl_res_m, lower.tail = FALSE)
c(estimativa = est, t = t_stat, p_valor = p_valor) estimativa t p_valor
4.368909e+00 1.248282e+01 5.841082e-13 Conferindo pela via do modelo restrito (se \(\beta_{Girth}=\beta_{Height}=b\), o modelo se reduz a Volume ~ I(Girth + Height)), e usando de novo a relação \(F = t^2\):
modelo_restrito <- lm(Volume ~ I(Girth + Height), data = trees)
anova(modelo_restrito, rlm)Analysis of Variance Table
Model 1: Volume ~ I(Girth + Height)
Model 2: Volume ~ Girth + Height
Res.Df RSS Df Sum of Sq F Pr(>F)
1 29 2769.92
2 28 421.92 1 2348 155.82 5.841e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1c(t2 = t_stat^2) t2
155.8207 O valor de \(F\) da tabela anova() deve coincidir com t2 acima, já que a restrição testada envolve uma única combinação linear (\(q=1\)).
Quando um modelo é um caso particular de outro (aninhado), o teste \(F\) parcial via anova() compara diretamente as somas de quadrados residuais. Vamos então considerar a comparação três modelos incrementais. Além dos dois (rls e rlm) já vistos vamos incluir com referência inicial o modelo sem covariáveis, ou seja, ajustar simplesmente e média que pode ser definido por
rl0 <- lm(Volume ~ 1, data = trees) # modelo "nulo"A análise a seguir mede a contribuição incremental de cada termo adicionado ao modelo.
anova(rl0, rls, rlm)Analysis of Variance Table
Model 1: Volume ~ 1
Model 2: Volume ~ Girth
Model 3: Volume ~ Girth + Height
Res.Df RSS Df Sum of Sq F Pr(>F)
1 30 8106.1
2 29 524.3 1 7581.8 503.1503 < 2e-16 ***
3 28 421.9 1 102.4 6.7943 0.01449 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Para comparar modelos não necessariamente aninhados, usam-se critérios que penalizam a complexidade do modelo, como AIC e BIC (quanto menor, melhor):
knitr::kable(
data.frame(
modelo = c("Volume ~ Girth", "Volume ~ Girth + Height"),
AIC = c(AIC(rls), AIC(rlm)),
BIC = c(BIC(rls), BIC(rlm))
),
digits = 2,
caption = "Comparação dos modelos por critérios de informação."
)| modelo | AIC | BIC |
|---|---|---|
| Volume ~ Girth | 181.64 | 185.95 |
| Volume ~ Girth + Height | 176.91 | 182.65 |
Procedimentos automáticos (stepwise, step()) podem ser úteis para explorar um grande número de covariáveis candidatas, mas têm limitações importantes:
Uma alternativa a escolher variáveis “dentro ou fora” do modelo é penalizar a magnitude dos coeficientes durante a própria estimação, em vez de minimizar apenas \(SQRes(\boldsymbol\beta)\). Isso é especialmente útil quando há muitas covariáveis candidatas, correlacionadas entre si (multicolinearidade), ou quando \(p\) é grande em relação a \(n\). Essas são situações em que \((\mathbf{X}'\mathbf{X})^{-1}\) existe mas é numericamente instável, e os estimadores de mínimos quadrados, embora não viesados, passam a ter variância muito alta.
Os métodos de regularização mais utilizados são:
Em todos os casos, \(\lambda\) controla o compromisso entre viés e variância (quanto maior \(\lambda\), mais encolhimento, mais viés e menos variância) e é escolhido por validação cruzada, minimizando o erro de previsão estimado. Nunca deve-se usar o ajuste dentro da própria amostra, que sempre favoreceria \(\lambda=0\) (mínimos quadrados sem penalização).
Vamos considerar o exemplo a seguir no qual expandimos a base de covariáveis para ilustrar o efeito da penalização. Adicionamos ao modelo termos quadráticos e de interação que são naturalmente bastante correlacionados com os termos lineares que os geram.
X_ext <- model.matrix(~ Girth + Height + I(Girth^2) + I(Height^2) + Girth:Height,
data = trees)[, -1] # remove a coluna de intercepto;
y <- trees$Volumeif (requireNamespace("glmnet", quietly = TRUE)) {
library(glmnet)
set.seed(2026) # para reprodutibilidade da validação cruzada
ridge.cv <- cv.glmnet(X_ext, y, alpha = 0) # alpha = 0: ridge
lasso.cv <- cv.glmnet(X_ext, y, alpha = 1) # alpha = 1: lasso
elastic.cv <- cv.glmnet(X_ext, y, alpha = 0.5) # alpha = 0.5: elastic net
coef(ridge.cv, s = "lambda.min")
coef(lasso.cv, s = "lambda.min")
coef(elastic.cv, s = "lambda.min")
}6 x 1 sparse Matrix of class "dgCMatrix"
lambda.min
(Intercept) -0.954780386
Girth -2.145261632
Height .
I(Girth^2) 0.197092735
I(Height^2) 0.001083918
Girth:Height 0.016492548Para comparar os quatro ajustes lado a lado — mínimos quadrados (MQ, sem penalização), ridge, lasso e elastic net —, ajustamos por MQ o mesmo conjunto de termos usado em X_ext e organizamos os coeficientes em uma única tabela:
if (requireNamespace("glmnet", quietly = TRUE)) {
# MQ sem penalização, nos mesmos termos usados em X_ext -------------------
rlm_ext <- lm(Volume ~ Girth * Height + I(Girth^2) + I(Height^2), data = trees)
# Extrai um vetor nomeado de coeficientes de um objeto cv.glmnet ----------
extrai_coef <- function(fit, s = "lambda.min") {
cf <- as.matrix(coef(fit, s = s))
setNames(cf[, 1], rownames(cf))
}
termos <- names(coef(rlm_ext)) # nomes em comum entre os quatro ajustes
tabela_coef <- data.frame(
termo = termos,
MQ = coef(rlm_ext)[termos],
ridge = extrai_coef(ridge.cv)[termos],
lasso = extrai_coef(lasso.cv)[termos],
elastic_net = extrai_coef(elastic.cv)[termos],
row.names = NULL
)
knitr::kable(
tabela_coef, digits = 3,
caption = "Coeficientes por mínimos quadrados, ridge, lasso e elastic net (mesmos termos; penalizados com $\\lambda$ escolhido por validação cruzada)."
)
}| termo | MQ | ridge | lasso | elastic_net |
|---|---|---|---|---|
| (Intercept) | 6.607 | -26.357 | -14.176 | -0.955 |
| Girth | -5.122 | 1.308 | 0.000 | -2.145 |
| Height | 0.295 | 0.069 | 0.000 | 0.000 |
| I(Girth^2) | 0.164 | 0.070 | 0.167 | 0.197 |
| I(Height^2) | -0.005 | 0.001 | 0.002 | 0.001 |
| Girth:Height | 0.066 | 0.016 | 0.000 | 0.016 |
Vale comparar os coeficientes obtidos por ridge e lasso, para o \(\lambda\) escolhido por validação cruzada, com os do modelo Volume ~ Girth * Height + I(Girth^2) + I(Height^2) ajustado por mínimos quadrados sem penalização. Espera-se que o lasso zere alguns dos termos quadráticos/interação, cuja contribuição é pequena diante da correlação com os termos lineares.
Alguns padrões esperados nessa comparação:
Os exercícios a seguir incluem implementações computacionais. Nestas, sempre que possível, faça o cálculo “na mão” por operações matriciais (model.matrix(), crossprod(), solve()) e confira o resultado com lm().
Exercício 1. Para o exemplo do texto com os dados de trees obtenha por operações matriciais os resultados reportados por summary() e anova() para o modelo de regressão múltipla.
Exercício 2. Usando o conjunto de dados cars (base R: speed, velocidade, e dist, distância de frenagem):
dist sobre speed, calculando \(\hat{\boldsymbol\beta}\) tanto pelas fórmulas escalares quanto pela forma matricial \((\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{y}\) (com model.matrix(), crossprod(), solve()). Confira que os dois métodos e lm() concordam.anova(). Confirme numericamente que \(F = t^2\).Exercício 3. Ainda com cars, crie a covariável speed2 = speed^2 e ajuste dist ~ speed + speed2.
lm().anova(). Interprete a soma de quadrados sequencial associada a speed2: o termo quadrático melhora significativamente o ajuste?dist ~ speed e dist ~ speed + speed2 por AIC, BIC e teste F (via anova() com os dois modelos aninhados). As conclusões dos três critérios concordam entre si?Exercício 4. Usando o modelo Volume ~ Girth + Height (dados trees):
anova(), usando a relação \(F = t^2\).Exercício 5. Usando o conjunto de dados mtcars (ou outro, à sua escolha, com várias covariáveis correlacionadas):
mpg sobre todas as demais variáveis (mpg ~ .) e observe o número de coeficientes e possíveis sinais de multicolinearidade (erros-padrão elevados, sinais contra-intuitivos, VIF alto).glmnet::cv.glmnet(), escolhendo \(\lambda\) por validação cruzada, e compare os coeficientes obtidos com os da regressão sem penalização. Quais variáveis o lasso zera?Exercício 6. Para cada uma das extensões do curso (splines/GAM, modelos não lineares, heterocedásticos, efeitos aleatórios, dados no intervalo unitário), identifique explicitamente qual suposição do modelo linear clássico revisado neste material está sendo relaxada, e dê um exemplo de situação prática (diferente dos exemplos do material) em que essa suposição seria violada.