Adequação do modelo, diagnósticos de influência e extensões
Prof. Dr. Juliano Bortolini
Depto. de Estatística — UFMT
Suposições do modelo de regressão
Até aqui, ao ajustar o modelo \(\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}\), assumimos que (Montgomery, Cap. 4.1):
A relação entre a resposta \(y\) e os regressores é linear (ao menos aproximadamente);
O erro \(\varepsilon\) tem média zero;
O erro \(\varepsilon\) tem variância constante\(\sigma^2\) (homocedasticidade);
Os erros são não correlacionados;
Os erros têm distribuição normal.
As suposições 4 e 5 juntas implicam que os erros são independentes. A suposição 5 é necessária para os testes de hipóteses e a construção de intervalos.
Devemos sempre considerar duvidosa a validade dessas suposições e conduzir análises para examinar a adequação do modelo proposto.
O que faremos nesta aula
Cap. 4 — Adequação do modelo: análise de resíduos, gráficos de resíduos, estatística PRESS e detecção de outliers;
Cap. 5 — Transformações: estabilização de variância e linearização;
Cap. 6 — Leverage e influência:hat matrix, distância de Cook, DFBETAS e DFFITS;
Cap. 7 — Regressão polinomial: modelos de 2ª ordem e interação;
Cap. 8 — Variáveis indicadoras: variáveis qualitativas na regressão.
Referência principal: Montgomery, D. C.; Peck, E. A.; Vining, G. G. Introduction to Linear Regression Analysis, 6ª ed. Wiley, 2021 (Cap. 4 a 8).
Definição de resíduo
O resíduo é a diferença entre o valor observado e o valor ajustado:
\[
e_i = y_i - \hat{y}_i, \qquad i = 1, 2, \ldots, n.
\]
O resíduo pode ser visto como o desvio entre o dado e o ajuste — logo, é uma medida da variabilidade da resposta não explicada pelo modelo. Qualquer desvio das suposições sobre os erros deve aparecer nos resíduos.
Os resíduos têm média zero e sua variância média aproximada é estimada por
Para comparar resíduos em escala comum e detectar observações atípicas, usamos versões escalonadas (Montgomery, 4.2.1):
Resíduo
Fórmula
Uso
Bruto
\(e_i = y_i - \hat{y}_i\)
base
Padronizado
\(d_i = \dfrac{e_i}{\sqrt{MS_{Res}}}\)
\(\lvert d_i\rvert > 3\) sugere outlier
Studentizado (interno)
\(r_i = \dfrac{e_i}{\sqrt{MS_{Res}(1-h_{ii})}}\)
variância constante
PRESS
\(e_{(i)} = \dfrac{e_i}{1-h_{ii}}\)
capacidade preditiva
R-Student (externo)
\(t_i = \dfrac{e_i}{\sqrt{S^2_{(i)}(1-h_{ii})}}\)
outlier
onde \(h_{ii}\) é o elemento \(i\) da diagonal da matriz \(\mathbf{H}\).
Resíduo padronizado
Como a variância média aproximada de um resíduo é estimada por \(MS_{Res}\), um escalonamento natural é o resíduo padronizado:
\[
d_i = \frac{e_i}{\sqrt{MS_{Res}}}, \qquad i = 1, 2, \ldots, n.
\]
Os resíduos padronizados têm média zero e variância aproximadamente unitária. Um resíduo padronizado grande (\(\lvert d_i \rvert > 3\), por exemplo) sinaliza um potencial outlier.
Resíduo studentizado
Usar \(MS_{Res}\) como variância do \(i\)-ésimo resíduo é apenas uma aproximação. Podemos melhorar o escalonamento usando o desvio-padrão exato do resíduo. Como \(\mathbf{e} = (\mathbf{I}-\mathbf{H})\mathbf{y}\),
\[
r_i = \frac{e_i}{\sqrt{MS_{Res}(1-h_{ii})}}, \qquad i = 1, 2, \ldots, n,
\]
que tem variância constante \(\mathrm{Var}(r_i)=1\) quando o modelo está correto.
Resíduo PRESS
Outra forma de tornar os resíduos úteis na detecção de outliers é examinar \(y_i - \hat{y}_{(i)}\), em que \(\hat{y}_{(i)}\) é o valor ajustado da \(i\)-ésima resposta com base em todas as observações, exceto a \(i\)-ésima:
O resíduo PRESS é o resíduo comum ponderado pelos elementos \(h_{ii}\) da hat matrix. Pontos com \(h_{ii}\) grande terão resíduos PRESS grandes — em geral, pontos de alta influência. Sua variância é
Usando essa estimativa obtemos o resíduo externamente studentizado (R-Student):
\[
t_i = \frac{e_i}{\sqrt{S^2_{(i)}(1-h_{ii})}}, \qquad i = 1, 2, \ldots, n.
\]
Se a \(i\)-ésima observação for influente, \(S^2_{(i)}\) difere bastante de \(MS_{Res}\), e o R-Student é mais sensível a esse ponto.
O conjunto de dados (Exemplo 3.1, p. 76)
Uma engarrafadora de refrigerantes analisa as rotas de atendimento a máquinas de venda. Deseja-se prever o tempo\(y\) (minutos) que o entregador gasta para abastecer as máquinas, em função do número de caixas\(x_1\) e da distância caminhada\(x_2\) (pés). São 25 observações.
Recomenda-se usar o resíduo studentizado ou R-Student. Espera-se que sigam, aproximadamente, uma distribuição \(t\) com \(n-p-1\) graus de liberdade (ou \(N(0,1)\) para \(n\) grande).
Por que representar graficamente os resíduos?
A análise gráfica dos resíduos é uma forma muito eficaz de investigar a adequação do ajuste e verificar as suposições. Em geral, representamos os resíduos externamente studentizados (que têm variância constante).
Três gráficos fundamentais:
Gráfico de probabilidade normal dos resíduos — verifica normalidade;
Resíduos vs. valores ajustados\(\hat{y}_i\) — verifica variância e linearidade;
Resíduos na ordem temporal — verifica independência (autocorrelação).
Gráfico de probabilidade normal (Figura 4.3)
Pequenos desvios da normalidade não afetam muito o modelo, mas não-normalidade grosseira é séria, pois as estatísticas \(t\) e \(F\) e os intervalos dependem dessa suposição.
Figura 1: Padrões típicos de gráficos de probabilidade normal (dados ilustrativos).
Resíduos vs. valores ajustados (Figura 4.5)
Um gráfico dos resíduos studentizados \(t_i\) contra os valores ajustados \(\hat{y}_i\) é útil para detectar várias inadequações do modelo.
Figura 2: Padrões de resíduos vs. ajustados: (a) satisfatório; (b) funil; (c) duplo arco; (d) não linear.
Os padrões (b), (c) e (d) indicam, respectivamente, variância não constante, variância que cresce e depois decresce, e não-linearidade (falta de termos no modelo).
Resíduos na ordem temporal — autocorrelação
Se a ordem de coleta dos dados é conhecida, convém representar os resíduos contra essa ordem. O ideal é uma faixa horizontal aleatória.
A correlação entre erros em períodos diferentes chama-se autocorrelação e é uma violação potencialmente séria das suposições básicas (tratada no Cap. 14).
Exemplo em R — gráficos de resíduos
qqnorm(res_3, pch =19, main ="Resíduos studentizados")qqline(res_3, col ="blue", lwd =2)
Figura 3: Gráfico de probabilidade normal dos resíduos studentizados (delivery time).
O PRESS mede quão bem o modelo prevê novos dados. Um modelo com PRESS pequeno é preferível.
Usando o PRESS para comparar modelos
No próprio delivery time, ao acrescentar \(x_2\) (distância) ao modelo que continha apenas \(x_1\) (caixas), o PRESS diminui — indicando melhor capacidade preditiva.
mod_x1 <-lm(y ~ x1, data = entrega) # só caixasmod # caixas + distância (já ajustado)#> #> Call:#> lm(formula = y ~ x1 + x2, data = entrega)#> #> Coefficients:#> (Intercept) x1 x2 #> 2.3412 1.6159 0.0144PRESS_x1 <-sum(rstandard(mod_x1, type ="pred")^2)PRESS_x12 <-sum(rstandard(mod, type ="pred")^2)c(PRESS_so_x1 = PRESS_x1, PRESS_x1_x2 = PRESS_x12)#> PRESS_so_x1 PRESS_x1_x2 #> 733.55 459.04
O PRESS cai de 733.55 (só \(x_1\)) para \(\approx\)459 (\(x_1+x_2\)) — o modelo de dois regressores prevê melhor. (Montgomery cita 733.55 e 457.40; a pequena diferença é de arredondamento.)
O que é um outlier
Um outlier é uma observação extrema, consideravelmente diferente das demais. Resíduos muito maiores em valor absoluto (três ou quatro desvios-padrão) indicam potenciais outliers no espaço de \(y\).
Os gráficos de resíduos vs. \(\hat{y}_i\) e o gráfico de probabilidade normal ajudam a identificá-los. Examinar resíduos escalonados (studentizado e R-Student) é uma excelente forma de detectar outliers.
Exemplo em R — identificando o outlier
plot(res_3, main ="Resíduos studentizados", pch =19,xlab ="Observações", ylab ="Resíduo studentizado")abline(h =c(-2, 2), col ="red", lty =2)which(abs(res_3) >2) # candidatos a outlier#> 9 #> 9
Figura 5: Resíduos studentizados com limites de referência ±2.
A observação 9 (\(t_9 \approx 4.3\)) destaca-se claramente.
O que fazer com o outlier?
Uma opção é reestimar o modelo sem a observação e comparar resultados.
# Comparação via estatística PRESSPRESS_completo <-sum(rstandard(mod, type ="pred")^2)PRESS_sem9 <-sum(rstandard(mod_sem9, type ="pred")^2)c(PRESS_completo = PRESS_completo, PRESS_sem9 = PRESS_sem9)#> PRESS_completo PRESS_sem9 #> 459.04 165.25
Remover a observação 9 reduz drasticamente\(\hat{\sigma}\) e o PRESS: ela é influente. A decisão de excluir deve, porém, considerar o contexto do problema.
Por que transformar?
Quando a análise de resíduos revela variância não constante ou não-linearidade, uma transformação de \(y\) e/ou dos regressores pode corrigir a inadequação e tornar o modelo linear adequado.
Dois objetivos principais:
Estabilizar a variância (Cap. 5.2);
Linearizar o modelo (Cap. 5.3).
Transformações estabilizadoras de variância (5.2)
Frequentemente a variância de \(y\) está ligada à sua média. Escolhemos a transformação conforme essa relação (Montgomery, Tabela 5.1):
p <-length(coef(mod)) # nº de parâmetros (b0, b1, b2)n <-nrow(entrega) # nº de observações2* p / n # limite de referência#> [1] 0.24h <-hatvalues(mod)which(h >2* p / n) # pontos de leverage#> 9 22 #> 9 22
As observações 9 e 22 têm valores de leverage acima de \(2p/n\).
Exemplo em R — distância de Cook
D <-cooks.distance(mod)plot(D, type ="h", lwd =2, main ="Distância de Cook",xlab ="Observação", ylab =expression(D[i]))abline(h =1, col ="red", lty =2)which(D >1) # observações influentes#> 9 #> 9
Figura 8: Distância de Cook por observação (limite de referência = 1).
A observação 9 (\(D_9 \approx 3.4\)) é claramente influente.
# DFFITS: influência de cada obs. na predição (corte 2*sqrt(p/n))dffits_val <-dffits(mod)2*sqrt(p / n)#> [1] 0.69282which(abs(dffits_val) >2*sqrt(p / n)) # observações destacadas#> 9 22 #> 9 22
Vários diagnósticos convergem para a observação 9 (e, em menor grau, a 22) como ponto de leverage/influência.
Modelos polinomiais
O modelo linear \(\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}\) serve para qualquer relação linear nos parâmetros — inclusive polinômios.
Segunda ordem em uma variável: \[
y = \beta_0 + \beta_1 x + \beta_2 x^2 + \varepsilon.
\]
Segunda ordem em duas variáveis (Eq. 7.7): \[
y = \beta_0 + \beta_1 x_1 + \beta_2 x_2
+ \beta_{11} x_1^2 + \beta_{22} x_2^2 + \beta_{12} x_1 x_2 + \varepsilon,
\]
com dois efeitos lineares, dois quadráticos e um de interação. Polinômios modelam bem respostas curvilíneas e superfícies de resposta.
Superfície de resposta (ilustração)
Figura 9: Superfície de resposta de um modelo de 2ª ordem em duas variáveis.
Exemplo em R — termos quadráticos
Usando o delivery time, ajustamos modelos polinomiais de 2ª ordem:
# Y = b0 + b1 x1 + b2 x1^2 + erromod_q1 <-lm(y ~ x1 +I(x1^2), data = entrega)summary(mod_q1)$coefficients#> Estimate Std. Error t value Pr(>|t|)#> (Intercept) 8.212645 2.309937 3.5554 0.0017721#> x1 1.135618 0.429848 2.6419 0.0148926#> I(x1^2) 0.034561 0.013785 2.5071 0.0200542# Y = b0 + b1 x2 + b2 x2^2 + erromod_q2 <-lm(y ~ x2 +I(x2^2), data = entrega)summary(mod_q2)$coefficients#> Estimate Std. Error t value Pr(>|t|)#> (Intercept) 1.1614e+01 2.7213e+00 4.26788 0.00031362#> x2 8.7999e-03 1.0450e-02 0.84214 0.40877202#> I(x2^2) 2.6646e-05 7.7082e-06 3.45679 0.00224496
Modelos encaixados.\(M_1\) está encaixado em \(M_2\) se todos os parâmetros de \(M_1\) estão contidos em \(M_2\), e \(M_2\) possui parâmetros adicionais. O modelo mod (só \(x_1, x_2\)) está encaixado em mod_int.
anova(mod, mod_int) # a interação melhora o ajuste?#> Analysis of Variance Table#> #> Model 1: y ~ x1 + x2#> Model 2: y ~ x1 + x2 + x1:x2#> Res.Df RSS Df Sum of Sq F Pr(>F) #> 1 22 234 #> 2 21 126 1 108 18 0.00037 ***#> ---#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod_nulo <-lm(y ~1, data = entrega)anova(mod_nulo, mod) # o modelo é significativo?#> Analysis of Variance Table#> #> Model 1: y ~ 1#> Model 2: y ~ x1 + x2#> Res.Df RSS Df Sum of Sq F Pr(>F) #> 1 24 5785 #> 2 22 234 2 5551 261 4.7e-16 ***#> ---#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Variáveis qualitativas na regressão
Muitas vezes precisamos usar variáveis qualitativas (categóricas) como preditoras: operador, turno, sexo, tipo de ferramenta, etc. Elas não têm escala natural de medida.
Atribuímos níveis por meio de variáveis indicadoras (ou dummy). Para dois níveis:
\[
x_2 = \begin{cases} 0 & \text{se a observação é do tipo A}\\ 1 & \text{se a observação é do tipo B} \end{cases}
\]
Em geral, uma variável qualitativa com \(a\) níveis é representada por \(a-1\) variáveis indicadoras.
Exemplo 8.1 — vida de ferramenta
Deseja-se relacionar a vida útil\(y\) (horas) de uma ferramenta de corte à velocidade do torno\(x_1\) (rpm) e ao tipo de ferramenta (A ou B).
Modelo de primeira ordem: \[
y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \varepsilon.
\]
Tipo A (\(x_2=0\)): \(\;y = \beta_0 + \beta_1 x_1 + \varepsilon\);
Tipo B (\(x_2=1\)): \(\;y = (\beta_0+\beta_2) + \beta_1 x_1 + \varepsilon\).
Duas retas paralelas (mesma inclinação \(\beta_1\)), com interceptos separados por \(\beta_2\).
ggplot(ferramenta, aes(x1, y, shape = tipo, color = tipo)) +geom_point(size =3) +labs(x ="Velocidade (rpm)", y ="Vida útil (horas)", color ="Tipo", shape ="Tipo")
Figura 10: Vida útil y vs. velocidade x1, por tipo de ferramenta.
Os pontos do tipo B situam-se acima dos do tipo A: a ferramenta B dura mais.
Ajuste do modelo com indicadora
mod_ind <-lm(y ~ x1 + tipo, data = ferramenta)summary(mod_ind)#> #> Call:#> lm(formula = y ~ x1 + tipo, data = ferramenta)#> #> Residuals:#> Min 1Q Median 3Q Max #> -5.553 -1.787 -0.002 1.839 4.984 #> #> Coefficients:#> Estimate Std. Error t value Pr(>|t|) #> (Intercept) 36.98560 3.51038 10.54 7.2e-09 ***#> x1 -0.02661 0.00452 -5.89 1.8e-05 ***#> tipoB 15.00425 1.35967 11.04 3.6e-09 ***#> ---#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#> #> Residual standard error: 3.04 on 17 degrees of freedom#> Multiple R-squared: 0.9, Adjusted R-squared: 0.889 #> F-statistic: 76.7 on 2 and 17 DF, p-value: 3.09e-09
Compare com a Tabela 8.2 de Montgomery: \[
\hat{y} = 36.986 - 0.027\,x_1 + 15.004\,x_2, \qquad R^2 = 0.9003.
\]
As duas retas ajustadas
b <-coef(mod_ind)ggplot(ferramenta, aes(x1, y, color = tipo, shape = tipo)) +geom_point(size =3) +geom_abline(intercept = b[1], slope = b[2], color ="#F8766D", linewidth =1) +geom_abline(intercept = b[1] + b[3], slope = b[2], color ="#00BFC4", linewidth =1) +labs(x ="Velocidade (rpm)", y ="Vida útil (horas)", color ="Tipo", shape ="Tipo")
Figura 11: Retas paralelas ajustadas para os tipos A e B.
O coeficiente \(\hat{\beta}_2 = 15.004\) é a diferença de vida útil entre B e A (para a mesma velocidade), e é altamente significativo (\(p < 10^{-8}\)).
Mais de dois níveis
Para três tipos A, B e C, usamos duas indicadoras \(x_2, x_3\):
Indicadoras também podem substituir um regressor quantitativo difícil de medir com precisão (ex.: faixas de renda) — Cap. 8.2.2.
Regressão linear simples — adsorção de corante
Artigo “Ecofriendly Dyeing of Silk with Extract of Yerba Mate” (Textile Res. J., 2017): efeito da concentração de corante (mg/L) sobre a adsorção (mg/g).
