15  Modelos de Regressão

Neste capítulo estudamos os principais modelos para problemas em que a variável resposta é contínua (regressão). Partimos do modelo clássico de regressão linear — simples e múltipla —, com a estimação por mínimos quadrados, sua geometria (projeções ortogonais), o teorema de Gauss–Markov, a inferência sob normalidade, a interpretação dos coeficientes, as medidas de qualidade do ajuste e o diagnóstico por resíduos, influência e heterocedasticidade. Em seguida, tratamos de seleção de modelos (\(C_p\), AIC e BIC), regressão polinomial, splines e modelos aditivos, variáveis indicadoras e interações, dos métodos de regularização (Ridge, Lasso e Elastic Net, incluindo sua interpretação bayesiana) e da regressão logística, ponte natural entre regressão e classificação, que generalizamos para os Modelos Lineares Generalizados. Por fim, aprendemos a aplicar Árvores de Decisão, Florestas Aleatórias e o k-NN quando a resposta é contínua, avaliando os modelos por métricas como o erro médio absoluto, o erro médio quadrático e o coeficiente de determinação.

15.1 Diferença entre Classificação e Regressão

Conforme estudado na introdução ao aprendizado supervisionado, a diferença crucial entre a classificação e a regressão reside no tipo da variável de resposta. Formalmente, seja \(\mathbf{X}\) uma matriz de dados de treinamento com \(n\) observações e \(p\) variáveis, e \(\mathbf{y}\) o vetor da variável resposta (ou alvo).

• Na Classificação, a resposta \(\mathbf{y}\) assume valores categóricos ou rótulos discretos (ex: \(y \in \{0, 1\}\)).
• Na Regressão, a resposta \(\mathbf{y}\) assume valores reais quantitativos contínuos (ex: \(y \in \mathbb{R}\) ou \(y \in \mathbb{R}_{+}\)).

Em ambos os casos, supomos que os dados são gerados por um modelo da forma

\[ Y = f(\mathbf{X}) + \varepsilon, \]

em que \(f\) é uma função fixa e desconhecida e \(\varepsilon\) é um erro aleatório com média zero, independente de \(\mathbf{X}\). Em regressão, sob a perda quadrática, a melhor previsão possível de \(Y\) dado \(\mathbf{X} = \mathbf{x}\) é a esperança condicional \(f(\mathbf{x}) = \mathbb{E}[Y \mid \mathbf{X} = \mathbf{x}]\), chamada função de regressão (Hastie, Tibshirani, e Friedman 2009). Todos os modelos deste capítulo são maneiras diferentes de estimar essa função:

• Modelos paramétricos (regressão linear, polinomial, logística) assumem uma forma funcional para \(f\) e estimam um número finito de parâmetros.
• Modelos não paramétricos (árvores, florestas, k-NN) não fixam a forma de \(f\) e deixam que os dados a determinem, ao custo de exigirem mais observações e serem menos interpretáveis.

Embora algoritmos como Árvores de Decisão e Florestas Aleatórias sejam intuitivamente conhecidos para classificação, eles podem ser facilmente adaptados para realizar tarefas de regressão com alta eficiência preditiva (James et al. 2013).

15.1.1 Um Pouco de História

O termo regressão foi cunhado por Francis Galton no século XIX, ao estudar a relação entre a altura de pais e filhos (Galton 1886). Galton observou que filhos de pais muito altos tendiam a ser altos, mas menos que os pais — havia uma “regressão em direção à média” (regression towards mediocrity). O método de mínimos quadrados, que usamos para estimar a reta de regressão, é ainda mais antigo, publicado por Legendre (1805) e Gauss (1809) para problemas de astronomia. No século XX, a regressão linear tornou-se a ferramenta central da Estatística aplicada, e seus refinamentos — regressão logística (Cox 1958), modelos lineares generalizados (Nelder e Wedderburn 1972), Ridge (Hoerl e Kennard 1970) e Lasso (Tibshirani 1996) — formam a base de boa parte do aprendizado de máquina moderno.

---

15.2 Regressão Linear Simples

15.2.1 O Modelo

A regressão linear simples relaciona uma resposta \(Y\) a um único preditor \(X\):

\[ Y_i = \beta_0 + \beta_1 x_i + \varepsilon_i, \qquad i = 1, \dots, n, \]

em que \(\beta_0\) é o intercepto (valor esperado de \(Y\) quando \(x = 0\)), \(\beta_1\) é a inclinação (variação esperada em \(Y\) para um aumento unitário em \(x\)) e \(\varepsilon_i\) é o erro aleatório.

15.2.2 Suposições do Modelo

As suposições clássicas sobre os erros são:

1. Linearidade: \(\mathbb{E}[Y_i \mid x_i] = \beta_0 + \beta_1 x_i\), isto é, \(\mathbb{E}[\varepsilon_i] = 0\).
2. Homocedasticidade: \(\operatorname{Var}(\varepsilon_i) = \sigma^2\), constante para todo \(i\).
3. Independência: os erros \(\varepsilon_1, \dots, \varepsilon_n\) são não correlacionados (ou independentes).
4. Normalidade (necessária apenas para inferência exata — testes \(t\), intervalos de confiança): \(\varepsilon_i \sim N(0, \sigma^2)\).

Para fins puramente preditivos, as suposições 2 a 4 são menos críticas; para inferência sobre os coeficientes, elas são essenciais e devem ser verificadas pela análise de resíduos.

15.2.3 Estimação por Mínimos Quadrados

O método de Mínimos Quadrados Ordinários (MQO; Ordinary Least Squares - OLS) escolhe \(\hat\beta_0\) e \(\hat\beta_1\) que minimizam a soma dos quadrados dos resíduos:

\[ S(\beta_0, \beta_1) = \sum_{i=1}^n (y_i - \beta_0 - \beta_1 x_i)^2. \]

Derivação. Igualamos as derivadas parciais a zero:

\[ \begin{aligned} \frac{\partial S}{\partial \beta_0} &= -2\sum_{i=1}^n (y_i - \beta_0 - \beta_1 x_i) = 0, \\ \frac{\partial S}{\partial \beta_1} &= -2\sum_{i=1}^n x_i(y_i - \beta_0 - \beta_1 x_i) = 0. \end{aligned} \]

Essas são as equações normais. Da primeira, dividindo por \(n\): \(\bar{y} = \beta_0 + \beta_1\bar{x}\), logo \(\hat\beta_0 = \bar{y} - \hat\beta_1\bar{x}\). Substituindo na segunda e usando \(\sum x_i(y_i - \bar y) = \sum (x_i - \bar x)(y_i - \bar y)\) e \(\sum x_i (x_i - \bar x) = \sum (x_i - \bar x)^2\), obtemos

\[ \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{S_{xy}}{S_{xx}}, \qquad \hat\beta_0 = \bar{y} - \hat\beta_1 \bar{x}. \]

Como \(S\) é uma função quadrática convexa (a matriz hessiana é positiva definida sempre que os \(x_i\) não são todos iguais), esse ponto crítico é o mínimo global. Note ainda que \(\hat\beta_1 = r_{xy}\, s_y / s_x\), em que \(r_{xy}\) é o coeficiente de correlação de Pearson: a inclinação é a correlação “reescalada”.

Propriedades importantes da reta ajustada \(\hat y_i = \hat\beta_0 + \hat\beta_1 x_i\) e dos resíduos \(e_i = y_i - \hat y_i\):

• A reta passa pelo ponto médio \((\bar x, \bar y)\).
• \(\sum_i e_i = 0\) e \(\sum_i x_i e_i = 0\) (consequência direta das equações normais).
• Sob as suposições 1 a 3, \(\hat\beta_0\) e \(\hat\beta_1\) são não viesados e, pelo Teorema de Gauss-Markov, possuem a menor variância entre todos os estimadores lineares não viesados (são os BLUE — Best Linear Unbiased Estimators).
• \(\operatorname{Var}(\hat\beta_1) = \sigma^2 / S_{xx}\), e \(\sigma^2\) é estimado sem viés por \(\hat\sigma^2 = \frac{1}{n-2}\sum_i e_i^2\).

Exemplo numérico resolvido: reta de mínimos quadrados

Considere \(n = 5\) observações:

\(i\)	\(x_i\)	\(y_i\)	\(x_i - \bar x\)	\(y_i - \bar y\)	\((x_i-\bar x)(y_i - \bar y)\)	\((x_i - \bar x)^2\)
1	1	2	\(-2\)	\(-2\)	4	4
2	2	4	\(-1\)	0	0	1
3	3	5	0	1	0	0
4	4	4	1	0	0	1
5	5	5	2	1	2	4
Soma	15	20			6	10

Médias: \(\bar x = 3\) e \(\bar y = 4\). Então

\[ \hat\beta_1 = \frac{6}{10} = 0{,}6, \qquad \hat\beta_0 = 4 - 0{,}6 \times 3 = 2{,}2, \qquad \hat y = 2{,}2 + 0{,}6\,x. \]

Valores ajustados e resíduos:

\(x_i\)	\(y_i\)	\(\hat y_i\)	\(e_i\)	\(e_i^2\)	\((y_i - \bar y)^2\)
1	2	2,8	\(-0{,}8\)	0,64	4
2	4	3,4	0,6	0,36	0
3	5	4,0	1,0	1,00	1
4	4	4,6	\(-0{,}6\)	0,36	0
5	5	5,2	\(-0{,}2\)	0,04	1
Soma			0	2,40	6

Verifique que \(\sum e_i = 0\). Temos \(SQ_{Res} = 2{,}4\), \(SQ_{Tot} = 6\), \(\hat\sigma^2 = 2{,}4/3 = 0{,}8\) e \(\operatorname{EP}(\hat\beta_1) = \sqrt{0{,}8/10} \approx 0{,}283\). A estatística \(t = 0{,}6/0{,}283 \approx 2{,}12\) (com 3 graus de liberdade) não é suficiente para rejeitar \(H_0: \beta_1 = 0\) ao nível de 5% — com tão poucos dados, há muita incerteza. Usaremos esses números novamente nas seções de métricas e de árvores.

---

15.3 Regressão Linear Múltipla

15.3.1 Formulação Matricial

Com \(p\) preditores, o modelo é

\[ Y_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \dots + \beta_p x_{ip} + \varepsilon_i, \]

ou, em notação matricial, \(\mathbf{y} = \mathbf{X}\boldsymbol\beta + \boldsymbol\varepsilon\), em que

\[ \mathbf{y} = \begin{pmatrix} y_1 \\ \vdots \\ y_n \end{pmatrix}, \quad \mathbf{X} = \begin{pmatrix} 1 & x_{11} & \cdots & x_{1p} \\ \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n1} & \cdots & x_{np} \end{pmatrix}, \quad \boldsymbol\beta = \begin{pmatrix} \beta_0 \\ \vdots \\ \beta_p \end{pmatrix}, \]

\(\mathbf{X}\) é a matriz de delineamento (\(n \times (p+1)\), com uma coluna de uns para o intercepto) e \(\mathbb{E}[\boldsymbol\varepsilon] = \mathbf{0}\), \(\operatorname{Var}(\boldsymbol\varepsilon) = \sigma^2\mathbf{I}_n\).

15.3.2 O Estimador de Mínimos Quadrados

Queremos minimizar \(S(\boldsymbol\beta) = (\mathbf{y} - \mathbf{X}\boldsymbol\beta)^\top(\mathbf{y} - \mathbf{X}\boldsymbol\beta) = \mathbf{y}^\top\mathbf{y} - 2\boldsymbol\beta^\top\mathbf{X}^\top\mathbf{y} + \boldsymbol\beta^\top\mathbf{X}^\top\mathbf{X}\boldsymbol\beta\). Usando as regras de derivação matricial \(\partial(\mathbf{a}^\top\boldsymbol\beta)/\partial\boldsymbol\beta = \mathbf{a}\) e \(\partial(\boldsymbol\beta^\top\mathbf{A}\boldsymbol\beta)/\partial\boldsymbol\beta = 2\mathbf{A}\boldsymbol\beta\) (para \(\mathbf{A}\) simétrica):

\[ \frac{\partial S}{\partial \boldsymbol\beta} = -2\mathbf{X}^\top\mathbf{y} + 2\mathbf{X}^\top\mathbf{X}\boldsymbol\beta = \mathbf{0} \quad\Longrightarrow\quad \mathbf{X}^\top\mathbf{X}\,\hat{\boldsymbol\beta} = \mathbf{X}^\top\mathbf{y}. \]

Se \(\mathbf{X}\) tem posto coluna completo (colunas linearmente independentes), \(\mathbf{X}^\top\mathbf{X}\) é invertível e

\[ \boxed{\hat{\boldsymbol\beta} = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{y}.} \]

Os valores ajustados são \(\hat{\mathbf{y}} = \mathbf{X}\hat{\boldsymbol\beta} = \mathbf{H}\mathbf{y}\), com \(\mathbf{H} = \mathbf{X}(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\), a matriz chapéu (hat matrix). Geometricamente, \(\hat{\mathbf{y}}\) é a projeção ortogonal de \(\mathbf{y}\) sobre o espaço gerado pelas colunas de \(\mathbf{X}\), e o vetor de resíduos \(\mathbf{e} = (\mathbf{I} - \mathbf{H})\mathbf{y}\) é ortogonal a esse espaço (\(\mathbf{X}^\top\mathbf{e} = \mathbf{0}\)).

Sob as suposições do modelo:

\[ \mathbb{E}[\hat{\boldsymbol\beta}] = \boldsymbol\beta, \qquad \operatorname{Var}(\hat{\boldsymbol\beta}) = \sigma^2(\mathbf{X}^\top\mathbf{X})^{-1}, \qquad \hat\sigma^2 = \frac{\mathbf{e}^\top\mathbf{e}}{n - p - 1}. \]

Os elementos diagonais \(h_{ii}\) de \(\mathbf{H}\) são as alavancas (leverage) e indicam o quanto cada observação “puxa” o próprio ajuste — úteis no diagnóstico de pontos influentes.

Verificação matricial do exemplo simples

Para os dados do exemplo anterior, \(\mathbf{X}\) tem uma coluna de uns e a coluna \(x = (1,2,3,4,5)\):

\[ \mathbf{X}^\top\mathbf{X} = \begin{pmatrix} 5 & 15 \\ 15 & 55 \end{pmatrix}, \qquad \mathbf{X}^\top\mathbf{y} = \begin{pmatrix} 20 \\ 66 \end{pmatrix}, \qquad (\mathbf{X}^\top\mathbf{X})^{-1} = \frac{1}{50}\begin{pmatrix} 55 & -15 \\ -15 & 5 \end{pmatrix}. \]

Logo \(\hat\beta_0 = (55\cdot 20 - 15\cdot 66)/50 = 110/50 = 2{,}2\) e \(\hat\beta_1 = (-15\cdot 20 + 5\cdot 66)/50 = 30/50 = 0{,}6\), idênticos ao cálculo pela fórmula escalar.

15.3.3 Geometria dos Mínimos Quadrados: a Matriz Chapéu

A interpretação geométrica do MQO é uma das ideias mais férteis da Estatística. Considere \(\mathbf{y}\) como um vetor de \(\mathbb{R}^n\) e seja \(\mathcal{C}(\mathbf{X}) = \{\mathbf{X}\mathbf{b} : \mathbf{b} \in \mathbb{R}^{p+1}\}\) o espaço coluna de \(\mathbf{X}\), um subespaço de dimensão \(p+1\). Minimizar \(\|\mathbf{y} - \mathbf{X}\boldsymbol\beta\|^2\) é procurar o ponto de \(\mathcal{C}(\mathbf{X})\) mais próximo de \(\mathbf{y}\), e o teorema da projeção da Álgebra Linear garante que esse ponto é único e caracterizado pela ortogonalidade do erro: \(\mathbf{X}^\top(\mathbf{y} - \mathbf{X}\hat{\boldsymbol\beta}) = \mathbf{0}\) — que são exatamente as equações normais. O nome “normais” vem justamente daí: o resíduo é normal (perpendicular) ao espaço coluna.

A matriz \(\mathbf{H} = \mathbf{X}(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\) é o operador de projeção ortogonal sobre \(\mathcal{C}(\mathbf{X})\) e possui as seguintes propriedades:

1. Simetria: \(\mathbf{H}^\top = \mathbf{H}\) (imediato, pois \((\mathbf{X}^\top\mathbf{X})^{-1}\) é simétrica).
2. Idempotência: \(\mathbf{H}^2 = \mathbf{X}(\mathbf{X}^\top\mathbf{X})^{-1}\underbrace{\mathbf{X}^\top\mathbf{X}(\mathbf{X}^\top\mathbf{X})^{-1}}_{\mathbf{I}}\mathbf{X}^\top = \mathbf{H}\). Projetar duas vezes é o mesmo que projetar uma vez.
3. \(\mathbf{H}\mathbf{X} = \mathbf{X}\) e \((\mathbf{I} - \mathbf{H})\mathbf{X} = \mathbf{0}\): vetores que já estão em \(\mathcal{C}(\mathbf{X})\) não se alteram; \(\mathbf{I} - \mathbf{H}\) projeta no complemento ortogonal e também é simétrica e idempotente.
4. Autovalores: como \(\mathbf{H}\mathbf{v} = \lambda\mathbf{v}\) implica \(\lambda\mathbf{v} = \mathbf{H}^2\mathbf{v} = \lambda^2\mathbf{v}\), todo autovalor satisfaz \(\lambda \in \{0, 1\}\).
5. Traço: usando \(\operatorname{tr}(\mathbf{A}\mathbf{B}) = \operatorname{tr}(\mathbf{B}\mathbf{A})\), \[\operatorname{tr}(\mathbf{H}) = \operatorname{tr}\big((\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{X}\big) = \operatorname{tr}(\mathbf{I}_{p+1}) = p + 1.\] Para uma matriz idempotente, traço = posto = número de autovalores iguais a 1. Por isso \(\operatorname{tr}(\mathbf{H})\) é interpretado como o número de graus de liberdade do ajuste, ideia que generalizaremos para Ridge e splines.
6. Alavancas: como \(\mathbf{H}\) é idempotente e simétrica, \(h_{ii} = \sum_j h_{ij}^2 \geq h_{ii}^2\), logo \(0 \leq h_{ii} \leq 1\); com intercepto vale ainda \(h_{ii} \geq 1/n\). A alavanca média é \(\bar h = (p+1)/n\), e observações com \(h_{ii} > 2(p+1)/n\) costumam ser sinalizadas como de alta alavanca. Na regressão simples, \(h_{ii} = \dfrac{1}{n} + \dfrac{(x_i - \bar x)^2}{S_{xx}}\): a alavanca cresce com a distância de \(x_i\) ao centro dos dados.

Como \(\mathbf{e} = (\mathbf{I} - \mathbf{H})\mathbf{y}\) e \(\hat{\mathbf{y}} = \mathbf{H}\mathbf{y}\) são ortogonais, vale o teorema de Pitágoras \(\|\mathbf{y}\|^2 = \|\hat{\mathbf{y}}\|^2 + \|\mathbf{e}\|^2\); aplicado a \(\mathbf{y} - \bar y\mathbf{1}\), ele fornece a decomposição \(SQ_{Tot} = SQ_{Reg} + SQ_{Res}\) apresentada adiante.

15.3.4 Propriedades do Estimador e o Teorema de Gauss–Markov

Escrevendo \(\hat{\boldsymbol\beta} = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top(\mathbf{X}\boldsymbol\beta + \boldsymbol\varepsilon) = \boldsymbol\beta + (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\boldsymbol\varepsilon\), obtemos imediatamente \(\mathbb{E}[\hat{\boldsymbol\beta}] = \boldsymbol\beta\) e

\[ \operatorname{Var}(\hat{\boldsymbol\beta}) = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\,(\sigma^2\mathbf{I})\,\mathbf{X}(\mathbf{X}^\top\mathbf{X})^{-1} = \sigma^2(\mathbf{X}^\top\mathbf{X})^{-1}. \]

O resultado central da teoria clássica afirma que nenhum outro estimador linear e não viesado faz melhor.

Teorema (Gauss–Markov). Sob o modelo \(\mathbf{y} = \mathbf{X}\boldsymbol\beta + \boldsymbol\varepsilon\), com \(\mathbb{E}[\boldsymbol\varepsilon] = \mathbf{0}\), \(\operatorname{Var}(\boldsymbol\varepsilon) = \sigma^2\mathbf{I}\) e \(\mathbf{X}\) de posto completo, seja \(\tilde{\boldsymbol\beta} = \mathbf{C}\mathbf{y}\) qualquer estimador linear não viesado de \(\boldsymbol\beta\). Então \(\operatorname{Var}(\tilde{\boldsymbol\beta}) - \operatorname{Var}(\hat{\boldsymbol\beta})\) é positiva semidefinida. Em particular, para todo vetor \(\mathbf{a}\), \(\operatorname{Var}(\mathbf{a}^\top\tilde{\boldsymbol\beta}) \geq \operatorname{Var}(\mathbf{a}^\top\hat{\boldsymbol\beta})\): o MQO é o BLUE.

Prova. Como \(\mathbb{E}[\mathbf{C}\mathbf{y}] = \mathbf{C}\mathbf{X}\boldsymbol\beta\) deve ser igual a \(\boldsymbol\beta\) para todo \(\boldsymbol\beta\), precisamos de \(\mathbf{C}\mathbf{X} = \mathbf{I}\). Escreva \(\mathbf{C} = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top + \mathbf{D}\); então \(\mathbf{C}\mathbf{X} = \mathbf{I} + \mathbf{D}\mathbf{X}\), o que força \(\mathbf{D}\mathbf{X} = \mathbf{0}\). Assim,

\[ \operatorname{Var}(\tilde{\boldsymbol\beta}) = \sigma^2\mathbf{C}\mathbf{C}^\top = \sigma^2\Big[(\mathbf{X}^\top\mathbf{X})^{-1} + (\mathbf{X}^\top\mathbf{X})^{-1}\underbrace{\mathbf{X}^\top\mathbf{D}^\top}_{\mathbf{0}} + \underbrace{\mathbf{D}\mathbf{X}}_{\mathbf{0}}(\mathbf{X}^\top\mathbf{X})^{-1} + \mathbf{D}\mathbf{D}^\top\Big] = \operatorname{Var}(\hat{\boldsymbol\beta}) + \sigma^2\mathbf{D}\mathbf{D}^\top. \]

Como \(\mathbf{a}^\top\mathbf{D}\mathbf{D}^\top\mathbf{a} = \|\mathbf{D}^\top\mathbf{a}\|^2 \geq 0\), a diferença é positiva semidefinida, com igualdade apenas se \(\mathbf{D} = \mathbf{0}\). \(\blacksquare\)

Três observações sobre o alcance do teorema:

• Ele não exige normalidade, apenas média zero, homocedasticidade e ausência de correlação.
• Ele compara o MQO apenas com estimadores lineares e não viesados. Estimadores viesados, como o Ridge, podem ter erro quadrático médio menor — é exatamente isso que a regularização explora.
• Se a homocedasticidade falha, o MQO continua não viesado, mas deixa de ser o BLUE; o melhor estimador passa a ser o de mínimos quadrados ponderados (ver adiante).

15.3.5 Estimação de \(\sigma^2\): por que dividir por \(n - p - 1\)?

Como \((\mathbf{I} - \mathbf{H})\mathbf{X} = \mathbf{0}\), o vetor de resíduos depende apenas dos erros: \(\mathbf{e} = (\mathbf{I} - \mathbf{H})\mathbf{y} = (\mathbf{I} - \mathbf{H})\boldsymbol\varepsilon\). Usando a idempotência e o “truque do traço” (\(\mathbb{E}[\boldsymbol\varepsilon^\top\mathbf{A}\boldsymbol\varepsilon] = \operatorname{tr}(\mathbf{A}\operatorname{Var}(\boldsymbol\varepsilon))\) quando \(\mathbb{E}[\boldsymbol\varepsilon] = \mathbf{0}\)):

\[ \mathbb{E}[\mathbf{e}^\top\mathbf{e}] = \mathbb{E}\big[\boldsymbol\varepsilon^\top(\mathbf{I} - \mathbf{H})\boldsymbol\varepsilon\big] = \sigma^2\operatorname{tr}(\mathbf{I}_n - \mathbf{H}) = \sigma^2(n - p - 1). \]

Portanto \(s^2 = \hat\sigma^2 = \mathbf{e}^\top\mathbf{e}/(n - p - 1)\) é não viesado para \(\sigma^2\). Intuitivamente, perdemos \(p+1\) graus de liberdade porque os resíduos são forçados a satisfazer as \(p+1\) restrições \(\mathbf{X}^\top\mathbf{e} = \mathbf{0}\). Note também que \(\operatorname{Var}(e_i) = \sigma^2(1 - h_{ii})\): os resíduos não têm variância constante, mesmo quando os erros têm — observações de alta alavanca têm resíduos artificialmente pequenos, o que motiva a padronização dos resíduos.

15.3.6 Interpretação dos Coeficientes

Na regressão múltipla, \(\beta_j\) é a variação esperada em \(Y\) para um aumento de uma unidade em \(x_j\), mantendo todos os demais preditores constantes (ceteris paribus). Isso difere da inclinação de uma regressão simples de \(Y\) em \(x_j\), que mistura o efeito de \(x_j\) com o de variáveis correlacionadas a ela. Por exemplo, ao regredir o consumo de combustível (mpg) de automóveis no peso (wt) e na potência (hp), o coeficiente de wt mede o efeito do peso entre carros de mesma potência.

Cuidados na interpretação:

• Associação não implica causalidade: em dados observacionais, \(\beta_j\) pode refletir variáveis de confusão omitidas.
• O intercepto pode não ter interpretação prática quando \(x = \mathbf{0}\) está fora do domínio dos dados.
• Coeficientes dependem das unidades; para comparar a importância relativa, padronize os preditores.

15.3.7 Inferência

Distribuição amostral sob normalidade. Se, além das suposições anteriores, \(\boldsymbol\varepsilon \sim N_n(\mathbf{0}, \sigma^2\mathbf{I})\), valem três fatos fundamentais:

1. \(\hat{\boldsymbol\beta} = \boldsymbol\beta + (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\boldsymbol\varepsilon \sim N_{p+1}\big(\boldsymbol\beta,\; \sigma^2(\mathbf{X}^\top\mathbf{X})^{-1}\big)\), por ser transformação linear de um vetor normal.
2. \((n-p-1)\,s^2/\sigma^2 = \boldsymbol\varepsilon^\top(\mathbf{I} - \mathbf{H})\boldsymbol\varepsilon/\sigma^2 \sim \chi^2_{n-p-1}\), pois uma forma quadrática de um vetor normal padrão em uma matriz simétrica idempotente de posto \(r\) tem distribuição qui-quadrado com \(r\) graus de liberdade (teorema de Cochran), e o posto de \(\mathbf{I} - \mathbf{H}\) é igual ao seu traço, \(n - p - 1\).
3. \(\hat{\boldsymbol\beta}\) e \(s^2\) são independentes: \(\operatorname{Cov}\big(\hat{\boldsymbol\beta}, \mathbf{e}\big) = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top(\sigma^2\mathbf{I})(\mathbf{I} - \mathbf{H}) = \mathbf{0}\), e, para vetores conjuntamente normais, covariância nula implica independência.

Combinando os três fatos com a definição da distribuição \(t\) de Student (normal padrão dividida pela raiz de uma qui-quadrado independente sobre seus graus de liberdade), para cada coeficiente:

\[ t_j = \frac{\hat\beta_j - \beta_j}{\operatorname{EP}(\hat\beta_j)} \sim t_{n-p-1}, \qquad \operatorname{EP}(\hat\beta_j) = s\sqrt{[(\mathbf{X}^\top\mathbf{X})^{-1}]_{jj}}. \]

Daí seguem o teste \(t\) de \(H_0: \beta_j = 0\) (estatística \(\hat\beta_j/\operatorname{EP}(\hat\beta_j)\)) e o intervalo de confiança de nível \(1 - \alpha\):

\[ \hat\beta_j \pm t_{n-p-1;\, 1-\alpha/2}\;\operatorname{EP}(\hat\beta_j). \]

Teste \(F\) para modelos encaixados. Mais geralmente, queremos testar se um modelo reduzido (R), obtido impondo \(q\) restrições lineares — por exemplo, \(\beta_{j_1} = \dots = \beta_{j_q} = 0\) —, é suficiente em comparação com o modelo completo (C). A estatística é

\[ F = \frac{(SQ_{Res}^{R} - SQ_{Res}^{C})/q}{SQ_{Res}^{C}/(n-p-1)} \sim F_{q,\, n-p-1} \quad \text{sob } H_0. \]

Esboço da derivação. Como \(\mathcal{C}(\mathbf{X}_R) \subset \mathcal{C}(\mathbf{X}_C)\), a diferença \(\mathbf{H}_C - \mathbf{H}_R\) também é uma projeção ortogonal, de posto \(q\), e \(SQ_{Res}^R - SQ_{Res}^C = \|(\mathbf{H}_C - \mathbf{H}_R)\mathbf{y}\|^2\). Sob \(H_0\), \(\mathbf{X}\boldsymbol\beta \in \mathcal{C}(\mathbf{X}_R)\) é anulado por \(\mathbf{H}_C - \mathbf{H}_R\), de modo que esse termo é \(\sigma^2\chi^2_q\); ele é independente de \(SQ_{Res}^C = \|(\mathbf{I} - \mathbf{H}_C)\mathbf{y}\|^2\) porque as duas projeções são ortogonais entre si. A razão de duas qui-quadrado independentes divididas pelos graus de liberdade é uma \(F\).

O teste \(F\) global, de \(H_0: \beta_1 = \dots = \beta_p = 0\), é o caso particular em que o modelo reduzido contém apenas o intercepto (\(SQ_{Res}^R = SQ_{Tot}\), \(q = p\)):

\[ F = \frac{(SQ_{Tot} - SQ_{Res})/p}{SQ_{Res}/(n-p-1)} \sim F_{p,\, n-p-1}. \]

Quando \(q = 1\), vale \(F = t_j^2\): o teste \(F\) para um único coeficiente é equivalente ao teste \(t\) bilateral. No R, a comparação de modelos encaixados é feita com anova(modelo_reduzido, modelo_completo).

15.3.8 Intervalos de Confiança e de Predição

Dado um novo vetor de preditores \(\mathbf{x}_0\) (incluindo o 1 do intercepto), há duas perguntas diferentes:

• Qual é a resposta média \(\mathbb{E}[Y \mid \mathbf{x}_0] = \mathbf{x}_0^\top\boldsymbol\beta\)? A incerteza vem apenas da estimação de \(\boldsymbol\beta\): \(\operatorname{Var}(\mathbf{x}_0^\top\hat{\boldsymbol\beta}) = \sigma^2\,\mathbf{x}_0^\top(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{x}_0\).
• Qual será o valor de uma nova observação \(Y_0 = \mathbf{x}_0^\top\boldsymbol\beta + \varepsilon_0\)? Além da incerteza de estimação, há o erro \(\varepsilon_0\), independente da amostra: \(\operatorname{Var}(Y_0 - \mathbf{x}_0^\top\hat{\boldsymbol\beta}) = \sigma^2\big[1 + \mathbf{x}_0^\top(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{x}_0\big]\).

Daí os intervalos de nível \(1 - \alpha\):

\[ \underbrace{\hat y_0 \pm t_{n-p-1;\,1-\alpha/2}\; s\sqrt{\mathbf{x}_0^\top(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{x}_0}}_{\text{IC para a média}} \qquad\text{e}\qquad \underbrace{\hat y_0 \pm t_{n-p-1;\,1-\alpha/2}\; s\sqrt{1 + \mathbf{x}_0^\top(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{x}_0}}_{\text{intervalo de predição}}. \]

Na regressão simples, \(\mathbf{x}_0^\top(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{x}_0 = \frac{1}{n} + \frac{(x_0 - \bar x)^2}{S_{xx}}\). Os dois intervalos são mais estreitos em \(x_0 = \bar x\) e se alargam à medida que nos afastamos do centro dos dados. Quando \(n \to \infty\), o IC para a média encolhe até um ponto, mas o intervalo de predição converge para \(\mathbf{x}_0^\top\boldsymbol\beta \pm z_{1-\alpha/2}\,\sigma\): nenhuma quantidade de dados elimina o erro irredutível.

Exemplo resolvido: testes e intervalos na regressão simples

Retomemos os dados \(x = (1,2,3,4,5)\), \(y = (2,4,5,4,5)\), com \(\hat y = 2{,}2 + 0{,}6x\), \(s^2 = 0{,}8\), \(S_{xx} = 10\) e \(n - 2 = 3\) graus de liberdade. O quantil é \(t_{3;\,0{,}975} \approx 3{,}182\).

Coeficientes. \(\operatorname{EP}(\hat\beta_1) = \sqrt{0{,}8/10} \approx 0{,}283\) e \(\operatorname{EP}(\hat\beta_0) = \sqrt{0{,}8\,(1/5 + 9/10)} = \sqrt{0{,}88} \approx 0{,}938\). O IC de 95% para \(\beta_1\) é

\[ 0{,}6 \pm 3{,}182 \times 0{,}283 = [-0{,}300;\; 1{,}500], \]

que contém o zero, em acordo com o teste \(t = 2{,}12\) (valor-\(p \approx 0{,}124\)).

Teste \(F\). \(SQ_{Reg} = 6 - 2{,}4 = 3{,}6\) e \(F = (3{,}6/1)/(2{,}4/3) = 4{,}5 = 2{,}121^2 = t^2\), com o mesmo valor-\(p\) de 0,124.

Intervalos em \(x_0 = 4\). \(\hat y_0 = 4{,}6\) e \(\frac{1}{5} + \frac{(4-3)^2}{10} = 0{,}3\).

• IC para a média: \(4{,}6 \pm 3{,}182\sqrt{0{,}8 \times 0{,}3} = 4{,}6 \pm 3{,}182 \times 0{,}490 = [3{,}04;\; 6{,}16]\).
• Intervalo de predição: \(4{,}6 \pm 3{,}182\sqrt{0{,}8 \times 1{,}3} = 4{,}6 \pm 3{,}182 \times 1{,}020 = [1{,}35;\; 7{,}85]\).

Extrapolação para \(x_0 = 6\). Agora \(\frac{1}{5} + \frac{9}{10} = 1{,}1\), \(\hat y_0 = 5{,}8\), e o IC para a média é \([2{,}81;\; 8{,}79]\), quase o dobro da largura do IC em \(x_0 = 4\) — além de depender da hipótese, não verificável, de que a reta continua válida fora do intervalo observado.

15.3.9 Coeficiente de Determinação \(R^2\) e \(R^2\) Ajustado

Quando o modelo tem intercepto, vale a decomposição da soma de quadrados:

\[ \underbrace{\sum_{i}(y_i - \bar y)^2}_{SQ_{Tot}} = \underbrace{\sum_i (\hat y_i - \bar y)^2}_{SQ_{Reg}} + \underbrace{\sum_i (y_i - \hat y_i)^2}_{SQ_{Res}}. \]

(O termo cruzado se anula porque \(\mathbf{e}\) é ortogonal às colunas de \(\mathbf{X}\).) Daí

\[ R^2 = \frac{SQ_{Reg}}{SQ_{Tot}} = 1 - \frac{SQ_{Res}}{SQ_{Tot}}, \]

a proporção da variabilidade de \(Y\) explicada pelo modelo. Na regressão simples, \(R^2 = r_{xy}^2\).

O \(R^2\) nunca diminui quando se adiciona um preditor, mesmo que irrelevante. Para comparar modelos com números diferentes de variáveis, usa-se o \(R^2\) ajustado, que penaliza a complexidade:

\[ R^2_{aj} = 1 - \frac{SQ_{Res}/(n-p-1)}{SQ_{Tot}/(n-1)} = 1 - (1 - R^2)\,\frac{n-1}{n-p-1}. \]

No exemplo numérico: \(R^2 = 1 - 2{,}4/6 = 0{,}60\) e \(R^2_{aj} = 1 - 0{,}4 \times 4/3 \approx 0{,}467\). Critérios de informação como AIC e BIC (discutidos na seção de seleção de modelos, adiante) e, principalmente, a validação cruzada são alternativas mais robustas para seleção de modelos.

15.3.10 Análise de Resíduos

Os resíduos \(e_i\) são a principal ferramenta para verificar as suposições. Os gráficos clássicos (produzidos por plot(modelo) no R) são:

Gráfico	O que verificar	Sinal de problema
Resíduos vs. valores ajustados	Linearidade e homocedasticidade	Padrão curvo (não linearidade); formato de funil (heterocedasticidade)
Q-Q normal dos resíduos	Normalidade	Pontos afastados da reta, caudas pesadas
Escala-locação (\(\sqrt{\lvert r_i \rvert}\) vs. ajustados)	Homocedasticidade	Tendência crescente
Resíduos vs. alavanca (com distância de Cook)	Pontos influentes	Pontos com alta alavanca e grande resíduo

Costuma-se usar resíduos padronizados \(r_i = e_i / (\hat\sigma\sqrt{1 - h_{ii}})\), que têm variância aproximadamente 1; valores \(|r_i| > 3\) merecem atenção. Quando há heterocedasticidade, transformações da resposta (como \(\log y\)) frequentemente ajudam. Para dados ordenados no tempo, a independência pode ser verificada com o gráfico de resíduos vs. ordem ou o teste de Durbin-Watson.

15.3.11 Observações Influentes: Alavanca e Distância de Cook

Três conceitos relacionados, mas distintos, descrevem observações “atípicas”:

• Um outlier na resposta tem resíduo grande: \(y_i\) está longe do valor previsto pelo modelo.
• Um ponto de alta alavanca tem \(\mathbf{x}_i\) atípico (alto \(h_{ii}\)), isto é, está longe do centro da nuvem de preditores.
• Um ponto influente é aquele cuja remoção altera substancialmente as estimativas. Influência exige, em geral, alavanca e resíduo consideráveis ao mesmo tempo.

Como \(\operatorname{Var}(e_i) = \sigma^2(1 - h_{ii})\), comparamos os resíduos na escala adequada. Além dos resíduos padronizados \(r_i = e_i/(s\sqrt{1 - h_{ii}})\), usam-se os resíduos estudentizados externamente (studentized deleted residuals), \(t_i = e_i/(s_{(i)}\sqrt{1 - h_{ii}})\), em que \(s_{(i)}\) é a estimativa de \(\sigma\) obtida sem a observação \(i\); sob o modelo, \(t_i \sim t_{n-p-2}\), e vale a identidade \(t_i = r_i\sqrt{(n-p-2)/(n-p-1-r_i^2)}\), que dispensa reajustar o modelo \(n\) vezes.

A distância de Cook (Cook 1977) mede o deslocamento do vetor de coeficientes quando a observação \(i\) é removida, na métrica natural de \(\operatorname{Var}(\hat{\boldsymbol\beta})\):

\[ D_i = \frac{(\hat{\boldsymbol\beta} - \hat{\boldsymbol\beta}_{(i)})^\top\mathbf{X}^\top\mathbf{X}\,(\hat{\boldsymbol\beta} - \hat{\boldsymbol\beta}_{(i)})}{(p+1)\,s^2} = \frac{r_i^2}{p+1}\cdot\frac{h_{ii}}{1 - h_{ii}}. \]

A segunda igualdade (que decorre da fórmula de atualização de Sherman–Morrison para \((\mathbf{X}^\top\mathbf{X} - \mathbf{x}_i\mathbf{x}_i^\top)^{-1}\)) mostra que \(D_i\) é o produto de um termo de resíduo (\(r_i^2\)) por um termo de alavanca (\(h_{ii}/(1-h_{ii})\)). Regras práticas: \(D_i > 1\) é preocupante; \(D_i > 4/n\) merece inspeção.

Exemplo resolvido: alavanca e influência

Nos dados \(x = (1,2,3,4,5)\), \(y = (2,4,5,4,5)\), as alavancas são \(h_{ii} = \frac{1}{5} + \frac{(x_i - 3)^2}{10}\):

\(i\)	\(x_i\)	\(e_i\)	\(h_{ii}\)	\(r_i\)	\(t_i\)	\(D_i\)
1	1	\(-0{,}8\)	0,6	\(-1{,}414\)	\(-2{,}000\)	1,500
2	2	0,6	0,3	0,802	0,739	0,138
3	3	1,0	0,2	1,250	1,474	0,195
4	4	\(-0{,}6\)	0,3	\(-0{,}802\)	\(-0{,}739\)	0,138
5	5	\(-0{,}2\)	0,6	\(-0{,}354\)	\(-0{,}295\)	0,094

Verifique que \(\sum h_{ii} = 2 = p + 1\). Por exemplo, \(r_1 = -0{,}8/\sqrt{0{,}8 \times 0{,}4} = -1{,}414\) e \(D_1 = \frac{2}{2}\cdot\frac{0{,}6}{0{,}4} = 1{,}5\).

As observações 1 e 5 têm a mesma alavanca (0,6), mas influências muito diferentes. Removendo a observação 1, a reta passa a ser \(\hat y = 3{,}8 + 0{,}2x\) (a inclinação cai de 0,6 para 0,2); removendo a observação 5, obtemos \(\hat y = 2{,}0 + 0{,}7x\), uma mudança bem menor. Alavanca sem resíduo grande não gera influência.

15.3.12 Heterocedasticidade: Diagnóstico e Remédios

Quando \(\operatorname{Var}(\varepsilon_i) = \sigma_i^2\) varia entre observações (por exemplo, a variabilidade da renda cresce com a escolaridade), o MQO continua não viesado, mas (i) deixa de ser o BLUE e (ii) a fórmula \(s^2(\mathbf{X}^\top\mathbf{X})^{-1}\) para a variância fica errada, invalidando testes \(t\), \(F\) e intervalos.

Teste de Breusch–Pagan (Breusch e Pagan 1979). A ideia é verificar se os resíduos ao quadrado podem ser explicados pelos preditores (ou por outras variáveis \(\mathbf{z}_i\)): ajusta-se a regressão auxiliar de \(e_i^2\) sobre \(\mathbf{z}_i\) e calcula-se \(LM = n R^2_{aux}\), que, sob a hipótese nula de homocedasticidade, tem distribuição aproximada \(\chi^2_q\), com \(q\) o número de variáveis em \(\mathbf{z}_i\). (Essa é a versão “estudentizada” de Koenker, padrão da função bptest do pacote lmtest, robusta à não normalidade.)

Erros padrão robustos (sanduíche). Sem supor homocedasticidade, a variância exata do MQO é

\[ \operatorname{Var}(\hat{\boldsymbol\beta}) = (\mathbf{X}^\top\mathbf{X})^{-1}\,\mathbf{X}^\top\boldsymbol\Omega\mathbf{X}\,(\mathbf{X}^\top\mathbf{X})^{-1}, \qquad \boldsymbol\Omega = \operatorname{diag}(\sigma_1^2, \dots, \sigma_n^2), \]

uma forma de “sanduíche”: o “pão” \((\mathbf{X}^\top\mathbf{X})^{-1}\) envolve o “recheio” \(\mathbf{X}^\top\boldsymbol\Omega\mathbf{X}\). White (White 1980) mostrou que, substituindo \(\sigma_i^2\) por \(e_i^2\) (estimador HC0), obtém-se um estimador consistente, embora cada \(e_i^2\) isoladamente seja um péssimo estimador de \(\sigma_i^2\) — o que importa é a média ponderada \(\frac{1}{n}\sum_i e_i^2\mathbf{x}_i\mathbf{x}_i^\top\). Em amostras pequenas, recomenda-se a correção HC3, que usa \(e_i^2/(1 - h_{ii})^2\).

Mínimos quadrados ponderados (MQP). Se \(\operatorname{Var}(\varepsilon_i) = \sigma^2/w_i\) com pesos \(w_i\) conhecidos, multiplicamos cada equação por \(\sqrt{w_i}\): o modelo transformado \(\sqrt{w_i}\,y_i = \sqrt{w_i}\,\mathbf{x}_i^\top\boldsymbol\beta + \sqrt{w_i}\,\varepsilon_i\) tem erros homocedásticos de variância \(\sigma^2\). Aplicando o MQO a ele (com \(\mathbf{W} = \operatorname{diag}(w_1, \dots, w_n)\)):

\[ \hat{\boldsymbol\beta}_{MQP} = \arg\min_{\boldsymbol\beta}\sum_{i=1}^n w_i\,(y_i - \mathbf{x}_i^\top\boldsymbol\beta)^2 = (\mathbf{X}^\top\mathbf{W}\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{W}\mathbf{y}, \]

e, pelo teorema de Gauss–Markov aplicado ao modelo transformado, o MQP é o BLUE. Exemplos: se \(y_i\) é a média de \(m_i\) observações, \(w_i = m_i\); se o desvio padrão é proporcional a \(x_i\), \(w_i = 1/x_i^2\). Quando os pesos são desconhecidos, podem ser estimados em duas etapas (regressão de \(\log e_i^2\) nos preditores, por exemplo), procedimento chamado de MQ generalizados factíveis. O MQP reaparecerá adiante como o “motor” do algoritmo IRLS da regressão logística.

Qual remédio usar?

Se o interesse é inferência e a forma da heterocedasticidade é desconhecida, o caminho mais simples é manter o MQO e reportar erros padrão robustos (HC3). Se a forma da variância é conhecida ou bem modelada, o MQP é mais eficiente. Se a variância cresce com a média e a resposta é positiva, uma transformação (\(\log y\), \(\sqrt{y}\)) ou um modelo linear generalizado adequado (Poisson, Gama) costuma resolver o problema na raiz.

15.3.13 Multicolinearidade e VIF

Há multicolinearidade quando preditores são fortemente correlacionados entre si. Nesse caso, \(\mathbf{X}^\top\mathbf{X}\) fica próxima da singularidade e as variâncias \(\sigma^2[(\mathbf{X}^\top\mathbf{X})^{-1}]_{jj}\) “explodem”: os coeficientes tornam-se instáveis, podem ter sinais contraintuitivos e erros padrão enormes, embora as previsões possam continuar boas.

O diagnóstico usual é o Fator de Inflação da Variância (Variance Inflation Factor - VIF):

\[ \operatorname{VIF}_j = \frac{1}{1 - R_j^2}, \]

em que \(R_j^2\) é o coeficiente de determinação da regressão de \(x_j\) sobre todos os demais preditores. O nome vem da identidade \(\operatorname{Var}(\hat\beta_j) = \frac{\sigma^2}{(n-1)s_{x_j}^2}\operatorname{VIF}_j\). Regra prática: \(\operatorname{VIF} > 5\) indica multicolinearidade moderada e \(\operatorname{VIF} > 10\), severa. Por exemplo, se \(R_j^2 = 0{,}9\), então \(\operatorname{VIF}_j = 10\): a variância de \(\hat\beta_j\) é 10 vezes maior do que seria com preditores não correlacionados. Remédios: remover variáveis redundantes, combiná-las (por exemplo, por componentes principais) ou usar regularização (Ridge).

15.3.14 Variáveis Categóricas: Variáveis Indicadoras (Dummies)

Um preditor categórico com \(k\) níveis é incluído por meio de \(k-1\) variáveis indicadoras (dummies), escolhendo-se um nível de referência. Por exemplo, para a região de um imóvel com níveis {Norte, Sul, Centro} e referência “Centro”:

\[ Y = \beta_0 + \beta_1\, \text{area} + \beta_2 D_{\text{Norte}} + \beta_3 D_{\text{Sul}} + \varepsilon, \]

em que \(D_{\text{Norte}} = 1\) se o imóvel está no Norte e 0 caso contrário. Então \(\beta_2\) é a diferença esperada de preço entre um imóvel no Norte e outro no Centro de mesma área. Usar \(k\) indicadoras junto com o intercepto gera colinearidade perfeita (as indicadoras somam 1) — a chamada armadilha das dummies (dummy variable trap) —, pois \(\mathbf{X}^\top\mathbf{X}\) deixa de ser invertível. No R, lm cria as dummies automaticamente para variáveis do tipo factor; em Python, usa-se pandas.get_dummies(..., drop_first=True) ou OneHotEncoder(drop="first").

15.3.15 Interações

O modelo aditivo supõe que o efeito de \(x_1\) não depende do valor de \(x_2\). Para relaxar essa hipótese, inclui-se um termo de interação:

\[ Y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_1 x_2 + \varepsilon. \]

Agora, o efeito de aumentar \(x_1\) em uma unidade é \(\beta_1 + \beta_3 x_2\), que varia com \(x_2\). Se \(x_2\) é uma dummy (por exemplo, sexo), a interação permite inclinações diferentes para cada grupo. Pelo princípio hierárquico, ao incluir \(x_1x_2\) devem-se manter também os efeitos principais \(x_1\) e \(x_2\), mesmo que não significativos. No R, a fórmula y ~ x1 * x2 equivale a y ~ x1 + x2 + x1:x2.

---

15.4 Seleção de Modelos: \(C_p\), AIC e BIC

Com \(p\) preditores candidatos existem \(2^p\) subconjuntos possíveis. Como escolher? O erro de treinamento (\(SQ_{Res}\) ou \(R^2\)) não serve, pois sempre melhora com mais variáveis. O problema de fundo é o otimismo do erro de treino.

15.4.1 O Otimismo do Erro de Treinamento

Considere um modelo linear com \(d = p + 1\) parâmetros, corretamente especificado, e imagine novas respostas \(\mathbf{y}^{\text{novo}}\) observadas nos mesmos \(\mathbf{x}_i\). Usando os resultados de traço da seção anterior:

\[ \mathbb{E}[SQ_{Res}] = \sigma^2(n - d), \qquad \mathbb{E}\big[\|\mathbf{y}^{\text{novo}} - \hat{\mathbf{y}}\|^2\big] = \sigma^2 n + \sigma^2\operatorname{tr}(\mathbf{H}) = \sigma^2(n + d). \]

A diferença, \(2d\sigma^2\), é o otimismo: em média, o erro de treino subestima o erro em novos dados por \(2d\sigma^2\), e o subestima mais quanto mais parâmetros o modelo tem. Corrigir esse viés é a ideia comum dos critérios a seguir (Hastie, Tibshirani, e Friedman 2009).

15.4.2 \(C_p\) de Mallows, AIC e BIC

• \(C_p\) de Mallows (Mallows 1973): estimando \(\sigma^2\) pelo \(s^2\) do modelo com todos os preditores, \[C_p = \frac{SQ_{Res}(d)}{s^2_{\text{completo}}} - n + 2d.\] Se o modelo com \(d\) parâmetros não tem viés, \(\mathbb{E}[SQ_{Res}] \approx (n - d)\sigma^2\) e \(C_p \approx d\). Procuram-se modelos com \(C_p\) pequeno e próximo de \(d\). (James et al. usam a versão equivalente \((SQ_{Res} + 2d\hat\sigma^2)/n\), que é exatamente “erro de treino + otimismo estimado”.)
• Critério de Informação de Akaike (Akaike 1974): \[\operatorname{AIC} = -2\,\ell(\hat{\boldsymbol\theta}) + 2d,\] em que \(\ell(\hat{\boldsymbol\theta})\) é a log-verossimilhança maximizada. Akaike mostrou que \(-2\ell(\hat{\boldsymbol\theta})/n\) é um estimador viesado da divergência de Kullback–Leibler esperada entre o modelo ajustado e o verdadeiro mecanismo gerador, com viés aproximadamente \(2d/n\): o AIC é a correção desse viés. Na regressão normal, \(-2\ell(\hat{\boldsymbol\theta}) = n\log(SQ_{Res}/n) + \text{constante}\), e o AIC é essencialmente equivalente ao \(C_p\).
• Critério de Informação Bayesiano (Schwarz 1978): \[\operatorname{BIC} = -2\,\ell(\hat{\boldsymbol\theta}) + d\log n.\] Ele surge de uma aproximação de Laplace da verossimilhança marginal \(p(\mathbf{y} \mid \mathcal{M}) = \int p(\mathbf{y} \mid \boldsymbol\theta, \mathcal{M})\,\pi(\boldsymbol\theta)\,d\boldsymbol\theta\): para \(n\) grande, \(-2\log p(\mathbf{y} \mid \mathcal{M}) \approx \operatorname{BIC}\). Escolher o menor BIC aproxima, portanto, escolher o modelo de maior probabilidade a posteriori sob prioris iguais entre modelos.

Como \(\log n > 2\) para \(n \geq 8\), o BIC penaliza a complexidade mais fortemente que o AIC e tende a escolher modelos menores. As duas penalidades refletem objetivos distintos: o BIC é consistente (se o modelo verdadeiro está entre os candidatos, é escolhido com probabilidade tendendo a 1), enquanto o AIC é eficiente para previsão (minimiza assintoticamente o erro de previsão), mas tende a incluir variáveis em excesso.

Exemplo resolvido: comparando modelos para o consumo dos carros (mtcars)

Ajustando modelos para mpg na base mtcars (\(n = 32\)), com o \(s^2\) do modelo completo de seis preditores (wt, hp, qsec, am, disp, drat) para o \(C_p\), obtemos (AIC e BIC calculados pelas funções AIC() e BIC() do R, que contam \(\sigma\) como parâmetro):

Modelo	\(d\)	\(SQ_{Res}\)	\(R^2_{aj}\)	\(C_p\)	AIC	BIC
wt	2	278,32	0,745	18,36	166,03	170,43
wt + hp	3	195,05	0,815	6,49	156,65	162,52
wt + hp + qsec	4	186,06	0,817	6,99	157,14	164,47
wt + qsec + am	4	169,29	0,834	4,20	154,12	161,45
wt + hp + qsec + am	5	160,07	0,837	4,66	154,33	163,12
completo (6 preditores)	7	150,09	0,835	7,00	156,27	167,99

O \(SQ_{Res}\) diminui sempre, mas \(C_p\), AIC e BIC concordam no modelo wt + qsec + am, que também é o escolhido pela função step do R partindo do modelo nulo (direction = "both"), por eliminação regressiva a partir do modelo completo e por eliminação regressiva com penalidade BIC (k = log(n)). O \(R^2_{aj}\) prefere, por uma diferença mínima, o modelo com cinco parâmetros. Note que o modelo completo tem \(C_p = d = 7\) por construção.

15.4.3 Seleção do Melhor Subconjunto e Seleção Passo a Passo

• Melhor subconjunto (best subset): ajusta todos os \(2^p\) modelos e escolhe, para cada tamanho, o de menor \(SQ_{Res}\); depois compara os tamanhos por \(C_p\), AIC, BIC ou validação cruzada. É inviável para \(p\) maior que algumas dezenas (para \(p = 40\), há mais de \(10^{12}\) modelos).
• Seleção progressiva (forward): começa com o modelo nulo e adiciona, a cada passo, a variável que mais melhora o critério; ajusta no máximo \(1 + p(p+1)/2\) modelos.
• Eliminação regressiva (backward): começa com todos os preditores (exige \(n > p\)) e remove, a cada passo, a menos útil.
• Passo a passo (stepwise, ambas as direções): permite adicionar e remover a cada passo.

Cuidados com a seleção automática

Os métodos gulosos não garantem encontrar o melhor subconjunto. Mais grave: após selecionar variáveis com os próprios dados, os valores-\(p\), erros padrão e intervalos de confiança reportados pelo summary do modelo final são otimistas, porque ignoram a etapa de seleção (o problema da inferência pós-seleção). Para previsão, prefira avaliar todo o procedimento (incluindo a seleção) por validação cruzada; para seleção em alta dimensão, o Lasso é uma alternativa contínua e mais estável.

---

15.5 Regressão Polinomial

Quando a relação entre \(x\) e \(Y\) é curva, podemos incluir potências do preditor:

\[ Y = \beta_0 + \beta_1 x + \beta_2 x^2 + \dots + \beta_d x^d + \varepsilon. \]

Apesar de a relação com \(x\) ser não linear, o modelo continua linear nos parâmetros: basta tomar \(x, x^2, \dots, x^d\) como colunas de \(\mathbf{X}\) e usar mínimos quadrados. Pontos de atenção:

• O grau \(d\) controla a complexidade: graus altos produzem curvas que oscilam fortemente nas extremidades e sobreajustam (dilema viés-variância, Capítulo de Aprendizado de Máquina). Na prática, raramente se usa \(d > 3\) ou \(4\); o grau deve ser escolhido por validação cruzada.
• As potências \(x, x^2, x^3\) são altamente correlacionadas; polinômios ortogonais (poly(x, d) no R) evitam problemas numéricos.
• Extrapolar fora do intervalo observado de \(x\) é especialmente perigoso.
• Alternativas mais flexíveis e estáveis incluem splines e modelos aditivos generalizados (James et al. 2013), apresentados a seguir.

---

15.6 Splines e Modelos Aditivos Generalizados

O defeito dos polinômios globais é que cada coeficiente afeta a curva inteira: um ponto discrepante em uma extremidade altera o ajuste em todo o domínio. A ideia das splines é usar polinômios de grau baixo por partes, unidos de forma suave.

15.6.1 Splines de Regressão

Fixados nós \(\xi_1 < \xi_2 < \dots < \xi_K\), uma spline cúbica é uma função que é um polinômio cúbico em cada intervalo entre nós consecutivos e que tem \(f\), \(f'\) e \(f''\) contínuas nos nós. Uma base conveniente é a base de potências truncadas:

\[ f(x) = \beta_0 + \beta_1 x + \beta_2 x^2 + \beta_3 x^3 + \sum_{k=1}^K \beta_{3+k}\,(x - \xi_k)_+^3, \qquad (u)_+ = \max\{u, 0\}. \]

Cada termo \((x - \xi_k)^3_+\) é zero à esquerda de \(\xi_k\) e tem derivadas primeira e segunda contínuas em \(\xi_k\) (apenas a terceira derivada salta), o que garante a suavidade. São, portanto, \(K + 4\) parâmetros, estimados por mínimos quadrados comuns: uma spline de regressão é apenas um modelo linear com colunas construídas a partir de \(x\), e toda a teoria das seções anteriores (inferência, \(\mathbf{H}\), \(C_p\)) se aplica.

Na prática, prefere-se a base de B-splines (numericamente mais estável, função bs() no R) e as splines naturais (ns()), que impõem linearidade fora dos nós extremos. Essa restrição economiza 4 graus de liberdade e reduz a grande variância que as splines cúbicas têm nas bordas. O número de nós (ou de graus de liberdade) é escolhido por validação cruzada, e os nós costumam ser posicionados nos quantis de \(x\).

15.6.2 Splines de Suavização

Uma alternativa elegante evita escolher nós. Procura-se, entre todas as funções com segunda derivada integrável, aquela que minimiza

\[ \sum_{i=1}^n \big(y_i - f(x_i)\big)^2 + \lambda\int f''(t)^2\,dt. \]

O primeiro termo mede o ajuste; o segundo penaliza a rugosidade (curvatura) de \(f\). Com \(\lambda = 0\), qualquer função que interpole os dados é solução; com \(\lambda \to \infty\), \(f''\) é forçada a zero e obtém-se a reta de MQO. Um resultado notável (devido a Reinsch e Wahba) garante que, para qualquer \(\lambda > 0\), a solução é uma spline cúbica natural com nós em todos os \(x_i\) distintos — o problema de dimensão infinita reduz-se a um problema de dimensão \(n\), encolhido como no Ridge. Os valores ajustados são \(\hat{\mathbf{y}} = \mathbf{S}_\lambda\mathbf{y}\), e, por analogia com \(\operatorname{tr}(\mathbf{H}) = p + 1\), definem-se os graus de liberdade efetivos \(\operatorname{gl}(\lambda) = \operatorname{tr}(\mathbf{S}_\lambda)\), que decrescem de \(n\) a 2 à medida que \(\lambda\) cresce.

15.6.3 Modelos Aditivos Generalizados (GAMs)

Com vários preditores, os modelos aditivos (Hastie e Tibshirani 1986) substituem cada termo linear \(\beta_j x_j\) por uma função suave:

\[ Y = \beta_0 + f_1(x_1) + f_2(x_2) + \dots + f_p(x_p) + \varepsilon. \]

Eles preservam a interpretabilidade do modelo aditivo — cada \(f_j\) pode ser visualizada separadamente, “mantendo os demais preditores constantes” — e capturam não linearidades sem especificá-las de antemão. O ajuste clássico é feito pelo algoritmo de retroajuste (backfitting): repetidamente, cada \(f_j\) é reestimada suavizando os resíduos parciais \(y_i - \beta_0 - \sum_{k \neq j} f_k(x_{ik})\) em função de \(x_j\). Modernamente, o pacote mgcv do R (Wood 2017) representa cada \(f_j\) por splines penalizadas e escolhe as penalidades automaticamente (por validação cruzada generalizada ou máxima verossimilhança restrita). Os GAMs também se estendem a respostas binárias e de contagem (o “generalizados” do nome), combinando-se com os modelos lineares generalizados descritos adiante. Sua principal limitação é não capturar interações automaticamente; elas podem ser incluídas por termos do tipo \(f_{jk}(x_j, x_k)\).

---

15.7 Regularização: Ridge, Lasso e Elastic Net

Quando há muitos preditores (eventualmente \(p > n\)) ou forte multicolinearidade, o estimador de MQO tem variância alta — ou nem existe, se \(\mathbf{X}^\top\mathbf{X}\) for singular. A regularização (ou encolhimento, shrinkage) adiciona uma penalidade sobre o tamanho dos coeficientes, aceitando um pequeno viés em troca de grande redução de variância. Em todos os casos, o intercepto não é penalizado e os preditores devem ser padronizados antes do ajuste, pois a penalidade depende da escala.

15.7.1 Regressão Ridge

Proposta por Hoerl e Kennard (Hoerl e Kennard 1970), a regressão Ridge resolve

\[ \hat{\boldsymbol\beta}^{\text{ridge}} = \arg\min_{\boldsymbol\beta} \left\{ \sum_{i=1}^n \Big(y_i - \beta_0 - \sum_{j=1}^p \beta_j x_{ij}\Big)^2 + \lambda\sum_{j=1}^p \beta_j^2 \right\}, \qquad \lambda \geq 0. \]

Com preditores centrados (de modo que o intercepto possa ser tratado à parte), a solução tem forma fechada:

\[ \hat{\boldsymbol\beta}^{\text{ridge}} = (\mathbf{X}^\top\mathbf{X} + \lambda\mathbf{I})^{-1}\mathbf{X}^\top\mathbf{y}. \]

A soma \(\lambda\mathbf{I}\) torna a matriz sempre invertível para \(\lambda > 0\), mesmo com \(p > n\) ou colinearidade perfeita. Quando \(\lambda = 0\), recupera-se o MQO; quando \(\lambda \to \infty\), todos os coeficientes tendem a zero (mas nunca são exatamente zero). No caso de preditores ortonormais (\(\mathbf{X}^\top\mathbf{X} = \mathbf{I}\)), \(\hat\beta_j^{\text{ridge}} = \hat\beta_j^{\text{MQO}}/(1 + \lambda)\): um encolhimento proporcional.

Derivação. A função objetivo \(\|\mathbf{y} - \mathbf{X}\boldsymbol\beta\|^2 + \lambda\boldsymbol\beta^\top\boldsymbol\beta\) tem gradiente \(-2\mathbf{X}^\top\mathbf{y} + 2\mathbf{X}^\top\mathbf{X}\boldsymbol\beta + 2\lambda\boldsymbol\beta\), que se anula em \((\mathbf{X}^\top\mathbf{X} + \lambda\mathbf{I})\boldsymbol\beta = \mathbf{X}^\top\mathbf{y}\). A hessiana \(2(\mathbf{X}^\top\mathbf{X} + \lambda\mathbf{I})\) tem todos os autovalores \(\geq 2\lambda > 0\), logo a função é estritamente convexa e o ponto crítico é o único mínimo.

15.7.2 Ridge: Viés, Variância e Graus de Liberdade Efetivos

Escreva a decomposição em valores singulares \(\mathbf{X} = \mathbf{U}\mathbf{D}\mathbf{V}^\top\), com \(\mathbf{U}\) (\(n \times p\)) e \(\mathbf{V}\) (\(p \times p\)) de colunas ortonormais e \(\mathbf{D} = \operatorname{diag}(d_1 \geq \dots \geq d_p \geq 0)\). Substituindo na forma fechada:

\[ \mathbf{X}\hat{\boldsymbol\beta}^{\text{ridge}} = \sum_{j=1}^p \mathbf{u}_j\,\frac{d_j^2}{d_j^2 + \lambda}\,\mathbf{u}_j^\top\mathbf{y}, \qquad\text{enquanto}\qquad \mathbf{X}\hat{\boldsymbol\beta}^{\text{MQO}} = \sum_{j=1}^p \mathbf{u}_j\,\mathbf{u}_j^\top\mathbf{y}. \]

O Ridge projeta \(\mathbf{y}\) nas mesmas direções que o MQO (as componentes principais de \(\mathbf{X}\)), mas encolhe cada coordenada pelo fator \(d_j^2/(d_j^2 + \lambda) \in (0, 1)\). Direções de pequena variância nos preditores (pequeno \(d_j\)) — justamente as responsáveis pela instabilidade sob multicolinearidade — são as mais encolhidas. Por analogia com \(\operatorname{tr}(\mathbf{H}) = p\), os graus de liberdade efetivos do Ridge são

\[ \operatorname{gl}(\lambda) = \operatorname{tr}\big[\mathbf{X}(\mathbf{X}^\top\mathbf{X} + \lambda\mathbf{I})^{-1}\mathbf{X}^\top\big] = \sum_{j=1}^p \frac{d_j^2}{d_j^2 + \lambda}, \]

que decresce continuamente de \(p\) (em \(\lambda = 0\)) a 0 (quando \(\lambda \to \infty\)).

Com \(\mathbf{A}_\lambda = (\mathbf{X}^\top\mathbf{X} + \lambda\mathbf{I})^{-1}\), os momentos do estimador são

\[ \mathbb{E}[\hat{\boldsymbol\beta}^{\text{ridge}}] = \mathbf{A}_\lambda\mathbf{X}^\top\mathbf{X}\boldsymbol\beta = \boldsymbol\beta - \lambda\mathbf{A}_\lambda\boldsymbol\beta, \qquad \operatorname{Var}(\hat{\boldsymbol\beta}^{\text{ridge}}) = \sigma^2\mathbf{A}_\lambda\mathbf{X}^\top\mathbf{X}\mathbf{A}_\lambda. \]

O viés \(-\lambda\mathbf{A}_\lambda\boldsymbol\beta\) cresce com \(\lambda\), e a variância decresce. Hoerl e Kennard (Hoerl e Kennard 1970) provaram o teorema de existência: há sempre algum \(\lambda > 0\) para o qual o erro quadrático médio total do Ridge, \(\mathbb{E}\|\hat{\boldsymbol\beta}^{\text{ridge}} - \boldsymbol\beta\|^2\), é estritamente menor que o do MQO. A versão com um único preditor ilustra o porquê.

Exemplo resolvido: o \(\lambda\) ótimo do Ridge com um preditor

Com um preditor centrado, \(\hat\beta^{\text{ridge}} = S_{xy}/(S_{xx} + \lambda)\), e

\[ \operatorname{EQM}(\lambda) = \underbrace{\frac{\sigma^2 S_{xx}}{(S_{xx} + \lambda)^2}}_{\text{variância}} + \underbrace{\frac{\beta^2\lambda^2}{(S_{xx} + \lambda)^2}}_{\text{viés}^2}. \]

Derivando em relação a \(\lambda\) e igualando a zero, obtemos \(\beta^2\lambda(S_{xx} + \lambda) = \sigma^2 S_{xx} + \beta^2\lambda^2\), ou seja,

\[ \lambda^* = \frac{\sigma^2}{\beta^2}. \]

Além disso, \(\operatorname{EQM}'(0) = -2\sigma^2/S_{xx}^2 < 0\): partindo do MQO, um pouco de encolhimento sempre reduz o EQM. O encolhimento ótimo é grande quando o ruído é grande em relação ao sinal.

Para os dados do exemplo de regressão simples (\(S_{xx} = 10\), \(S_{xy} = 6\)) e supondo, para ilustração, que os valores verdadeiros fossem \(\beta = 0{,}6\) e \(\sigma^2 = 0{,}8\):

\(\lambda\)	\(\hat\beta^{\text{ridge}}\)	Variância	Viés\(^2\)	EQM
0 (MQO)	0,600	0,0800	0	0,0800
1	0,545	0,0661	0,0030	0,0691
\(\lambda^* = 2{,}22\)	0,491	0,0536	0,0119	0,0655
5	0,400	0,0356	0,0400	0,0756
10	0,300	0,0200	0,0900	0,1100

O EQM cai cerca de 18% no ótimo. Os graus de liberdade efetivos da inclinação em \(\lambda^*\) são \(10/12{,}22 \approx 0{,}82\). Na prática, \(\beta\) e \(\sigma^2\) são desconhecidos, e \(\lambda\) é escolhido por validação cruzada.

15.7.3 Lasso

O Lasso (Least Absolute Shrinkage and Selection Operator), de Tibshirani (Tibshirani 1996), troca a penalidade quadrática pela norma \(\ell_1\):

\[ \hat{\boldsymbol\beta}^{\text{lasso}} = \arg\min_{\boldsymbol\beta} \left\{ \sum_{i=1}^n \Big(y_i - \beta_0 - \sum_{j=1}^p \beta_j x_{ij}\Big)^2 + \lambda\sum_{j=1}^p |\beta_j| \right\}. \]

Não há forma fechada em geral, mas, no caso ortonormal, a solução é o limiar suave (soft thresholding):

\[ \hat\beta_j^{\text{lasso}} = \operatorname{sign}(\hat\beta_j^{\text{MQO}})\,\Big(|\hat\beta_j^{\text{MQO}}| - \tfrac{\lambda}{2}\Big)_+, \]

em que \((a)_+ = \max\{a, 0\}\). Assim, coeficientes pequenos são zerados exatamente: o Lasso faz seleção de variáveis automaticamente e produz modelos esparsos e mais interpretáveis.

Derivação do limiar suave. No caso ortonormal, \(\|\mathbf{y} - \mathbf{X}\boldsymbol\beta\|^2 = \|\mathbf{y} - \mathbf{X}\hat{\boldsymbol\beta}^{\text{MQO}}\|^2 + \|\hat{\boldsymbol\beta}^{\text{MQO}} - \boldsymbol\beta\|^2\), de modo que o problema se separa em \(p\) problemas unidimensionais da forma (escrevendo \(z = \hat\beta_j^{\text{MQO}}\))

\[ \min_{\beta}\; g(\beta) = (z - \beta)^2 + \lambda|\beta|. \]

A função \(|\beta|\) não é diferenciável em zero, mas é convexa; usamos o subgradiente \(\partial|\beta| = \{\operatorname{sign}(\beta)\}\) se \(\beta \neq 0\) e \(\partial|\beta| = [-1, 1]\) se \(\beta = 0\). A condição de otimalidade é \(0 \in -2(z - \beta) + \lambda\,\partial|\beta|\):

• se \(\beta > 0\): \(\beta = z - \lambda/2\), o que só é consistente quando \(z > \lambda/2\);
• se \(\beta < 0\): \(\beta = z + \lambda/2\), consistente quando \(z < -\lambda/2\);
• se \(\beta = 0\): precisamos de \(2z \in [-\lambda, \lambda]\), isto é, \(|z| \leq \lambda/2\).

As três condições juntas formam exatamente \(\hat\beta = \operatorname{sign}(z)(|z| - \lambda/2)_+\). O “intervalo morto” \(|z| \leq \lambda/2\) é a origem matemática da esparsidade; no Ridge, a derivada de \(\beta^2\) em zero é nula e não há esse intervalo.

Descida coordenada. Sem ortonormalidade, o Lasso é resolvido eficientemente por descida coordenada (Friedman, Hastie, e Tibshirani 2010): com os demais coeficientes fixos, a atualização de \(\beta_j\) é um problema unidimensional cuja solução é um limiar suave. Definindo o resíduo parcial \(\mathbf{r}^{(j)} = \mathbf{y} - \sum_{k \neq j}\mathbf{x}_k\beta_k\),

\[ \beta_j \leftarrow \frac{\operatorname{sign}\big(\mathbf{x}_j^\top\mathbf{r}^{(j)}\big)\big(|\mathbf{x}_j^\top\mathbf{r}^{(j)}| - \lambda/2\big)_+}{\mathbf{x}_j^\top\mathbf{x}_j}, \]

e percorrem-se as coordenadas ciclicamente até a convergência. Como o objetivo é convexo, o método converge para o mínimo global. Calculando as soluções para uma sequência decrescente de valores de \(\lambda\) e usando cada solução como ponto de partida da seguinte (warm starts), obtém-se todo o caminho de regularização a um custo pouco maior que o de um único ajuste — essa é a estratégia do glmnet.

Uma consequência útil da condição de otimalidade: se \(|\mathbf{x}_j^\top\mathbf{y}| \leq \lambda/2\) para todo \(j\) (com dados centrados), então \(\hat{\boldsymbol\beta}^{\text{lasso}} = \mathbf{0}\). Logo, \(\lambda_{\max} = 2\max_j|\mathbf{x}_j^\top\mathbf{y}|\) é o menor valor de \(\lambda\) que zera todos os coeficientes, e o caminho de regularização começa nele.

Exemplo numérico: Ridge vs. Lasso no caso ortonormal

Suponha estimativas de MQO \(\hat{\boldsymbol\beta}^{\text{MQO}} = (3{,}0;\; -1{,}5;\; 0{,}4)\) e \(\lambda = 1\).

	\(\beta_1\)	\(\beta_2\)	\(\beta_3\)
MQO	3,0	\(-1{,}5\)	0,4
Ridge: \(\hat\beta_j/(1+\lambda)\)	1,5	\(-0{,}75\)	0,2
Lasso: \(\operatorname{sign}(\hat\beta_j)(\lvert\hat\beta_j\rvert - 0{,}5)_+\)	2,5	\(-1{,}0\)	0

O Ridge encolhe todos os coeficientes pela mesma proporção; o Lasso subtrai uma quantidade fixa e zera o coeficiente pequeno.

Interpretação geométrica. Ambos os problemas são equivalentes a minimizar a soma de quadrados sujeita a uma restrição de orçamento: \(\sum\beta_j^2 \leq s\) (um disco, no Ridge) ou \(\sum|\beta_j| \leq s\) (um losango, no Lasso). As curvas de nível da soma de quadrados (elipses) tendem a tocar o losango em um vértice, onde algum coeficiente é zero; o disco, por não ter vértices, raramente produz zeros exatos (James et al. 2013; Hastie, Tibshirani, e Friedman 2009).

15.7.4 Interpretação Bayesiana da Regularização

A regularização tem uma leitura bayesiana natural. Suponha \(\mathbf{y} \mid \boldsymbol\beta \sim N(\mathbf{X}\boldsymbol\beta, \sigma^2\mathbf{I})\) e coloque uma distribuição a priori sobre os coeficientes. Pelo teorema de Bayes, a log-densidade a posteriori é

\[ \log p(\boldsymbol\beta \mid \mathbf{y}) = -\frac{1}{2\sigma^2}\|\mathbf{y} - \mathbf{X}\boldsymbol\beta\|^2 + \log\pi(\boldsymbol\beta) + \text{constante}. \]

• Priori gaussiana \(\beta_j \overset{\text{iid}}{\sim} N(0, \tau^2)\): \(\log\pi(\boldsymbol\beta) = -\frac{1}{2\tau^2}\sum_j\beta_j^2 + \text{cte}\). Maximizar a posteriori equivale a minimizar \(\|\mathbf{y} - \mathbf{X}\boldsymbol\beta\|^2 + \frac{\sigma^2}{\tau^2}\sum_j\beta_j^2\): o Ridge com \(\lambda = \sigma^2/\tau^2\). Como a posteriori é gaussiana, o estimador Ridge é simultaneamente a moda e a média a posteriori.
• Priori de Laplace (dupla exponencial) \(\pi(\beta_j) = \frac{1}{2b}e^{-|\beta_j|/b}\): \(\log\pi(\boldsymbol\beta) = -\frac{1}{b}\sum_j|\beta_j| + \text{cte}\), e a moda a posteriori é o Lasso com \(\lambda = 2\sigma^2/b\). A densidade de Laplace tem um “bico” em zero e caudas mais pesadas que a normal: acredita-se a priori que muitos coeficientes são praticamente nulos e alguns são grandes. A média a posteriori, porém, não é esparsa; o tratamento bayesiano completo é conhecido como Lasso bayesiano (Park e Casella 2008).

A leitura bayesiana deixa claro o papel de \(\lambda\): uma priori mais concentrada (menor \(\tau^2\) ou \(b\)) significa mais encolhimento. Ela também explica por que padronizar os preditores é importante: uma mesma priori para todos os \(\beta_j\) só faz sentido se eles estiverem na mesma escala.

15.7.5 Elastic Net

O Elastic Net (Zou e Hastie 2005) combina as duas penalidades. Na parametrização usada pelo pacote glmnet (Friedman, Hastie, e Tibshirani 2010):

\[ \hat{\boldsymbol\beta}^{\text{EN}} = \arg\min_{\boldsymbol\beta} \left\{ \frac{1}{2n}\sum_{i=1}^n \Big(y_i - \beta_0 - \mathbf{x}_i^\top\boldsymbol\beta\Big)^2 + \lambda\left[ \alpha\sum_{j=1}^p|\beta_j| + \frac{1-\alpha}{2}\sum_{j=1}^p\beta_j^2 \right] \right\}, \]

com \(\alpha \in [0, 1]\): \(\alpha = 1\) é o Lasso e \(\alpha = 0\) é o Ridge. O Elastic Net seleciona variáveis como o Lasso, mas lida melhor com grupos de preditores correlacionados (o Lasso tende a escolher arbitrariamente um deles e descartar os outros) e pode selecionar mais de \(n\) variáveis quando \(p > n\).

15.7.6 Escolha de \(\lambda\)

O hiperparâmetro \(\lambda\) é escolhido por validação cruzada: ajusta-se o modelo em uma grade de valores de \(\lambda\) (o “caminho de regularização”) e escolhe-se o que minimiza o erro de validação. Uma regra comum é a do “um erro padrão” (lambda.1se no glmnet): escolher o maior \(\lambda\) cujo erro esteja a até um erro padrão do mínimo, privilegiando modelos mais simples.

Método	Penalidade	Zera coeficientes?	Forma fechada?	Indicado quando
MQO	nenhuma	Não	Sim	\(n \gg p\), sem colinearidade
Ridge	\(\lambda\sum\beta_j^2\)	Não	Sim	Muitos preditores com efeitos pequenos; colinearidade
Lasso	\(\lambda\sum\lvert\beta_j\rvert\)	Sim	Não	Poucos preditores relevantes (esparsidade)
Elastic Net	combinação	Sim	Não	Grupos de preditores correlacionados; \(p > n\)

---

15.8 Regressão Logística: Ponte para a Classificação

E se a resposta for binária, \(Y \in \{0, 1\}\)? Ajustar uma regressão linear diretamente a \(Y\) (o “modelo de probabilidade linear”) pode gerar previsões fora do intervalo \([0,1]\). A regressão logística (Cox 1958) modela a probabilidade \(p(\mathbf{x}) = P(Y = 1 \mid \mathbf{X} = \mathbf{x})\) por meio da função logística (sigmoide):

\[ p(\mathbf{x}) = \frac{1}{1 + e^{-(\beta_0 + \boldsymbol\beta^\top\mathbf{x})}} = \frac{e^{\eta}}{1 + e^{\eta}}, \qquad \eta = \beta_0 + \beta_1x_1 + \dots + \beta_p x_p. \]

A função logística \(\sigma(\eta) = 1/(1+e^{-\eta})\) transforma qualquer número real em um valor entre 0 e 1, é crescente e vale \(0{,}5\) em \(\eta = 0\).

15.8.1 Odds e Log-Odds

As chances (odds) de \(Y = 1\) são \(p/(1-p)\). Invertendo a função logística:

\[ \log\left(\frac{p(\mathbf{x})}{1 - p(\mathbf{x})}\right) = \beta_0 + \beta_1x_1 + \dots + \beta_p x_p. \]

O lado esquerdo é a função logito (log-odds), que é linear nos preditores. Portanto, \(\beta_j\) é a variação no log das chances para um aumento unitário em \(x_j\), e \(e^{\beta_j}\) é a razão de chances (odds ratio): as chances são multiplicadas por \(e^{\beta_j}\) a cada unidade de \(x_j\). A regressão logística é um caso particular dos Modelos Lineares Generalizados (Nelder e Wedderburn 1972), com distribuição de Bernoulli e função de ligação logito.

15.8.2 Estimação por Máxima Verossimilhança

Não há mínimos quadrados aqui: os parâmetros são estimados por máxima verossimilhança. Como \(Y_i \sim \text{Bernoulli}(p_i)\), a log-verossimilhança é

\[ \ell(\boldsymbol\beta) = \sum_{i=1}^n \Big[ y_i \log p_i + (1 - y_i)\log(1 - p_i) \Big] = \sum_{i=1}^n \Big[ y_i\,\eta_i - \log\big(1 + e^{\eta_i}\big) \Big]. \]

Seu negativo é a entropia cruzada binária (binary cross-entropy, ou log-loss), a mesma função de perda usada em redes neurais para classificação. Derivando, o vetor gradiente é \(\mathbf{X}^\top(\mathbf{y} - \mathbf{p})\), que não admite solução fechada; o máximo é obtido numericamente pelo método de Newton-Raphson, que, neste contexto, é conhecido como mínimos quadrados reponderados iterativamente (IRLS). A função \(\ell\) é côncava, de modo que o máximo é global. As penalidades Ridge, Lasso e Elastic Net também se aplicam à regressão logística (basta somá-las a \(-\ell\)).

15.8.3 O Algoritmo de Newton–Raphson e o IRLS

Vamos derivar o algoritmo. Como \(\partial p_i/\partial\eta_i = p_i(1 - p_i)\) e \(\partial\eta_i/\partial\boldsymbol\beta = \mathbf{x}_i\), derivando \(\ell(\boldsymbol\beta) = \sum_i[y_i\eta_i - \log(1 + e^{\eta_i})]\):

\[ \nabla\ell(\boldsymbol\beta) = \sum_{i=1}^n (y_i - p_i)\,\mathbf{x}_i = \mathbf{X}^\top(\mathbf{y} - \mathbf{p}), \qquad \nabla^2\ell(\boldsymbol\beta) = -\sum_{i=1}^n p_i(1 - p_i)\,\mathbf{x}_i\mathbf{x}_i^\top = -\mathbf{X}^\top\mathbf{W}\mathbf{X}, \]

com \(\mathbf{W} = \operatorname{diag}\big(p_i(1-p_i)\big)\). Como \(\mathbf{W}\) tem diagonal positiva, a hessiana é negativa definida (se \(\mathbf{X}\) tem posto completo): \(\ell\) é estritamente côncava, e o máximo, se existir, é único. O método de Newton–Raphson aproxima \(\ell\) por sua expansão de Taylor de segunda ordem em torno da estimativa corrente e salta para o máximo dessa parábola:

\[ \boldsymbol\beta^{(t+1)} = \boldsymbol\beta^{(t)} + (\mathbf{X}^\top\mathbf{W}^{(t)}\mathbf{X})^{-1}\mathbf{X}^\top(\mathbf{y} - \mathbf{p}^{(t)}). \]

Reescrevendo, \(\boldsymbol\beta^{(t+1)} = (\mathbf{X}^\top\mathbf{W}\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{W}\mathbf{z}\), com a resposta de trabalho \(\mathbf{z} = \mathbf{X}\boldsymbol\beta^{(t)} + \mathbf{W}^{-1}(\mathbf{y} - \mathbf{p})\). Cada passo de Newton é, portanto, uma regressão por mínimos quadrados ponderados de \(\mathbf{z}\) em \(\mathbf{X}\) com pesos \(p_i(1 - p_i)\), e, como os pesos mudam a cada passo, o método recebe o nome de mínimos quadrados reponderados iterativamente (IRLS). Na convergência, \((\mathbf{X}^\top\hat{\mathbf{W}}\mathbf{X})^{-1}\) — a inversa da informação de Fisher observada — estima a matriz de covariância de \(\hat{\boldsymbol\beta}\), de onde vêm os erros padrão reportados pelo glm.

Exemplo resolvido: uma iteração do IRLS à mão

Considere \(x = (1, 2, 3, 4, 5, 6)\) e \(y = (0, 0, 1, 0, 1, 1)\). Partimos de \(\boldsymbol\beta^{(0)} = (0, 0)\).

Passo 1. Com \(\boldsymbol\beta^{(0)} = \mathbf{0}\), todos os \(p_i = 0{,}5\) e todos os pesos são \(p_i(1-p_i) = 0{,}25\). O gradiente é

\[ \mathbf{X}^\top(\mathbf{y} - \mathbf{p}) = \begin{pmatrix} \sum(y_i - 0{,}5) \\ \sum x_i(y_i - 0{,}5) \end{pmatrix} = \begin{pmatrix} 3 - 3 \\ 14 - 10{,}5 \end{pmatrix} = \begin{pmatrix} 0 \\ 3{,}5 \end{pmatrix}, \]

e a informação é

\[ \mathbf{X}^\top\mathbf{W}\mathbf{X} = 0{,}25\begin{pmatrix} 6 & 21 \\ 21 & 91 \end{pmatrix} = \begin{pmatrix} 1{,}5 & 5{,}25 \\ 5{,}25 & 22{,}75 \end{pmatrix}, \qquad (\mathbf{X}^\top\mathbf{W}\mathbf{X})^{-1} \approx \begin{pmatrix} 3{,}467 & -0{,}800 \\ -0{,}800 & 0{,}2286 \end{pmatrix}. \]

Passo 2. Atualização: \(\boldsymbol\beta^{(1)} = \mathbf{0} + (\mathbf{X}^\top\mathbf{W}\mathbf{X})^{-1}(0;\; 3{,}5)^\top = (-2{,}800;\; 0{,}800)\).

Iterações seguintes (calculadas numericamente):

\(t\)	\(\hat\beta_0^{(t)}\)	\(\hat\beta_1^{(t)}\)	\(\ell(\boldsymbol\beta^{(t)})\)
0	0	0	\(-4{,}1589\)
1	\(-2{,}8000\)	0,8000	\(-2{,}6065\)
2	\(-3{,}8842\)	1,1098	\(-2{,}4849\)
3	\(-4{,}2214\)	1,2061	\(-2{,}4780\)
4	\(-4{,}2489\)	1,2140	\(-2{,}4780\)
5	\(-4{,}2491\)	1,2140	\(-2{,}4780\)

A convergência é quadrática: o número de casas decimais corretas aproximadamente dobra a cada passo. O resultado coincide com glm(y ~ x, family = binomial), que reporta \(\hat\beta_1 = 1{,}214\) com erro padrão 0,913 e razão de chances \(e^{1{,}214} \approx 3{,}37\).

Separação completa

Se existe um hiperplano que separa perfeitamente os \(y_i = 0\) dos \(y_i = 1\) (por exemplo, \(y = (0,0,0,1,1,1)\) com \(x = (1, \dots, 6)\)), a log-verossimilhança cresce indefinidamente quando \(|\boldsymbol\beta| \to \infty\) e o EMV não existe. O IRLS diverge e o software emite avisos como “fitted probabilities numerically 0 or 1 occurred”. Remédios: penalização (Ridge/Lasso), correção de Firth ou reconsiderar o modelo.

15.8.4 Deviance e Teste da Razão de Verossimilhanças

Nos modelos estimados por máxima verossimilhança, o papel de \(SQ_{Res}\) é desempenhado pela deviance:

\[ D = 2\big[\ell(\text{modelo saturado}) - \ell(\hat{\boldsymbol\beta})\big], \]

em que o modelo saturado tem um parâmetro por observação e ajusta os dados perfeitamente. Para respostas binárias, \(\ell(\text{saturado}) = 0\) e \(D = -2\ell(\hat{\boldsymbol\beta})\); para o modelo normal linear, a deviance é exatamente \(SQ_{Res}/\sigma^2\). A deviance nula \(D_0\) corresponde ao modelo só com intercepto.

Para comparar modelos encaixados com \(q\) parâmetros de diferença, usa-se o teste da razão de verossimilhanças (TRV):

\[ \Lambda = D_{R} - D_{C} = 2\big[\ell(\hat{\boldsymbol\beta}_C) - \ell(\hat{\boldsymbol\beta}_R)\big] \;\overset{a}{\sim}\; \chi^2_q \quad \text{sob } H_0, \]

resultado assintótico conhecido como teorema de Wilks. Ele é o análogo do teste \(F\) e, em amostras pequenas, costuma ser mais confiável que o teste de Wald \(z = \hat\beta_j/\operatorname{EP}(\hat\beta_j)\), que depende de uma aproximação quadrática da verossimilhança. Uma medida descritiva análoga ao \(R^2\) é o pseudo-\(R^2\) de McFadden, \(1 - \ell(\hat{\boldsymbol\beta})/\ell_0\).

Exemplo resolvido: deviance e TRV

No exemplo anterior, a proporção amostral de uns é \(1/2\), então \(\ell_0 = 6\log 0{,}5 = -4{,}159\) e \(D_0 = 8{,}318\). Para o modelo ajustado, \(D = -2(-2{,}478) = 4{,}956\). O TRV para \(H_0: \beta_1 = 0\) é

\[ \Lambda = 8{,}318 - 4{,}956 = 3{,}362, \qquad \text{valor-}p = P(\chi^2_1 > 3{,}362) \approx 0{,}067, \]

enquanto o teste de Wald dá \(z = 1{,}214/0{,}913 = 1{,}33\) (valor-\(p \approx 0{,}18\)). Com \(n = 6\), os dois testes discordam bastante — um lembrete de que ambos são aproximações assintóticas. O pseudo-\(R^2\) de McFadden é \(1 - 2{,}478/4{,}159 \approx 0{,}40\), e o AIC é \(D + 2 \times 2 = 8{,}956\).

15.8.5 Da Probabilidade à Classe

Para classificar, escolhe-se um limiar \(c\) (usualmente \(0{,}5\)): \(\hat y = 1\) se \(\hat p(\mathbf{x}) \geq c\). Como \(\hat p \geq 0{,}5 \iff \hat\eta \geq 0\), a fronteira de decisão \(\beta_0 + \boldsymbol\beta^\top\mathbf{x} = 0\) é um hiperplano — assim como no SVM linear (Capítulo de SVM). O limiar pode ser ajustado conforme o custo dos erros e avaliado com as métricas de classificação do Capítulo de Aprendizado de Máquina (matriz de confusão, precisão, sensibilidade, curva ROC).

Exemplo numérico: interpretação da regressão logística

Suponha que um modelo para a probabilidade de inadimplência em função do número de parcelas em atraso \(x\) tenha \(\hat\beta_0 = -3\) e \(\hat\beta_1 = 0{,}5\).

• Para \(x = 4\): \(\hat\eta = -3 + 2 = -1\) e \(\hat p = 1/(1 + e^{1}) \approx 0{,}269\).
• Para \(x = 6\): \(\hat\eta = 0\) e \(\hat p = 0{,}5\) (ponto da fronteira de decisão com \(c = 0{,}5\)).
• Razão de chances: \(e^{0{,}5} \approx 1{,}65\). Cada parcela adicional em atraso multiplica as chances de inadimplência por 1,65 (aumento de 65%), qualquer que seja o valor de \(x\).

Note que o efeito sobre a probabilidade não é constante: de \(x=4\) para \(x=5\), \(\hat p\) passa de 0,269 para \(1/(1+e^{0{,}5}) \approx 0{,}378\); de \(x=10\) para \(x=11\), passa de \(0{,}881\) para \(0{,}924\).

---

15.9 Modelos Lineares Generalizados

A regressão linear normal e a regressão logística são casos particulares de uma mesma família, os Modelos Lineares Generalizados (MLGs) de Nelder e Wedderburn (Nelder e Wedderburn 1972; McCullagh e Nelder 1989). Um MLG tem três componentes:

1. Componente aleatório: \(Y_i\) tem distribuição na família exponencial, \[f(y; \theta, \phi) = \exp\left\{\frac{y\,\theta - b(\theta)}{a(\phi)} + c(y, \phi)\right\},\] com parâmetro canônico \(\theta\) e de dispersão \(\phi\). Derivando a identidade \(\int f = 1\), obtém-se \(\mathbb{E}[Y] = \mu = b'(\theta)\) e \(\operatorname{Var}(Y) = a(\phi)\,b''(\theta)\): a variância é uma função da média, \(V(\mu)\).
2. Componente sistemático: o preditor linear \(\eta_i = \mathbf{x}_i^\top\boldsymbol\beta\).
3. Função de ligação: uma função monótona \(g\) com \(g(\mu_i) = \eta_i\). A ligação que faz \(\theta = \eta\) é chamada canônica.

Distribuição	Suporte de \(Y\)	\(V(\mu)\)	Ligação canônica	Uso típico
Normal	\(\mathbb{R}\)	1	identidade, \(\mu\)	Respostas contínuas
Bernoulli/Binomial	\(\{0,1\}\) ou proporções	\(\mu(1-\mu)\)	logito, \(\log\frac{\mu}{1-\mu}\)	Respostas binárias
Poisson	\(\{0, 1, 2, \dots\}\)	\(\mu\)	log, \(\log\mu\)	Contagens
Gama	\((0, \infty)\)	\(\mu^2\)	inversa, \(1/\mu\) (na prática, usa-se log)	Valores positivos assimétricos (custos, tempos)

Todos os MLGs são ajustados por máxima verossimilhança via o mesmo algoritmo IRLS (escore de Fisher) da seção anterior, com pesos \(w_i = 1/\big[V(\mu_i)\,g'(\mu_i)^2\big]\), e comparados por deviance, TRV e AIC.

15.9.1 Regressão de Poisson

Para contagens — número de sinistros, de visitas a um site, de defeitos —, a regressão linear é inadequada (prevê valores negativos e supõe variância constante). A regressão de Poisson supõe \(Y_i \sim \text{Poisson}(\mu_i)\) com

\[ \log\mu_i = \beta_0 + \beta_1 x_{i1} + \dots + \beta_p x_{ip}. \]

A log-verossimilhança é \(\ell(\boldsymbol\beta) = \sum_i[y_i\eta_i - e^{\eta_i} - \log y_i!]\), e o escore \(\mathbf{X}^\top(\mathbf{y} - \boldsymbol\mu)\) tem exatamente a mesma forma do caso logístico — uma propriedade de todas as ligações canônicas. A interpretação é multiplicativa: \(e^{\beta_j}\) é a razão de taxas (rate ratio), o fator pelo qual a contagem esperada é multiplicada a cada unidade de \(x_j\).

Dois cuidados práticos:

• Exposição (offset): se as contagens foram observadas em tempos ou populações \(t_i\) diferentes, modela-se a taxa: \(\log\mu_i = \log t_i + \mathbf{x}_i^\top\boldsymbol\beta\), em que \(\log t_i\) entra com coeficiente fixo igual a 1.
• Superdispersão: a Poisson impõe \(\operatorname{Var}(Y) = \mathbb{E}[Y]\). Se a estatística \(\hat\phi = \frac{1}{n-p-1}\sum_i (y_i - \hat\mu_i)^2/\hat\mu_i\) é bem maior que 1, os erros padrão estão subestimados; usam-se então o modelo quase-Poisson (que multiplica as variâncias por \(\hat\phi\)) ou a regressão binomial negativa.

Exemplo: rupturas de fios na tecelagem (warpbreaks)

A base warpbreaks do R registra o número de rupturas de fios (breaks) em teares, segundo o tipo de lã (A ou B) e a tensão (L, M ou H). O ajuste glm(breaks ~ wool + tension, family = poisson) fornece

Termo	\(\hat\beta\)	EP	\(e^{\hat\beta}\)
Intercepto	3,692	0,045	40,12
lã B	\(-0{,}206\)	0,052	0,814
tensão M	\(-0{,}321\)	0,060	0,725
tensão H	\(-0{,}518\)	0,064	0,595

Interpretação: para a lã A sob tensão baixa, esperam-se cerca de 40,1 rupturas; a lã B reduz a taxa de rupturas em cerca de 19% (\(1 - 0{,}814\)), e a tensão alta, em cerca de 40%, mantidos os demais fatores. A deviance residual é 210,4 com 50 graus de liberdade, e \(\hat\phi \approx 4{,}26\): há forte superdispersão, de modo que os erros padrão acima são otimistas — com family = quasipoisson, eles seriam multiplicados por \(\sqrt{4{,}26} \approx 2{,}06\).

---

15.10 Árvores de Regressão

O funcionamento de uma Árvore de Regressão é análogo ao de classificação (Capítulo de Árvores de Decisão), porém com modificações essenciais em seus critérios matemáticos de divisão e de previsão de nós. A árvore particiona o espaço dos preditores em \(M\) regiões retangulares disjuntas \(R_1, \dots, R_M\) e ajusta uma constante em cada uma:

\[ \hat f(\mathbf{x}) = \sum_{m=1}^M c_m\, \mathbf{1}\{\mathbf{x} \in R_m\}. \]

15.10.1 Critério de Divisão

Em vez de usar a Entropia ou o Índice de Gini (que medem a impureza de classes), a árvore de regressão escolhe a variável e o ponto de corte que minimizam a Soma dos Quadrados dos Resíduos (RSS - Residual Sum of Squares) dos nós filhos resultantes:

\[RSS = \sum_{i \in R_1} (y_i - \hat{y}_{R_1})^2 + \sum_{i \in R_2} (y_i - \hat{y}_{R_2})^2\] Onde \(R_1\) e \(R_2\) são as duas regiões formadas após a divisão do nó, e \(\hat{y}_{R_1}\) é a média dos valores reais das amostras que caíram na região 1.

Mais precisamente, para cada variável \(j\) e ponto de corte \(s\), definem-se \(R_1(j, s) = \{\mathbf{x}: x_j \leq s\}\) e \(R_2(j,s) = \{\mathbf{x}: x_j > s\}\), e procura-se o par \((j, s)\) que minimiza o RSS. Como é inviável examinar todas as partições possíveis, a árvore é construída de forma gulosa (greedy) e recursiva — a chamada partição binária recursiva do algoritmo CART (Breiman et al. 1984): escolhe-se a melhor divisão no nó atual e o procedimento se repete em cada nó filho até que um critério de parada seja atingido (profundidade máxima, número mínimo de observações por folha etc.). Note que minimizar o RSS equivale a maximizar a redução da variância dentro dos nós.

Exemplo numérico: escolhendo o primeiro corte

Usando os dados \(x = (1, 2, 3, 4, 5)\) e \(y = (2, 4, 5, 4, 5)\) do exemplo de regressão linear, os cortes candidatos são os pontos médios entre valores consecutivos de \(x\):

Corte	\(R_1\) (\(y\))	\(\bar y_{R_1}\)	RSS\(_1\)	\(R_2\) (\(y\))	\(\bar y_{R_2}\)	RSS\(_2\)	RSS total
\(x \leq 1{,}5\)	{2}	2,00	0	{4, 5, 4, 5}	4,50	1,00	1,00
\(x \leq 2{,}5\)	{2, 4}	3,00	2,00	{5, 4, 5}	4,67	0,67	2,67
\(x \leq 3{,}5\)	{2, 4, 5}	3,67	4,67	{4, 5}	4,50	0,50	5,17
\(x \leq 4{,}5\)	{2, 4, 5, 4}	3,75	4,75	{5}	5,00	0	4,75

O melhor corte é \(x \leq 1{,}5\), com RSS total de 1,00 — bem menor que o \(SQ_{Tot} = 6\) do nó raiz e que o \(SQ_{Res} = 2{,}4\) da reta de mínimos quadrados. A árvore com um único corte prevê \(\hat y = 2\) para \(x \leq 1{,}5\) e \(\hat y = 4{,}5\) caso contrário. (Com apenas 5 pontos, isso, evidentemente, não significa que a árvore generalize melhor.)

15.10.2 Previsão Final do Nó

Quando um novo dado de teste percorre a árvore e atinge um nó folha terminal \(t\), a predição do valor numérico predito \(\hat{y}(t)\) é dada simplesmente pela média aritmética dos valores numéricos reais das amostras de treinamento que compõem aquele nó folha \(N(t)\):

\[\hat{y}(t) = \sum_{i \in N(t)} \frac{y_i}{|N(t)|}\]

A média é a escolha ótima porque, dentro de uma região, a constante \(c\) que minimiza \(\sum_{i}(y_i - c)^2\) é justamente \(\bar y\).

Limitação inerente: É importante notar que as predições dentro de cada região de uma árvore são constantes. Todas as observações de teste que caírem na mesma folha terminal receberão exatamente a mesma predição numérica, independentemente de pequenas variações nos valores de seus preditores individuais. A função estimada é, portanto, uma função em degraus, que aproxima mal relações suaves e lineares e não extrapola além do intervalo dos dados de treino.

15.10.3 Poda por Custo-Complexidade

Árvores muito profundas sobreajustam. Assim como na classificação, usa-se a poda por custo-complexidade (cost-complexity pruning): para cada \(\alpha \geq 0\), procura-se a subárvore \(T\) que minimiza

\[ \sum_{m=1}^{|T|}\sum_{i:\, \mathbf{x}_i \in R_m}(y_i - \hat y_{R_m})^2 + \alpha\,|T|, \]

em que \(|T|\) é o número de folhas. O parâmetro \(\alpha\) (o cp do pacote rpart, em escala relativa) é escolhido por validação cruzada. Note a semelhança com a regularização: perda + penalidade pela complexidade.

---

15.11 Florestas Aleatórias para Regressão

Assim como na classificação, a Floresta Aleatória para Regressão (Breiman 2001) combina o poder preditivo de centenas de árvores de regressão independentes geradas por meio de amostragem por bootstrap e seleção aleatória de atributos.

A previsão numérica final do conjunto da floresta é obtida calculando-se a média aritmética das previsões individuais de todas as árvores de regressão presentes no comitê:

\[ \hat f_{RF}(\mathbf{x}) = \frac{1}{B}\sum_{b=1}^B \hat T_b(\mathbf{x}). \]

Esse processo de média reduz significativamente a variância inerente das árvores individuais, resultando em previsões numéricas altamente estáveis e robustas. Formalmente, se cada árvore tem variância \(\sigma^2\) e a correlação entre pares de árvores é \(\rho\), a variância da média é

\[ \operatorname{Var}\big(\hat f_{RF}(\mathbf{x})\big) = \rho\,\sigma^2 + \frac{1 - \rho}{B}\,\sigma^2. \]

Aumentar \(B\) elimina o segundo termo, mas não o primeiro; é por isso que a floresta sorteia um subconjunto de mtry variáveis em cada divisão: isso descorrelaciona as árvores (reduz \(\rho\)) (Hastie, Tibshirani, e Friedman 2009). Para regressão, o padrão recomendado é mtry \(= \lfloor p/3 \rfloor\) (no pacote randomForest), contra \(\sqrt{p}\) na classificação.

Outros recursos úteis da floresta para regressão:

• Erro out-of-bag (OOB): cada árvore deixa de fora cerca de \(1/e \approx 36{,}8\%\) das observações; prever cada observação apenas com as árvores que não a viram fornece uma estimativa do erro de teste sem precisar de validação cruzada.
• Importância das variáveis: medida pelo aumento do MSE OOB quando os valores de uma variável são permutados aleatoriamente (%IncMSE), ou pela redução total do RSS proporcionada pelas divisões nessa variável (IncNodePurity).
• Limitação: como média de funções em degraus, a floresta também não extrapola: suas previsões ficam sempre dentro do intervalo dos \(y\) de treino.

---

15.12 k-NN para Regressão

O algoritmo k-NN (Capítulo de k-NN) também se adapta à regressão: em vez da votação majoritária, a previsão é a média das respostas dos \(k\) vizinhos mais próximos:

\[ \hat f(\mathbf{x}_0) = \frac{1}{k}\sum_{i \in \mathcal{N}_k(\mathbf{x}_0)} y_i, \]

em que \(\mathcal{N}_k(\mathbf{x}_0)\) é o conjunto dos índices dos \(k\) pontos de treino mais próximos de \(\mathbf{x}_0\). Uma variante pondera os vizinhos pelo inverso da distância. Valem as mesmas observações do caso de classificação: \(k\) pequeno gera alta variância, \(k\) grande gera alto viés, e a padronização das variáveis é obrigatória. Em alta dimensão, o k-NN sofre com a maldição da dimensionalidade, e modelos paramétricos costumam levar vantagem (James et al. 2013).

Exemplo: com os dados \(x = (1,2,3,4,5)\), \(y = (2,4,5,4,5)\) e \(k = 2\), para \(x_0 = 2{,}4\) os vizinhos são \(x = 2\) (distância 0,4) e \(x = 3\) (distância 0,6), e \(\hat f(2{,}4) = (4 + 5)/2 = 4{,}5\).

---

15.13 Métricas de Avaliação de Regressão

Para avaliar a acurácia dos modelos numéricos de regressão, comparamos os valores reais \(y_i\) com as previsões numéricas \(\hat{y}_i\) por meio das seguintes métricas estatísticas clássicas. Para avaliar a capacidade de generalização, elas devem ser calculadas no conjunto de teste ou por validação cruzada, nunca apenas no treino (Capítulo de Aprendizado de Máquina).

15.13.1 Erro Médio Absoluto (MAE - Mean Absolute Error)

Mapeia a média das diferenças absolutas (sem sinal) entre a previsão e o real, sendo de fácil interpretação física na mesma unidade de medida dos dados originais: \[MAE = \frac{1}{n} \sum_{i=1}^{n} |y_i - \hat{y}_i|\]

O MAE é mais robusto a outliers que as métricas quadráticas. A previsão constante que minimiza o MAE é a mediana de \(y\).

15.13.2 Erro Médio Quadrático (MSE - Mean Squared Error)

É a média dos erros ao quadrado, a função de perda minimizada pelos mínimos quadrados: \[MSE = \frac{1}{n} \sum_{i=1}^{n} (y_i - \hat{y}_i)^2\]

Sua unidade é o quadrado da unidade de \(y\), o que dificulta a interpretação direta. A previsão constante que minimiza o MSE é a média de \(y\). O MSE de teste esperado admite a famosa decomposição viés-variância: \(\mathbb{E}[(Y - \hat f(\mathbf{x}_0))^2] = \operatorname{Var}(\hat f(\mathbf{x}_0)) + [\operatorname{Vi\acute{e}s}(\hat f(\mathbf{x}_0))]^2 + \operatorname{Var}(\varepsilon)\) (James et al. 2013).

15.13.3 Raiz do Erro Médio Quadrático (RMSE - Root Mean Squared Error)

Mapeia a raiz quadrada da média das diferenças ao quadrado. Como eleva os desvios individuais ao quadrado antes de tirar a média, o RMSE penaliza de forma muito mais rigorosa grandes erros ou desvios isolados (outliers): \[RMSE = \sqrt{\frac{1}{n} \sum_{i=1}^{n} (y_i - \hat{y}_i)^2}\]

Sempre vale \(MAE \leq RMSE\); uma grande diferença entre ambos sinaliza a presença de alguns erros muito grandes.

15.13.4 Erro Percentual Absoluto Médio (MAPE - Mean Absolute Percentage Error)

Expressa o erro em termos relativos (percentuais), útil para comunicar resultados a públicos não técnicos e comparar séries de escalas diferentes: \[MAPE = \frac{100\%}{n}\sum_{i=1}^n \left|\frac{y_i - \hat y_i}{y_i}\right|\]

Limitações: não é definido quando algum \(y_i = 0\), explode quando \(y_i\) é próximo de zero e penaliza de forma assimétrica previsões acima e abaixo do valor real.

15.13.5 Coeficiente de Determinação (\(R^2\))

Mede a proporção da variância total da variável de resposta \(y\) que é explicada e capturada pelas colunas preditoras do modelo. Varia idealmente de 0 a 1, onde quanto mais próximo de 1, melhor o poder explanativo das variáveis do modelo: \[R^2 = 1 - \frac{SQ_{Res}}{SQ_{Tot}} = 1 - \frac{\sum_{i=1}^n (y_i - \hat y_i)^2}{\sum_{i=1}^n (y_i - \bar y)^2}.\]

Interpretação: o \(R^2\) compara o modelo com a previsão “ingênua” que sempre prevê \(\bar y\). Fora da amostra de treino (ou para modelos sem intercepto), o \(R^2\) pode ser negativo, indicando que o modelo é pior do que simplesmente prever a média.

Exemplo numérico: calculando as métricas

Para a reta \(\hat y = 2{,}2 + 0{,}6x\) ajustada aos dados do exemplo (\(e = (-0{,}8;\; 0{,}6;\; 1{,}0;\; -0{,}6;\; -0{,}2)\)):

• \(MAE = (0{,}8 + 0{,}6 + 1{,}0 + 0{,}6 + 0{,}2)/5 = 3{,}2/5 = 0{,}64\).
• \(MSE = 2{,}4/5 = 0{,}48\) e \(RMSE = \sqrt{0{,}48} \approx 0{,}693\).
• \(MAPE = \frac{100\%}{5}\left(\frac{0{,}8}{2} + \frac{0{,}6}{4} + \frac{1{,}0}{5} + \frac{0{,}6}{4} + \frac{0{,}2}{5}\right) = \frac{100\%}{5}(0{,}94) = 18{,}8\%\).
• \(R^2 = 1 - 2{,}4/6 = 0{,}60\).

Métrica	Unidade	Sensível a outliers?	Melhor valor	Observação
MAE	mesma de \(y\)	Pouco	0	Robusta, fácil de interpretar
MSE	quadrado de \(y\)	Muito	0	Diferenciável; base do MQO
RMSE	mesma de \(y\)	Muito	0	Mais usada em competições
MAPE	%	Pouco	0	Indefinida se \(y_i = 0\)
\(R^2\)	adimensional	Muito	1	Comparação com a média

---

15.14 Implementação Prática em R

15.14.1 Regressão Linear, Diagnóstico e VIF

# install.packages(c("car", "glmnet", "rpart", "randomForest", "MASS"))
data(mtcars)

# Regressão linear simples: consumo (mpg) em função do peso (wt, em 1000 lb)
m1 <- lm(mpg ~ wt, data = mtcars)
summary(m1)
#> Coeficientes: (Intercept) 37.285, wt -5.344; R2 = 0.753

# Verificação "à mão" da forma matricial beta = (X'X)^{-1} X'y
X <- model.matrix(m1)
y <- mtcars$mpg
solve(t(X) %*% X, t(X) %*% y)    # mesmos coeficientes de coef(m1)

# Regressão múltipla
m2 <- lm(mpg ~ wt + hp + qsec, data = mtcars)
summary(m2)                       # testes t, teste F, R2 e R2 ajustado
confint(m2)                       # intervalos de confiança de 95%

# Diagnóstico gráfico dos resíduos (4 gráficos)
par(mfrow = c(2, 2)); plot(m2); par(mfrow = c(1, 1))

# Multicolinearidade
library(car)
vif(m2)

# Variável categórica (dummies automáticas) e interação
mtcars$cyl_f <- factor(mtcars$cyl)          # referência: 4 cilindros
mtcars$am_f  <- factor(mtcars$am, labels = c("automatico", "manual"))
m3 <- lm(mpg ~ wt * am_f + cyl_f, data = mtcars)  # wt + am_f + wt:am_f + cyl_f
summary(m3)

# Regressão polinomial com polinômios ortogonais e comparação por ANOVA
m4a <- lm(mpg ~ hp, data = mtcars)
m4b <- lm(mpg ~ poly(hp, 2), data = mtcars)
anova(m4a, m4b)                   # o termo quadrático melhora o ajuste?

# Previsão com intervalo de predição para um carro novo
predict(m2, newdata = data.frame(wt = 3, hp = 150, qsec = 18),
        interval = "prediction")

15.14.2 Regularização com glmnet

library(glmnet)
library(MASS)                     # base Boston: valor mediano de imóveis (medv)
data(Boston)

set.seed(123)
idx    <- sample(nrow(Boston), size = 0.7 * nrow(Boston))
X      <- model.matrix(medv ~ ., data = Boston)[, -1]   # remove o intercepto
y      <- Boston$medv
X_tr <- X[idx, ];  y_tr <- y[idx]
X_te <- X[-idx, ]; y_te <- y[-idx]

# glmnet padroniza os preditores internamente (standardize = TRUE)
cv_ridge <- cv.glmnet(X_tr, y_tr, alpha = 0)     # Ridge
cv_lasso <- cv.glmnet(X_tr, y_tr, alpha = 1)     # Lasso
cv_enet  <- cv.glmnet(X_tr, y_tr, alpha = 0.5)   # Elastic Net

plot(cv_lasso)                                   # erro de VC vs log(lambda)
plot(cv_lasso$glmnet.fit, xvar = "lambda")       # caminho dos coeficientes
coef(cv_lasso, s = "lambda.1se")                 # repare nos zeros

rmse <- function(y, yhat) sqrt(mean((y - yhat)^2))
sapply(list(ridge = cv_ridge, lasso = cv_lasso, enet = cv_enet),
       function(m) rmse(y_te, predict(m, X_te, s = "lambda.min")))

15.14.3 Regressão Logística com glm

# Probabilidade de câmbio manual (am = 1) em função do peso
m_log <- glm(am ~ wt, data = mtcars, family = binomial)
summary(m_log)
exp(coef(m_log))                  # razões de chances (odds ratios)

p_hat <- predict(m_log, type = "response")       # probabilidades
y_hat <- ifelse(p_hat >= 0.5, 1, 0)
table(Real = mtcars$am, Previsto = y_hat)        # matriz de confusão
logLik(m_log)                                    # log-verossimilhança maximizada

15.14.4 Inferência, Influência, Heterocedasticidade e Seleção

# install.packages(c("lmtest", "sandwich"))
library(lmtest); library(sandwich)
data(mtcars)
# Exemplo do texto: inferência, intervalos e influência
x <- 1:5; y <- c(2, 4, 5, 4, 5)
m <- lm(y ~ x)
confint(m)                                      # IC de 95% para beta0 e beta1
predict(m, data.frame(x = 4), interval = "confidence")   # [3,04; 6,16]
predict(m, data.frame(x = 4), interval = "prediction")   # [1,35; 7,85]
hatvalues(m)                                    # 0,6 0,3 0,2 0,3 0,6
rstandard(m); rstudent(m)                       # resíduos padronizados e estudentizados
cooks.distance(m)                               # D_1 = 1,5
influence.measures(m)                           # quadro completo de diagnósticos

# Teste F para modelos encaixados
m_red  <- lm(mpg ~ wt, data = mtcars)
m_comp <- lm(mpg ~ wt + qsec + am, data = mtcars)
anova(m_red, m_comp)                            # H0: beta_qsec = beta_am = 0

# Heterocedasticidade: Breusch-Pagan e erros padrão robustos (HC3)
m2 <- lm(mpg ~ wt + hp, data = mtcars)
bptest(m2)                                      # BP = 0,88; p = 0,64
coeftest(m2, vcov = vcovHC(m2, type = "HC3"))
# Mínimos quadrados ponderados (ex.: desvio padrão proporcional a wt)
m_mqp <- lm(mpg ~ wt + hp, data = mtcars, weights = 1 / wt^2)
summary(m_mqp)

# Seleção de modelos: AIC, BIC e stepwise
completo <- lm(mpg ~ wt + hp + qsec + am + disp + drat, data = mtcars)
AIC(m_red, m_comp, completo); BIC(m_red, m_comp, completo)
step(lm(mpg ~ 1, data = mtcars), scope = formula(completo),
     direction = "both")                        # stepwise por AIC
step(completo, direction = "backward", k = log(nrow(mtcars)))   # por BIC

15.14.5 Logística via IRLS, MLGs, Splines e GAMs

library(MASS); library(splines); library(mgcv)
# IRLS "à mão" para a regressão logística do exemplo do texto
x <- 1:6; y <- c(0, 0, 1, 0, 1, 1)
X <- cbind(1, x); beta <- c(0, 0)
for (t in 1:10) {
  p    <- as.vector(1 / (1 + exp(-X %*% beta)))
  W    <- diag(p * (1 - p))
  beta <- beta + solve(t(X) %*% W %*% X, t(X) %*% (y - p))
}
print(drop(beta))                               # -4,249  1,214
m_log <- glm(y ~ x, family = binomial)
anova(m_log, test = "Chisq")                    # TRV: 3,362; p = 0,067

# Regressão de Poisson e superdispersão
m_pois <- glm(breaks ~ wool + tension, family = poisson, data = warpbreaks)
summary(m_pois)
exp(coef(m_pois))                               # razões de taxas
sum(residuals(m_pois, type = "pearson")^2) / df.residual(m_pois)   # phi = 4,26
m_qp <- glm(breaks ~ wool + tension, family = quasipoisson, data = warpbreaks)
m_nb <- glm.nb(breaks ~ wool + tension, data = warpbreaks)        # binomial negativa
AIC(m_pois, m_nb)

# Splines e GAM (base Boston: medv em função de lstat e rm)
m_ns  <- lm(medv ~ ns(lstat, df = 4) + rm, data = Boston)        # spline natural
m_gam <- gam(medv ~ s(lstat) + s(rm), data = Boston)             # mgcv
summary(m_gam)                                  # gl efetivos de cada termo
plot(m_gam, pages = 1, residuals = TRUE)        # funções f_j estimadas
AIC(lm(medv ~ lstat + rm, data = Boston), m_ns, m_gam)

15.14.6 Árvores e Florestas para Regressão

library(rpart)
library(randomForest)

treino <- Boston[idx, ]; teste <- Boston[-idx, ]

# Árvore de regressão (method = "anova" usa o critério de RSS)
arv <- rpart(medv ~ ., data = treino, method = "anova",
             control = rpart.control(cp = 0.001))
printcp(arv)                                     # erro de VC por cp
cp_otimo <- arv$cptable[which.min(arv$cptable[, "xerror"]), "CP"]
arv_podada <- prune(arv, cp = cp_otimo)

# Floresta aleatória (mtry padrão para regressão: p/3)
set.seed(123)
rf <- randomForest(medv ~ ., data = treino, ntree = 500, importance = TRUE)
rf                                              # MSE OOB e % de variância explicada
importance(rf)                                  # %IncMSE e IncNodePurity
varImpPlot(rf)

# Comparação no teste
metricas <- function(y, yhat) c(
  MAE  = mean(abs(y - yhat)),
  RMSE = sqrt(mean((y - yhat)^2)),
  MAPE = 100 * mean(abs((y - yhat) / y)),
  R2   = 1 - sum((y - yhat)^2) / sum((y - mean(y))^2)
)
lin <- lm(medv ~ ., data = treino)
rbind(
  linear = metricas(teste$medv, predict(lin, teste)),
  arvore = metricas(teste$medv, predict(arv_podada, teste)),
  floresta = metricas(teste$medv, predict(rf, teste))
)

---

15.15 Implementação Prática em Python

15.15.1 Regressão Linear e Regularizada

import numpy as np
import pandas as pd
from sklearn.datasets import load_diabetes
from sklearn.model_selection import train_test_split, cross_val_score, KFold
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression, RidgeCV, LassoCV, ElasticNetCV
from sklearn.metrics import (mean_absolute_error, mean_squared_error,
                             mean_absolute_percentage_error, r2_score)

# Dados: progressão do diabetes após um ano (442 pacientes, 10 preditores)
X, y = load_diabetes(return_X_y=True, as_frame=True)
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.3, random_state=42)

# Mínimos quadrados "à mão": beta = (X'X)^{-1} X'y  (usando solve, mais estável)
Xd = np.column_stack([np.ones(len(X_tr)), X_tr.to_numpy()])
beta = np.linalg.solve(Xd.T @ Xd, Xd.T @ y_tr.to_numpy())

ols = LinearRegression().fit(X_tr, y_tr)
print(np.allclose(beta, np.r_[ols.intercept_, ols.coef_]))   #> True

alphas = np.logspace(-3, 3, 100)   # no sklearn, "alpha" é o lambda do texto
modelos = {
    "MQO":         make_pipeline(StandardScaler(), LinearRegression()),
    "Ridge":       make_pipeline(StandardScaler(), RidgeCV(alphas=alphas)),
    "Lasso":       make_pipeline(StandardScaler(), LassoCV(cv=5, random_state=0)),
    "Elastic Net": make_pipeline(StandardScaler(),
                                 ElasticNetCV(l1_ratio=[0.2, 0.5, 0.8], cv=5,
                                              random_state=0)),
}

def metricas(y, yhat):
    return {"MAE": mean_absolute_error(y, yhat),
            "RMSE": np.sqrt(mean_squared_error(y, yhat)),
            "MAPE": 100 * mean_absolute_percentage_error(y, yhat),
            "R2": r2_score(y, yhat)}

resultados = {}
for nome, m in modelos.items():
    m.fit(X_tr, y_tr)
    resultados[nome] = metricas(y_te, m.predict(X_te))
print(pd.DataFrame(resultados).T.round(3))

# Coeficientes do Lasso: quais variáveis foram zeradas?
lasso = modelos["Lasso"][-1]
print("lambda escolhido:", lasso.alpha_)
print(pd.Series(lasso.coef_, index=X.columns).round(2))

Para uma saída no estilo “tabela de regressão” (erros padrão, testes \(t\), \(R^2\) ajustado) e para o VIF, o pacote statsmodels é mais adequado:

import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor

Xc = sm.add_constant(X_tr)
ajuste = sm.OLS(y_tr, Xc).fit()
print(ajuste.summary())

vif = pd.Series([variance_inflation_factor(Xc.values, j)
                 for j in range(1, Xc.shape[1])], index=X_tr.columns)
print(vif.round(2))

# Gráfico de resíduos vs. ajustados
import matplotlib.pyplot as plt
plt.scatter(ajuste.fittedvalues, ajuste.resid, alpha=0.6)
plt.axhline(0, color="red", linestyle="--")
plt.xlabel("Valores ajustados"); plt.ylabel("Resíduos")
plt.show()

15.15.2 Inferência, Diagnóstico e MLGs com statsmodels

import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf
from statsmodels.stats.diagnostic import het_breuschpagan

# Exemplo do texto: inferência, intervalos e influência
x = np.array([1, 2, 3, 4, 5]); y = np.array([2, 4, 5, 4, 5])
m = sm.OLS(y, sm.add_constant(x)).fit()
print(m.conf_int())                              # IC de 95% para beta0 e beta1
prev = m.get_prediction(np.array([[1, 4]])).summary_frame(alpha=0.05)
print(prev[["mean", "mean_ci_lower", "mean_ci_upper",
            "obs_ci_lower", "obs_ci_upper"]])    # IC da média e de predição
infl = m.get_influence()
print(infl.hat_matrix_diag)                      # alavancas
print(infl.cooks_distance[0])                    # distâncias de Cook

# Heterocedasticidade e erros padrão robustos
mtcars = sm.datasets.get_rdataset("mtcars").data
m2 = smf.ols("mpg ~ wt + hp", data=mtcars).fit()
lm_stat, lm_p, _, _ = het_breuschpagan(m2.resid, m2.model.exog)
print(f"Breusch-Pagan: LM = {lm_stat:.3f}, p = {lm_p:.3f}")
print(m2.get_robustcov_results(cov_type="HC3").bse)   # EPs robustos
m_mqp = smf.wls("mpg ~ wt + hp", data=mtcars, weights=1 / mtcars.wt**2).fit()
print(m_mqp.params)

# Critérios de informação
for f in ["mpg ~ wt", "mpg ~ wt + hp", "mpg ~ wt + qsec + am"]:
    r = smf.ols(f, data=mtcars).fit()
    print(f"{f:25s} AIC = {r.aic:7.2f}  BIC = {r.bic:7.2f}")
# Obs.: o statsmodels não conta sigma como parâmetro; seus AIC/BIC ficam
# 2 (ou log n) unidades abaixo dos do R, sem alterar a ordenação dos modelos.

# Regressão logística via IRLS (GLM) e teste da razão de verossimilhanças
xl = np.arange(1, 7); yl = np.array([0, 0, 1, 0, 1, 1])
glm_log = sm.GLM(yl, sm.add_constant(xl), family=sm.families.Binomial()).fit()
print(glm_log.params, glm_log.deviance, glm_log.null_deviance)
trv = glm_log.null_deviance - glm_log.deviance
from scipy.stats import chi2
print(f"TRV = {trv:.3f}, p = {chi2.sf(trv, 1):.3f}")

# Regressão de Poisson
wb = sm.datasets.get_rdataset("warpbreaks").data
pois = smf.glm("breaks ~ wool + tension", data=wb,
               family=sm.families.Poisson()).fit()
print(np.exp(pois.params))   # razões de taxas (referência: tensão H, a 1a em ordem alfabética)
print("phi =", pois.pearson_chi2 / pois.df_resid)    # superdispersão

15.15.3 Regressão Polinomial com Escolha do Grau por Validação Cruzada

from sklearn.preprocessing import PolynomialFeatures

rng = np.random.default_rng(1)
x = rng.uniform(-3, 3, 150).reshape(-1, 1)
y_poly = 1 + 0.5 * x.ravel() - 0.8 * x.ravel() ** 2 + rng.normal(0, 1, 150)

cv = KFold(n_splits=5, shuffle=True, random_state=0)
for grau in range(1, 9):
    modelo = make_pipeline(PolynomialFeatures(grau, include_bias=False),
                           StandardScaler(), LinearRegression())
    mse = -cross_val_score(modelo, x, y_poly, cv=cv,
                           scoring="neg_mean_squared_error").mean()
    print(f"grau {grau}: MSE de VC = {mse:.3f}")
# Espera-se o mínimo em torno do grau 2 (o verdadeiro)

15.15.4 Regressão Logística

from sklearn.datasets import load_breast_cancer
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import log_loss, accuracy_score

Xb, yb = load_breast_cancer(return_X_y=True)
Xb_tr, Xb_te, yb_tr, yb_te = train_test_split(Xb, yb, test_size=0.3,
                                              stratify=yb, random_state=0)
logit = make_pipeline(StandardScaler(), LogisticRegression(max_iter=1000))
logit.fit(Xb_tr, yb_tr)

p = logit.predict_proba(Xb_te)[:, 1]          # P(Y = 1 | x)
print("Log-loss (entropia cruzada):", round(log_loss(yb_te, p), 4))
print("Acurácia (limiar 0,5):", round(accuracy_score(yb_te, p >= 0.5), 4))
print("Razões de chances (5 primeiras):",
      np.exp(logit[-1].coef_[0][:5]).round(3))

15.15.5 Árvore, Floresta e k-NN para Regressão

from sklearn.tree import DecisionTreeRegressor
from sklearn.ensemble import RandomForestRegressor
from sklearn.neighbors import KNeighborsRegressor
from sklearn.model_selection import GridSearchCV

arvore = GridSearchCV(DecisionTreeRegressor(random_state=0),
                      {"max_depth": [2, 3, 4, 5, 6, 8, None],
                       "min_samples_leaf": [1, 5, 10, 20]},
                      cv=5, scoring="neg_root_mean_squared_error")

floresta = RandomForestRegressor(n_estimators=500, max_features=1/3,
                                 oob_score=True, random_state=0, n_jobs=-1)

knn = GridSearchCV(make_pipeline(StandardScaler(), KNeighborsRegressor()),
                   {"kneighborsregressor__n_neighbors": range(1, 41)},
                   cv=5, scoring="neg_root_mean_squared_error")

for nome, m in {"Árvore": arvore, "Floresta": floresta, "k-NN": knn}.items():
    m.fit(X_tr, y_tr)
    resultados[nome] = metricas(y_te, m.predict(X_te))

print(pd.DataFrame(resultados).T.round(3))
print("R2 out-of-bag da floresta:", round(floresta.oob_score_, 3))
print(pd.Series(floresta.feature_importances_, index=X.columns)
        .sort_values(ascending=False).round(3))

Qual modelo escolher?

Não existe um vencedor universal. Em bases pequenas com relações aproximadamente lineares (como diabetes), modelos lineares regularizados costumam ser tão bons quanto ou melhores que florestas, além de muito mais interpretáveis. Em bases grandes com interações e não linearidades complexas, florestas e outros ensembles tendem a vencer. Compare sempre os candidatos com a mesma validação cruzada ou o mesmo conjunto de teste.

---

15.16 Resumo do Capítulo

• Regressão estima a função \(f(\mathbf{x}) = \mathbb{E}[Y \mid \mathbf{X} = \mathbf{x}]\); modelos paramétricos assumem uma forma para \(f\) e os não paramétricos a aprendem dos dados.
• Na regressão linear simples, \(\hat\beta_1 = S_{xy}/S_{xx}\) e \(\hat\beta_0 = \bar y - \hat\beta_1\bar x\); na múltipla, \(\hat{\boldsymbol\beta} = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{y}\), a projeção ortogonal de \(\mathbf{y}\) no espaço das colunas de \(\mathbf{X}\).
• Cada \(\beta_j\) é o efeito de \(x_j\) mantidas as demais variáveis constantes; categorias entram como \(k-1\) dummies e efeitos que dependem de outra variável, como interações.
• \(R^2 = 1 - SQ_{Res}/SQ_{Tot}\) mede a variabilidade explicada; o \(R^2\) ajustado penaliza o número de preditores. Resíduos verificam as suposições; o VIF detecta multicolinearidade.
• Ridge (\(\ell_2\)), Lasso (\(\ell_1\)) e Elastic Net encolhem coeficientes em troca de menor variância; o Lasso zera coeficientes e seleciona variáveis; \(\lambda\) é escolhido por validação cruzada.
• A regressão logística modela o log das chances como função linear dos preditores, é estimada por máxima verossimilhança e produz uma fronteira de decisão linear.
• Árvores de regressão dividem o espaço minimizando o RSS e preveem a média da folha; florestas reduzem a variância pela média de árvores descorrelacionadas; o k-NN prevê a média dos vizinhos.
• Métricas: MAE, MSE, RMSE, MAPE e \(R^2\), sempre avaliadas fora da amostra de treinamento.
• A matriz chapéu \(\mathbf{H}\) é uma projeção ortogonal (simétrica, idempotente, com traço \(p+1\)); daí seguem o não viés de \(s^2 = SQ_{Res}/(n-p-1)\), as alavancas \(h_{ii}\) e a noção de graus de liberdade.
• Pelo teorema de Gauss–Markov, o MQO é o melhor estimador linear não viesado; sob normalidade, \(\hat{\boldsymbol\beta}\) é normal e independente de \(s^2\), o que fundamenta os testes \(t\) e \(F\) e os intervalos de confiança e de predição (este último sempre mais largo).
• A distância de Cook combina resíduo e alavanca para detectar pontos influentes; heterocedasticidade é diagnosticada por Breusch–Pagan e tratada por erros padrão robustos (sanduíche), MQP ou transformações.
• \(C_p\), AIC e BIC corrigem o otimismo do erro de treino; o BIC penaliza mais e é consistente, o AIC visa à previsão. Seleção automática invalida a inferência clássica no modelo final.
• O Ridge encolhe as direções de menor variância dos preditores e sempre admite um \(\lambda > 0\) com EQM menor que o do MQO; o Lasso é resolvido por limiar suave e descida coordenada. Ridge e Lasso são estimadores MAP sob prioris gaussiana e de Laplace.
• A regressão logística é ajustada por Newton–Raphson/IRLS e comparada por deviance e teste da razão de verossimilhanças; ela, a regressão normal e a de Poisson são casos dos Modelos Lineares Generalizados.
• Splines e GAMs estendem o modelo linear a relações suaves não lineares, mantendo interpretabilidade e o controle da complexidade por graus de liberdade efetivos.

---

15.17 Exercícios Propostos

1. Para os dados \(x = (2, 4, 6, 8)\) e \(y = (3, 7, 8, 12)\), calcule à mão \(\hat\beta_0\), \(\hat\beta_1\), os resíduos, \(SQ_{Res}\), \(SQ_{Tot}\) e \(R^2\). Verifique que \(\sum e_i = 0\) e \(\sum x_i e_i = 0\).
2. Mostre que, na regressão linear simples, \(\hat\beta_1 = r_{xy}\,s_y/s_x\) e que \(R^2 = r_{xy}^2\).
3. Derive as equações normais \(\mathbf{X}^\top\mathbf{X}\hat{\boldsymbol\beta} = \mathbf{X}^\top\mathbf{y}\) e mostre que a matriz chapéu \(\mathbf{H}\) é simétrica e idempotente (\(\mathbf{H}^2 = \mathbf{H}\)). Interprete geometricamente.
4. Mostre que \(\hat{\boldsymbol\beta}\) é não viesado e que \(\operatorname{Var}(\hat{\boldsymbol\beta}) = \sigma^2(\mathbf{X}^\top\mathbf{X})^{-1}\) sob as suposições do modelo linear.
5. Um modelo com \(n = 50\) e \(p = 6\) preditores tem \(R^2 = 0{,}72\). Calcule o \(R^2\) ajustado. Se acrescentarmos 4 preditores irrelevantes e o \(R^2\) subir para 0,73, o que acontece com o \(R^2\) ajustado?
6. Se a regressão de \(x_3\) sobre os demais preditores tem \(R_3^2 = 0{,}95\), qual é o VIF de \(x_3\)? Por qual fator o erro padrão de \(\hat\beta_3\) fica inflado?
7. Em um modelo \(\text{salario} = \beta_0 + \beta_1\,\text{exp} + \beta_2 D + \beta_3\,(\text{exp}\times D) + \varepsilon\), em que \(D = 1\) para pós-graduados, interprete cada coeficiente e escreva as retas para os dois grupos.
8. Mostre que, no caso ortonormal, a solução Ridge é \(\hat\beta_j^{\text{MQO}}/(1+\lambda)\) e que a do Lasso é o limiar suave. Aplique para \(\hat{\boldsymbol\beta}^{\text{MQO}} = (2;\; -0{,}3;\; 1)\) e \(\lambda = 1\).
9. Em uma regressão logística, \(\hat\beta_0 = -4\) e \(\hat\beta_1 = 0{,}8\). Calcule \(\hat p\) para \(x = 3\) e \(x = 7\), o valor de \(x\) em que \(\hat p = 0{,}5\) e a razão de chances associada a um aumento de 2 unidades em \(x\).
10. Mostre que o gradiente da log-verossimilhança da regressão logística é \(\mathbf{X}^\top(\mathbf{y} - \mathbf{p})\).
11. Para \(x = (1, 2, 3, 4, 5, 6)\) e \(y = (1, 2, 2, 8, 9, 10)\), encontre o melhor primeiro corte de uma árvore de regressão calculando o RSS de todos os cortes candidatos.
12. Calcule MAE, RMSE, MAPE e \(R^2\) para \(y = (10, 12, 15, 20)\) e \(\hat y = (11, 11, 16, 17)\).
13. 18. Com a base Boston (pacote MASS), compare por validação cruzada de 10 partições a regressão linear, o Lasso, uma árvore podada e uma floresta aleatória para prever medv. Analise os resíduos da regressão linear e o VIF das variáveis.
14. (Python) Com load_diabetes, trace o caminho dos coeficientes do Lasso em função de \(\lambda\) (use sklearn.linear_model.lasso_path) e identifique a ordem em que as variáveis entram no modelo.
15. (Python) Simule \(y = \sin(2x) + \varepsilon\) para \(x \in [0, 3]\) e compare regressão polinomial, árvore, floresta e k-NN. Em seguida, preveja para \(x \in [3, 4]\) (extrapolação) e comente o comportamento de cada modelo.
16. Mostre que \(\operatorname{tr}(\mathbf{H}) = p + 1\), que os autovalores de \(\mathbf{H}\) são 0 ou 1 e que \(1/n \leq h_{ii} \leq 1\) em um modelo com intercepto. Na regressão simples, prove que \(h_{ii} = 1/n + (x_i - \bar x)^2/S_{xx}\).
17. Reproduza a prova do teorema de Gauss–Markov. Em seguida, mostre por um exemplo que o estimador Ridge (que é linear em \(\mathbf{y}\)) não contradiz o teorema.
18. Para os dados do exercício 1, construa o IC de 95% para \(\beta_1\), o IC para a resposta média em \(x_0 = 5\) e o intervalo de predição em \(x_0 = 5\). Explique por que o segundo é mais largo.
19. Para os dados do exercício 1, calcule as alavancas, os resíduos padronizados e as distâncias de Cook. Alguma observação é influente?
20. Mostre que, sob \(\operatorname{Var}(\varepsilon_i) = \sigma^2/w_i\), o estimador de mínimos quadrados ponderados é \((\mathbf{X}^\top\mathbf{W}\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{W}\mathbf{y}\) e que ele é o BLUE. Qual é a variância do MQO nesse caso?
21. Em um problema com \(n = 100\), compare as penalidades por parâmetro do AIC e do BIC. A partir de que tamanho amostral o BIC penaliza mais que o AIC?
22. Mostre que a moda a posteriori de \(\boldsymbol\beta\) sob a priori \(N(\mathbf{0}, \tau^2\mathbf{I})\) é o estimador Ridge com \(\lambda = \sigma^2/\tau^2\). Com um preditor, \(S_{xx} = 20\), \(\beta = 1\) e \(\sigma^2 = 4\), calcule \(\lambda^*\) e o EQM do Ridge e do MQO.
23. Usando o subgradiente, derive a solução de \(\min_\beta\, (z - \beta)^2 + \lambda|\beta|\) e esboce o gráfico de \(\hat\beta\) em função de \(z\) para o Lasso e para o Ridge.
24. Para \(x = (1, 2, 3, 4)\) e \(y = (0, 1, 0, 1)\), faça à mão a primeira iteração do IRLS a partir de \(\boldsymbol\beta^{(0)} = \mathbf{0}\). Em seguida, use glm para obter o EMV, a deviance e o TRV para \(H_0: \beta_1 = 0\).
25. Escreva as distribuições Poisson e Bernoulli na forma da família exponencial, identificando \(\theta\), \(b(\theta)\) e \(a(\phi)\), e verifique que \(b'(\theta) = \mu\) e \(b''(\theta) = V(\mu)\).
26. 18. Com a base Boston, compare por validação cruzada uma regressão linear em lstat, uma polinomial de grau 3, uma spline natural com 4 graus de liberdade e um GAM s(lstat). Discuta o comportamento de cada ajuste nas extremidades de lstat.

---

15.18 Leituras Recomendadas

• James et al. (2013) (James et al. 2013), Capítulos 3, 4, 6 e 8: regressão linear, regressão logística, regularização e árvores, com laboratórios em R.
• Hastie, Tibshirani e Friedman (2009) (Hastie, Tibshirani, e Friedman 2009), Capítulos 3, 4, 9 e 15: tratamento aprofundado de métodos lineares, árvores e florestas.
• Montgomery, Peck e Vining (2012) (Montgomery, Peck, e Vining 2012): referência clássica de análise de regressão, com ênfase em diagnóstico e multicolinearidade.
• Morettin e Singer (2021) (Morettin e Singer 2021): Estatística e Ciência de Dados em português, com capítulos sobre regressão e regularização.
• Hoerl e Kennard (1970) (Hoerl e Kennard 1970), Tibshirani (1996) (Tibshirani 1996) e Zou e Hastie (2005) (Zou e Hastie 2005): artigos originais de Ridge, Lasso e Elastic Net.
• Friedman, Hastie e Tibshirani (2010) (Friedman, Hastie, e Tibshirani 2010): algoritmo e pacote glmnet.
• Breiman et al. (1984) (Breiman et al. 1984) e Breiman (2001) (Breiman 2001): árvores CART e florestas aleatórias.
• Galton (1886) (Galton 1886): a origem histórica do termo “regressão”.
• Seber e Lee (2003) (Seber e Lee 2003): tratamento rigoroso da teoria do modelo linear (projeções, Gauss–Markov, distribuições de formas quadráticas, testes \(F\)).
• Cook (1977) (Cook 1977), White (1980) (White 1980) e Breusch e Pagan (1979) (Breusch e Pagan 1979): artigos originais da distância de Cook, dos erros padrão robustos e do teste de heterocedasticidade.
• Mallows (1973) (Mallows 1973), Akaike (1974) (Akaike 1974) e Schwarz (1978) (Schwarz 1978): as origens de \(C_p\), AIC e BIC.
• McCullagh e Nelder (1989) (McCullagh e Nelder 1989): a referência clássica sobre Modelos Lineares Generalizados.
• Hastie e Tibshirani (1986) (Hastie e Tibshirani 1986) e Wood (2017) (Wood 2017): modelos aditivos generalizados, da proposta original ao pacote mgcv.
• Hastie, Tibshirani e Wainwright (2015) (Hastie, Tibshirani, e Wainwright 2015): teoria e algoritmos do Lasso e de métodos esparsos; Park e Casella (2008) (Park e Casella 2008) para a visão bayesiana.

Akaike, Hirotugu. 1974. “A New Look at the Statistical Model Identification”. IEEE Transactions on Automatic Control 19 (6): 716–23.

Breiman, Leo. 2001. “Random Forests”. Machine Learning 45 (1): 5–32.

Breiman, Leo, Jerome H. Friedman, Richard A. Olshen, e Charles J. Stone. 1984. Classification and Regression Trees. Belmont, CA: Wadsworth.

Breusch, Trevor S., e Adrian R. Pagan. 1979. “A Simple Test for Heteroscedasticity and Random Coefficient Variation”. Econometrica 47 (5): 1287–94.

Cook, R. Dennis. 1977. “Detection of Influential Observation in Linear Regression”. Technometrics 19 (1): 15–18.

Cox, David R. 1958. “The regression analysis of binary sequences”. Journal of the Royal Statistical Society: Series B (Methodological) 20 (2): 215–42.

Friedman, Jerome, Trevor Hastie, e Robert Tibshirani. 2010. “Regularization paths for generalized linear models via coordinate descent”. Journal of Statistical Software 33 (1): 1–22.

Galton, Francis. 1886. “Regression towards mediocrity in hereditary stature”. The Journal of the Anthropological Institute of Great Britain and Ireland 15: 246–63.

Hastie, Trevor, e Robert Tibshirani. 1986. “Generalized Additive Models”. Statistical Science 1 (3): 297–310.

Hastie, Trevor, Robert Tibshirani, e Jerome Friedman. 2009. The Elements of Statistical Learning: Data Mining, Inference, and Prediction. 2º ed. New York: Springer.

Hastie, Trevor, Robert Tibshirani, e Martin Wainwright. 2015. Statistical Learning with Sparsity: The Lasso and Generalizations. Boca Raton: Chapman; Hall/CRC.

Hoerl, Arthur E., e Robert W. Kennard. 1970. “Ridge Regression: Biased Estimation for Nonorthogonal Problems”. Technometrics 12 (1): 55–67.

James, Gareth, Daniela Witten, Trevor Hastie, e Robert Tibshirani. 2013. An Introduction to Statistical Learning: with Applications in R. Springer.

Mallows, Colin L. 1973. “Some Comments on \(C_p\)”. Technometrics 15 (4): 661–75.

McCullagh, Peter, e John A. Nelder. 1989. Generalized Linear Models. 2º ed. London: Chapman; Hall.

Montgomery, Douglas C., Elizabeth A. Peck, e G. Geoffrey Vining. 2012. Introduction to Linear Regression Analysis. 5º ed. Hoboken, NJ: Wiley.

Morettin, Pedro A, e Julio M Singer. 2021. Estatística e Ciência de Dados. Universidade de São Paulo.

Nelder, John A., e Robert W. M. Wedderburn. 1972. “Generalized linear models”. Journal of the Royal Statistical Society: Series A (General) 135 (3): 370–84.

Park, Trevor, e George Casella. 2008. “The Bayesian Lasso”. Journal of the American Statistical Association 103 (482): 681–86.

Schwarz, Gideon. 1978. “Estimating the Dimension of a Model”. The Annals of Statistics 6 (2): 461–64.

Seber, George A. F., e Alan J. Lee. 2003. Linear Regression Analysis. 2º ed. Hoboken: Wiley.

Tibshirani, Robert. 1996. “Regression Shrinkage and Selection via the Lasso”. Journal of the Royal Statistical Society: Series B (Methodological) 58 (1): 267–88.

White, Halbert. 1980. “A Heteroskedasticity-Consistent Covariance Matrix Estimator and a Direct Test for Heteroskedasticity”. Econometrica 48 (4): 817–38.

Wood, Simon N. 2017. Generalized Additive Models: An Introduction with R. 2º ed. Boca Raton: Chapman; Hall/CRC.

Zou, Hui, e Trevor Hastie. 2005. “Regularization and variable selection via the elastic net”. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 67 (2): 301–20.