\[E_i\sim N\left(0,\sigma^2\right)\iff Y_i\sim N\left(\beta_0+\beta_1 x_i,\sigma^2\right)\]
$
O que é um modelo estatístico?
Um modelo probabilístico que estabelece relações matemáticas entre algumas variáveis estatísticas
Uma forma de introduzir suposições sobre a observação de dados amostrais
O contexto
\(Y\) – uma variável resposta quantitativa
\({\bf x}=(x_1,\ldots,x_p)'\) – um qualquer número de preditores, covariáveis ou variáveis explicativas
Nota
Diferentes tipos de modelos (modelos de regressão, ANOVA, . . .) são necessários para incluir diferentes tipos de covariáveis
Algumas questões importantes
Qual o grau de associação entre essas covariáveis e a variável resposta?
Qual a precisão com que podemos predizer valores não observados de \(Y\), futuros ou passados?
Princípio fundamental da modelação
Variação observada em \(Y\)
=
Variação previsível + Variação aleatória
Forma funcional
\[Y = g(\mathbf{x}, \boldsymbol{\beta}) + E\]
onde \(E\) é chamado o erro aleatório
O modelo de regressão linear é caraterizado por:
\[g(\mathbf{x}, \boldsymbol{\beta})=\mathbf{x}' \boldsymbol{\beta}\]
e \(E\sim N(0, \sigma^2)\), não correlacionados \(\forall \mathbf{x}\)
Nota
O modelo é linear nos parâmetros!
Os dados
Um conjunto de pontos observáveis \((x_i,Y_i)\), \(i=1,\ldots,n\)
Modelo de Regressão Linear Simples (MRLS) \[Y_i=\beta_0+\beta_1x_i+E_i,\, i=1,\ldots,n\] \(\beta_0,\, \beta_1\) – parâmetros do MRLS
\(E_i\) – erro aleatório associado a \(Y_i=(Y\mid x=x_i)\)
Pressupostos usuais do MRLS
\(Var[E_i]=\sigma^2\) \(\iff Var[Y_i]=\sigma^2\)
\(E_i\)’s não correlacionadas \(\iff\) \(Y_i\)’s não correlacionadas
\(E_i\sim\) normal
\[E_i\sim N\left(0,\sigma^2\right)\iff Y_i\sim N\left(\beta_0+\beta_1 x_i,\sigma^2\right)\]
Para que servem os pressupostos?
. Tornam o modelo interpretável
\(E[Y\mid x]=\beta_0+\beta_1 x\) — curva de regressão
\(\beta_0=E[Y\mid x=0]\)
\(\beta_1=E[Y\mid x=x_0+1]-E[Y\mid x=x_0]\)
. Tornam simples a inferência estatística
Pelo método da máxima verossimilhança:
\[\left\lbrace\begin{array}{lcl} \hat {\beta }_1 &=& \frac{\sum_{i = 1}^n {x_i Y_i } - n\bar {x}\bar {Y}}{\sum_{i = 1}^n {x_i^2 } - n\bar {x}^2}\\ \hat{\beta }_0 &=& \bar{Y} - \hat{\beta }_1 \bar{x}\\ \hat{\sigma}_{MV}^2 & = & \frac{\sum_{i=1}^{n}{\left(Y_i-(\hat{\beta}_0+\hat{\beta}_1 x_i)\right)^2}}{n} \end{array}\right.\]
Nota
\(E[\hat{\beta}_i]=\beta_i\), \(i=1,2\)
\(E[\hat{\sigma}_{MV}^2]=\frac{n-2}{n}\sigma^2\implies \hat{\sigma}^2=\frac{n}{n-2}\hat{\sigma}_{MV}^2\) é um estimador centrado de \(\sigma^2\) que passaremos a usar
Mais ainda, temos:
\[\hat{\beta}_1 \sim N\left(\beta_1, \frac{\sigma^2}{Sxx}\right)\] \[\hat{\beta}_0 \sim N\left(\beta_0, \sigma^2 \left(\frac{1}{n} + \frac{\bar{x}^2}{Sxx} \right) \right)\]
em que \(Sxx = \sum_{i=1}^n (x_i-\bar{x})^2\)
e
\[\frac{(n-2)\hat{\sigma}^2}{\sigma^2}\sim \chi_{(n-2)}^2\]
\(\hat{\sigma}^2\) é independente de \(\hat{\beta}_0\) e \(\hat{\beta}_1\)
Os resultados anteriores conduzem a:
\[\frac{\hat{\beta}_1-\beta_1}{\hat{\sigma}\sqrt{\frac{1}{Sxx}}} = \frac{\hat{\beta}_1-\beta_1}{se(\hat{\beta}_1)} \sim t_{(n-2)}\]
e
\[\frac{\hat{\beta}_0-\beta_0}{\hat{\sigma}\sqrt{\frac{1}{n} + \frac{\bar{x}^2}{Sxx}}} = \frac{\hat{\beta}_0-\beta_0}{se(\hat{\beta}_0)} \sim t_{(n-2)}\]
que permitem construir intervalos de confiança ou testar hipóteses sobre \(\beta_1\) e \(\beta_0\).
Para que um modelo ajustado possa ser útil é necessário avaliar:
a adequação dos pressupostos
É possível encontrar algo que ponha em causa algum dos pressupostos?
a qualidade do ajustamento aos dados observados
O modelo permite explicar grande parte da variação observada em \(Y\)?
a capacidade preditiva do modelo
É possível obter predições precisas?
Os dados
Um conjunto de pontos observáveis \((\mathbf{x}_i,y_i)\), \(i=1,\ldots,n\) com \(\mathbf{x}_i=(x_{i,1}, \ldots,x_{i,p-1})\).
Vamos, por agora, considerar que todas as \(p-1\) covariáveis são VN.
Modelo de Regressão Linear Múltipla (MRLM) \[y_i = \beta_0 + \sum_{j=1}^{p-1}\beta_j x_{i,j} + E_i,\,i=1,\ldots,n\]
\(E_i\) — erro aleatório associado a \(y_i=(y\mid \mathbf{x}=\mathbf{x}_i)\) com \(E_i\sim N(0, \sigma^2)\), não correlacionadas.
O MRLM pode ser escrito em forma matricial:
\[\underset{n\times 1}{\mathbf{y}} = \underset{n\times p}{\mathbf{X}}\quad \underset{p\times 1}{\boldsymbol{\beta}}+\underset{n\times 1}{\mathbf{E}}\]
\[\begin{pmatrix}y_1\\ \vdots\\ y_n\end{pmatrix}=\begin{pmatrix} 1 & x_{1,1} & \cdots & x_{1,p-1}\\ \vdots & \vdots & \ddots & \vdots\\ 1 & x_{n,1} & \cdots & x_{n,p-1} \end{pmatrix} \begin{pmatrix}\beta_0\\ \vdots\\ \beta_{p-1}\end{pmatrix}+ \begin{pmatrix}E_1\\ \vdots\\ E_n\end{pmatrix}\]
\(\mathbf{X}\) — a matriz de delineamento
Notas
A distribuição normal multivariada
\(\mathbf{Y}\sim N_m(\boldsymbol{\mu}, \boldsymbol{\Sigma}),\,\mathbf{y}\in\mathbb{R}^m\)
\[f_{\mathbf{Y}}(\mathbf{y}) = (2\pi)^{-\frac{n}{2}}|\boldsymbol{\Sigma}|^{-\frac{1}{2}}\exp\left\{-\frac{1}{2}(\mathbf{y}-\boldsymbol{\mu})'\boldsymbol{\Sigma}^{-1}(\mathbf{y}-\boldsymbol{\mu})\right\}\]
\(E[\mathbf{y}]=\boldsymbol{\mu}\)
\(Cov[\mathbf{y}]=\boldsymbol{\Sigma}\)
Regressão de 1ª ordem
\[y_i = \beta_0 + \sum_{j=1}^{p-1}\beta_j x_{i,j} + E_i,\,i=1,\ldots,n\]
\[\begin{pmatrix}y_1\\ \vdots\\ y_n\end{pmatrix}=\begin{pmatrix} 1 & x_{1,1} & \cdots & x_{1,p-1}\\ \vdots & \vdots & \ddots & \vdots\\ 1 & x_{n,1} & \cdots & x_{n,p-1} \end{pmatrix} \begin{pmatrix}\beta_0\\ \vdots\\ \beta_{p-1}\end{pmatrix}+ \begin{pmatrix}E_1\\ \vdots\\ E_n\end{pmatrix}\]
Frequentemente é útil incluir termos de maior ordem:
| 1ª ordem | potências-(1) | \(x_1,\ldots, x_{p-1}\) |
| 2ª ordem | potências-(2) | \(x_1^2,\ldots, x_{p-1}^2\) |
| interações-(1,1) | \(x_1x_2,\ldots, x_1x_{p-1}, x_2x_3,\ldots, x_2x_{p-1},\ldots\) | |
| 3ª ordem | potências-(3) | \(x_1^3,\ldots, x_{p-1}^3\) |
| interações-(1,2) | \(x_1x_2^2,\ldots, x_1x_{p-1}^2, x_2x_3^2,\ldots, x_2x_{p-1}^2,\ldots\) | |
| interações-(2,1) | \(x_1^2x_2,\ldots, x_1^2x_{p-1}, x_2^2x_3,\ldots, x_2^2x_{p-1},\ldots\) | |
| . . . | . . . | . . . |
ou mesmo outras transformações das covariáveis.
Regressão polinomial
\[y_i = \sum_{j=0}^{p-1}\beta_j x_{i}^j + E_i,\,i=1,\ldots,n\]
\[\begin{pmatrix}y_1\\ \vdots\\ y_n\end{pmatrix}=\begin{pmatrix} 1 & x_{1} & x_{1}^2 & \cdots & x_1^{p-1}\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 1 & x_{n} & x_{n}^2 & \cdots & x_n^{p-1} \end{pmatrix} \begin{pmatrix}\beta_0\\ \vdots\\ \beta_{p-1}\end{pmatrix}+ \begin{pmatrix}E_1\\ \vdots\\ E_n\end{pmatrix}\]
Estimação pontual
Maximizando \[\log \mathcal{L}(\boldsymbol{\beta},\sigma^2)=-\frac{n}{2}\log \sigma^2-\frac{1}{2\sigma^2}(\mathbf{y}-\mathbf{X}\boldsymbol{\beta})'(\mathbf{y}-\mathbf{X}\boldsymbol{\beta}) + C\] obtêm-se os estimadores de MV:
\(\hat{\boldsymbol{\beta}}=(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{y}\)
\(\hat{\sigma}^2_{MV}=\frac{(\mathbf{y}-\mathbf{X}\hat{\boldsymbol{\beta}})'(\mathbf{y}-\mathbf{X}\hat{\boldsymbol{\beta}})}{n}\)
\(\hat{\mathbf{y}} = \hat{E}[\mathbf{y}\mid \mathbf{x}] = \mathbf{X}\hat{\boldsymbol{\beta}}\) — valores ajustados
Usando a expressão do estimador \(\hat{\boldsymbol{\beta}}\) temos:
\[\hat{\mathbf{y}}=\mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{y}= \mathbf{H}\mathbf{y}\]
\(\mathbf{H}_{n\times n}\) é chamada a matriz de projeção (ou matriz chapéu)
\(\mathbf{e} = \mathbf{y} - \hat{\mathbf{y}} = \mathbf{y} - \mathbf{H}\mathbf{y} = (\mathbf{I}_n-\mathbf{H})\mathbf{y}\) — resíduos
\(\mathbf{e}'\mathbf{e}=(\mathbf{y}-\mathbf{X}\hat{\boldsymbol{\beta}})'(\mathbf{y}-\mathbf{X}\hat{\boldsymbol{\beta}}) = SSE\) — soma dos quadrados dos resíduos
Sob os pressupostos do modelo:
\(\hat{\boldsymbol{\beta}} \sim N_p\left(\boldsymbol{\beta}, \sigma^2(\mathbf{X}'\mathbf{X})^{-1}\right) \implies \hat{\beta}_k \sim N\left(\beta_k, \sigma^2(\mathbf{X}'\mathbf{X})^{-1}_{k,k}\right)\)
\(\frac{SSE}{\sigma^2}=\frac{n\hat{\sigma}^2_{MV}}{\sigma^2}\sim \chi^2_{(n-p)}\)
\(\hat{\boldsymbol{\beta}}\) e \(\hat{\sigma}^2_{MV}\) são independentes
De 2. temos que \(E[SSE]=\sigma^2(n-p)\) o que mostra que \[E[\hat{\sigma}^2_{MV}]=E\left[\frac{SSE}{n}\right]=\frac{n-p}{n}\sigma^2\]
Assim, passaremos a usar o estimador \[\hat{\sigma}^2=\frac{SSE}{n-p}=MSE\]
Em particular, os resultados anteriores implicam
\[\frac{\hat{\boldsymbol{\beta}}_k-\boldsymbol{\beta}_k}{se(\hat{\boldsymbol{\beta}}_k)}\sim t_{(n-p)}\] com \(se(\hat{\boldsymbol{\beta}}_k)=\hat{\sigma}\sqrt{(\mathbf{X}'\mathbf{X})^{-1}_{kk}}\)
Estes resultados distribucionais podem ser usados para calcular intervalos de confiança ou testar hipóteses sobre cada um dos parâmetros de regressão individualmente.
Exemplo
tibble [234 × 11] (S3: tbl_df/tbl/data.frame)
$ manufacturer: chr [1:234] "audi" "audi" "audi" "audi" ...
$ model : chr [1:234] "a4" "a4" "a4" "a4" ...
$ displ : num [1:234] 1.8 1.8 2 2 2.8 2.8 3.1 1.8 1.8 2 ...
$ year : int [1:234] 1999 1999 2008 2008 1999 1999 2008 1999 1999 2008 ...
$ cyl : int [1:234] 4 4 4 4 6 6 6 4 4 4 ...
$ trans : chr [1:234] "auto(l5)" "manual(m5)" "manual(m6)" "auto(av)" ...
$ drv : chr [1:234] "f" "f" "f" "f" ...
$ cty : int [1:234] 18 21 20 21 16 18 18 18 16 20 ...
$ hwy : int [1:234] 29 29 31 30 26 26 27 26 25 28 ...
$ fl : chr [1:234] "p" "p" "p" "p" ...
$ class : chr [1:234] "compact" "compact" "compact" "compact" ...
Call:
lm(formula = hwy ~ displ + cyl, data = ggplot2::mpg)
Residuals:
Min 1Q Median 3Q Max
-7.510 -2.195 -0.205 1.902 14.922
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 38.216 1.048 36.46 <2e-16 ***
displ -1.960 0.519 -3.77 0.0002 ***
cyl -1.354 0.416 -3.25 0.0013 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.76 on 231 degrees of freedom
Multiple R-squared: 0.605, Adjusted R-squared: 0.601
F-statistic: 177 on 2 and 231 DF, p-value: <2e-16A significância de um modelo de regressão
Consideremos as hipóteses
\[H_0: \forall i>0 : \beta_i = 0\] contra \[H_1: \exists i>0 : \beta_i \neq 0\]
Nota
\(H_0\) significa que nenhuma dos termos contribui para explicar a variação de \(y\) sob o modelo ou
\[H_0: E[y\mid \mathbf{x}] = \beta_0\]
Sob \(H_0\), \(\hat{\beta}_0=\bar{y}\).
Pode mostrar-se que, para qualquer modelo de regressão:
\[\sum\left(y_i-\bar{y}\right)^2=\sum\left(y_i-\hat{y}_i\right)^2+\sum\left(\hat{y}_i-\bar{y}\right)^2\]
\[\iff SST=SSE+SSR\]
Nota
Cada modelo de regressão conduz a diferentes valores ajustados (\(\hat{\mathbf{y}}\)), e portanto, diferentes decomposições da variação observada na variável resposta.
Sob \(H_0\),
\[F=\frac{\frac{SSR}{p-1}}{\frac{SSE}{n-p}} = \frac{MSR}{MSE}\sim \mathcal{F}_{(p-1, n-p)}\]
e \(H_0\) pode ser rejeitada se \(F_o > F^{-1}_{\mathcal{F}_{(p-1, n-p)}}(1-\alpha)\).
Relacionado com o teste anterior, já conhecíamos o coeficiente de determinação:
\[R^2=\frac{SSR}{SST}=1-\frac{SSE}{SST}\]
Note-se que:
\(F=\frac{R^2}{1-R^2} \frac{n-p}{p-1}\)
\(SST\) é fixo e não depende do modelo
\(SSR\) não pode diminuir se se incluirem mais termos
\(R^2\) nunca diminui com a inclusão de termos!
Por causa desta inflação do \(R^2\), é preferível usar o:
Coeficiente de determinação ajustado
\[R^2_a=1-\frac{MSE}{MST}=1-\frac{(n-1) SSE}{(n-p) SST}\]
Notas
perde-se a interpretação mas \(R^2_a\) penaliza o aumento da complexidade do modelo
\(R^2_a\) pode ser negativo, ao contrário de \(R^2\)
Call:
lm(formula = hwy ~ displ + cyl, data = ggplot2::mpg)
Residuals:
Min 1Q Median 3Q Max
-7.510 -2.195 -0.205 1.902 14.922
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 38.216 1.048 36.46 <2e-16 ***
displ -1.960 0.519 -3.77 0.0002 ***
cyl -1.354 0.416 -3.25 0.0013 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.76 on 231 degrees of freedom
Multiple R-squared: 0.605, Adjusted R-squared: 0.601
F-statistic: 177 on 2 and 231 DF, p-value: <2e-16Comparação de modelos encaixados
Um modelo R (modelo reduzido) diz-se encaixado no modelo A (modelo alargado) se o conjunto de termos de R é um subconjunto ou um subespaço do conjunto de termos de A.
Para comparar modelos encaixados
\[H_0: \text{Modelo R}\ \ \text{contra}\ \ H_1: \text{Modelo A}\] podemos usar a estatística de teste
\[F = \frac{\frac{SSE(R) - SSE(A)}{gl_R-gl_A}}{\frac{SSE(A)}{gl_A}}\stackrel{H_0}{\sim}\mathcal{F}_{(gl_R-gl_A,\ gl_A)}\]
e \(H_0\) pode ser rejeitada se \(F^*_o > F^{-1}_{\mathcal{F}_{\left(gl_R-gl_A,\ gl_A\right)}}(1-\alpha)\)
Testes F sequenciais
Comparação de 2 modelos encaixados
Inferências para a resposta média
O estimador de \(E[\mathbf{y}\mid \mathbf{x}_0]=\mathbf{x}'_0\boldsymbol{\beta}\) é naturalmente \(\mathbf{x}'_0\hat{\boldsymbol{\beta}}\), e:
\[\frac{\mathbf{x}'_0\hat{\boldsymbol{\beta}}-\mathbf{x}'_0\boldsymbol{\beta}}{se(\mathbf{x}'_0\hat{\boldsymbol{\beta}})}\sim t_{(n-p)}\]
com \(se(\mathbf{x}'_0\hat{\boldsymbol{\beta}})=\hat{\sigma}\sqrt{\mathbf{x}'_0(\mathbf{X}'\mathbf{X})^{-1}\mathbf{x}_0}\)
Um intervalo de confiança com nivel de confiança \(\gamma\) é dado por:
\[\mathbf{x}'_0\hat{\boldsymbol{\beta}}\pm F_{t_{(n-p)}}^{-1}\left(\frac{1+\gamma}{2}\right)\times se(\mathbf{x}'_0\hat{\boldsymbol{\beta}})\]
Cuidado com as extrapolações!
Predição de uma nova observação
O objetivo é predizer uma nova observação sob a validade do modelo, isto é
\[y_0\sim N(\mathbf{x}'_0\boldsymbol{\beta}, \sigma^2)\]
independente dos dados observados.
O estimador do seu valor esperado, \(\hat{y}_0=\mathbf{x}'_0\hat{\boldsymbol{\beta}}\), é agora chamado um preditor.
Temos ainda
\[\frac{\hat{y}_0-y_0}{se(\hat{y}_0-y_0)}\sim t_{(n-p)}\]
com \(se(\hat{y}_0-y_0)=\hat{\sigma}\sqrt{1+\mathbf{x}'_0(\mathbf{X}'\mathbf{X})^{-1}\mathbf{x}_0}\)
Um intervalo de predição ao nível de confiança \(\gamma\) é dado por
\[\hat{y}_0\pm F_{t_{(n-p)}}^{-1}\left(\frac{1+\gamma}{2}\right)\times se(\hat{y}_0-y_0)\]
Cuidado com as extrapolações!
Exemplo – banda de confiança a 95%
Exemplo – bandas de confiança e predição a 95%
Covariáveis categorizadas
Quando incluímos covariáveis categorizadas é útil analisar a superfície de resposta:
\[E[Y\mid \mathbf{x}] = \mathbf{X} \boldsymbol{\beta}\]
Covariáveis binárias (2 níveis)
Suponhamos que \(p=3\) e que \(x_2\) é binária (0/1). A superfície de resposta do modelo de 1ª ordem sem transformações é:
\[\begin{align*} E[Y\mid \mathbf{x}] & = \beta_0 + \beta_1 x_1 + \beta_2 x_2 =\\[0.25cm] & = \begin{cases}\beta_0+\beta_1 x_1, & x_2=0\\[0.25cm](\beta_0 + \beta_2)+\beta_1 x_1, & x_2=1\end{cases} \end{align*}\]
\(\leadsto\) estamos a ajustar 2 superfícies de resposta paralelas!
Com \(q\) covariáveis binárias teríamos \(2^q\) superfícies de resposta paralelas.
Para remover o paralelismo é necessário incluir interações:
\[\begin{align*} E[Y\mid \mathbf{x}] & = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_1 x_2 =\\[0.25cm] & = \begin{cases}\beta_0+\beta_1 x_1, & x_2=0\\[0.25cm](\beta_0 + \beta_2)+(\beta_1 + \beta_3) x_1, & x_2=1\end{cases} \end{align*}\]
Mais de 2 níveis
No exemplo anterior, suponhamos agora que \(x_2\) é uma covariável categorizada ou fator com 3 níveis. Codificando esses níveis como 0, 1 ou 2 temos:
\[\begin{align*} E[Y\mid \mathbf{x}] & = \beta_0 + \beta_1 x_1 + \beta_2 x_2 =\\[0.25cm] & = \begin{cases}\beta_0+\beta_1 x_1, & x_2=0\\[0.25cm](\beta_0 + \beta_2)+\beta_1 x_1, & x_2=1\\[0.25cm](\beta_0 + 2\beta_2)+\beta_1 x_1, & x_2=2\end{cases} \end{align*}\]
Além de paralelas, as superfícies de resposta estão igualmente espaçadas e ordenadas. Isto é, em geral, indesejável.
Poderíamos associar uma variável binária a cada nível do fator, \(x_2\rightsquigarrow (x_{2,0}, x_{2,1}, x_{2,2})\), mas isso não conduziria a uma matriz de delineamento de posto completo.
Em vez disso usam-se apenas 2 variáveis indicadoras ou variáveis mudas:
| Nível | \(x_{2,1}\) | \(x_{2,2}\) |
|---|---|---|
| 0 | 0 | 0 |
| 1 | 1 | 0 |
| 2 | 0 | 1 |
Isto é chamada a codificação de nível de referência em que o nível codificado com \(\mathbf{x}=\mathbf{0}\) é o nível de referência.
\[\begin{align*} E[Y\mid \mathbf{x}] & = \beta_0 + \beta_1 x_1 + \beta_{2,1} x_{2,1} + \beta_{2,2} x_{2,2} =\\[0.25cm] & = \begin{cases}\beta_0+\beta_1 x_1, & x_2=0\\[0.25cm](\beta_0 + \beta_{2,1})+\beta_1 x_1, & x_2=1\\[0.25cm](\beta_0 + \beta_{2,2})+\beta_1 x_1, & x_2=2\end{cases} \end{align*}\]
Também aqui, a inclusão de interações pode remover o paralelismo.
Notas
Avaliação de modelos de regressão linear
Potenciais problemas
O preditor linear não é adequado
A estrutura do erro aleatório não é adequada
Variância não constante (heterocedasticidade)
Correlação nos erros
Não-normalidade
\(\leadsto\) inclusão de novos termos no modelo, transformações de \(Y\) ou das covariáveis, tentar outros modelos.
\(\leadsto\) requer mais investigação: há uma explicação? É um erro? Ou uma fraude? Considerar a remoção com cautela.
\(\leadsto\) Estandardizar as covariáveis ou usar outros modelos de regressão. Considerar a remoção de covariáveis com cautela.
Muitos dos meios de diagnóstico são baseados nos resíduos \(\mathbf{e}=\mathbf{y}-\hat{\mathbf{y}}\):
Métodos gráficos
resíduos em função de \(\hat{\mathbf{y}}\) e \(\mathbf{x}_i\)
Histogramas, box-plots, QQ-plots, . . .
Medidas numéricas, testes de hipóteses, . . .
Estratégia usual
Começar com alguns gráficos e usar testes apenas para clarificar indicações suspeitas ou pouco claras dadas pelos gráficos.
Resíduos
Sob os pressupostos do modelo, um gráfico dos resíduos deve ter o aspeto de uma nuvem simétrica de pontos em torno de zero com variação constante.
No entanto,
\[\mathbf{e}=(\mathbf{I} - \mathbf{H})\mathbf{y}\sim N_n(\mathbf{0}, \sigma^2(\mathbf{I} - \mathbf{H}))\]
\(\implies e_i\sim N(0, \sigma^2(1-h_{ii}))\)
\(\therefore\) os resíduos são v.a.’s correlacionadas com diferentes variâncias
Resíduos Studentizados
\[r_i=\frac{e_i}{se(e_i)}=\frac{e_i}{\sqrt{MSE (1-h_{ii})}},\ i=1,\ldots, n\]
são chamados os resíduos internamente Studentizados ou estandardizados.
Nota: rstandard() in R
Seja \[d_i=y_i-\hat{y}_{-i}\] em que \(\hat{y}_{-i}\) é o valor ajustado após a remoção da iª observação (deleted residual) e \(r_i^*\) o correspondente resíduo Studentizado.
Aparentemente, o cálculo destes resíduos deverá ser dispendioso!
Felizmente, tem-se \(d_i=\frac{e_i}{1-h_{ii}}\) e,
\[\begin{align*} r_i^* & = e_i\sqrt{\frac{n-p-1}{SSE (1-h_{ii})-e_i^2}}=\\[0.25cm] & = r_i \sqrt{\frac{n-p-1}{n-p- r_i^2}}\sim t_{(n-p-1)},\ i=1,\ldots, n \end{align*}\]
são chamados os resíduos externamente Studentizados.
Nota: rstudent() in R
Exemplo
Num gráfico de resíduos devemos olhar para:
a existência de algum padrão ou um desvio sistemático de zero (A, B1 and B2)
assimetria (B3)
valores extremos (C)
alguma dependência em alguma covariável relativamente aos pontos anteriores
Exemplos
Exemplos
Outros gráficos
Observações atípicas e alavancagem (leverage)
Observações influentes: alavancagem \(\times\) resíduo
Lembremos que \(Var[e_i]=\sigma^2 (1- h_{ii})\).
\(\leadsto\) um valor elevado de \(h_{ii}\) torna a \(Var[e_i]\) pequena e, consequentemente, a superfície de regressão estimada será mais fortemente “atraída” por \((y_i, \mathbf{x}_i)\).
Chamamos a \(h_{ii}\) a leverage da iª observação.
A distância de Cook é definida como
\[D_i=\frac{(\hat{\mathbf{y}} - \hat{\mathbf{y}}_{-i})'(\hat{\mathbf{y}} - \hat{\mathbf{y}}_{-i})}{p MSE}=\frac{1}{p} (r^*_i)^2\frac{h_{ii}}{1-h_{ii}}\]
e é usada como uma medida agregada da influência de cada observação.
Nota Usualmente os valores de \(D_i\) são comparados com os quantis de uma distribuição \(\mathcal{F}_{(p, n-p)}\).
leverage std.residual cooks.d
High leverage, Low residual 0.4096 -0.163 0.00926
Low leverage, High residual 0.0478 3.616 0.32819
High leverage, High residual 0.4096 -3.916 5.31837Regressão logística
Consideremos a observação de uma variável resposta binária:
\[Y_i \sim \text {Ber}(\theta_i),\, \theta_i=P(Y_i=1)\in ]0,1[,\]
independentes para \(i=1,\ldots,n\).
Neste caso, não podemos adotar a estratégia seguida na regressão linear, tomando \[E[Y_i\mid \mathbf{x}_i]=\theta_i=\mathbf{x}_i'\boldsymbol{\beta}\]
Porquê?
Em vez disso, podemos fazer \[g(\theta_i)=\mathbf{x}_i'\boldsymbol{\beta}\]
em que \(g: ]0,1[ \rightarrow \mathbb{R}\) é uma função de ligação.
Uma escolha possível é a função logito:
\[\mathop{\mathrm{\text{logit}}}(\theta)=\log \left(\frac{\theta}{1-\theta}\right)\]
A inversa é a função logística:
\[g^{-1}(x)=\frac{1}{1 + e^{-x}}=\frac{e^x}{e^x+1}\]
com \(x\in\mathbb{R}\).
Notas
A interpretação dos parâmetros dos parâmetros de regressão está agora relacionada com o logaritmo da chance de sucesso:
logaritmo da chance de sucesso para \(\mathbf{x}=\mathbf{0}\)
variação do logaritmo da chance de sucesso quando \(x_i\) aumenta uma unidade e as restantes covariáveis permanecem constantes
Toda a restante inferência baseia-se em aproximações que requerem amostras grandes como, por exemplo:
\[\frac{\hat{\beta}_i-\beta_i}{se(\hat{\beta}_i)}\stackrel{a}{\sim} N(0,1)\]
em que \(se(\hat{\beta}_i)\) é calculado a partir da matriz Hessiana da função de verossimilhança.
Generalização (observações repetidas)
\[Y_i \sim \text {binomial}(n_i, \theta_i),\]
independentes para \(i=1,\ldots,k\).
\(E\left[\frac{Y_i}{n_i}\right] = \theta_i\) e \(\mathop{\mathrm{\text{logit}}}(\theta_i)=\mathbf{x}_i'\boldsymbol{\beta}\)
Exemplo
Na tabela seguinte apresenta-se o número de mortes, \(y_i\), em 8 grupos com \(n_i\) besouros expostos a uma dose \(x_i\) de \(CS_2\).
| \(i\) | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 |
|---|---|---|---|---|---|---|---|---|
| \(y_i\) | 6 | 13 | 18 | 28 | 52 | 53 | 61 | 60 |
| \(n_i\) | 59 | 60 | 62 | 56 | 63 | 59 | 62 | 60 |
| \(x_i\) | 1.6907 | 1.7242 | 1.7552 | 1.7842 | 1.8113 | 1.8369 | 1.8610 | 1.8839 |
data.frame(
y = c(6, 13, 18, 28, 52, 53, 61, 60),
n = c(59, 60, 62, 56, 63, 59, 62, 60),
x = c(1.6907, 1.7242, 1.7552, 1.7842, 1.8113, 1.8369, 1.8610, 1.8839)
) -> besouros
modelo <- glm(cbind(y, n - y) ~ x, family = "binomial", data = besouros)
estimativas <- coef(modelo)
estimativas(Intercept) x
-60.7 34.3
Call:
glm(formula = cbind(y, n - y) ~ x, family = "binomial", data = besouros)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -60.72 5.18 -11.7 <2e-16 ***
x 34.27 2.91 11.8 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 284.202 on 7 degrees of freedom
Residual deviance: 11.232 on 6 degrees of freedom
AIC: 41.43
Number of Fisher Scoring iterations: 4
\(\leadsto \widehat{\text{DL}}_{50}\approx 1.77\) e \(\widehat{\text{DL}}_{90}\approx 1.83\)
Qualidade do ajustamento
A desviância de um modelo \(M\) é dada por
\[-2\log \frac{\mathcal{L}\left(\hat{\boldsymbol{\beta}}_M\right)}{\mathcal{L}\left(\hat{\boldsymbol{\beta}}_S\right)}= 2\left(\log \mathcal{L}\left(\hat{\boldsymbol{\beta}}_S\right)- \log \mathcal{L}\left(\hat{\boldsymbol{\beta}}_M\right)\right)\]
em que \(S\) representa o modelo saturado (um parâmetro por observação).
No R temos:
Desviância nula: \(2\left(\log \mathcal{L}\left(\hat{\boldsymbol{\beta}}_S\right)- \log \mathcal{L}\left(\hat{\boldsymbol{\beta}}_0\right)\right)\)
(\(gl=n-1\))
Desviância residual: \(2\left(\log \mathcal{L}\left(\hat{\boldsymbol{\beta}}_S\right)- \log \mathcal{L}\left(\hat{\boldsymbol{\beta}}_M\right)\right)\)
(\(gl=n-p\))
A diferença é encarada como uma medida da capacidade explicativa de um modelo M relativamente ao modelo nulo.
A medida
\[\text{AIC} = 2p - 2\log \mathcal{L}\left(\hat{\boldsymbol{\beta}}_M\right)\]
é conhecida como o critério de informação de Akaike.
Tal como com a desviância, a diminuição desta medida é indicadora de uma melhor qualidade de ajustamento.
Ao contrário da desviância, o AIC penaliza o aumento de complexidade dos modelos.
Comparação de modelos encaixados
Para testar a hipótese
\(H_0: \text{Modelo reduzido (R)}\) contra \(H_1: \text{Modelo alargado (A)}\)
pode-se usar a estatística de teste
\[G^2 = -2\log \frac{\mathcal{L}\left(\hat{\boldsymbol{\beta}}_R\right)}{\mathcal{L}\left(\hat{\boldsymbol{\beta}}_A\right)}\underset{H_0}{\stackrel{a}{\sim}}\chi^2_{(p_A - p_R)}\]
Rejeita-se \(H_0\) se \(G^2_o > F^{-1}_{\chi^2_{\left(p_A - p_R\right)}}(1-\alpha)\).
Regressão de Poisson
A ideia anterior pode ser aplicada a outras distribuições. Por exemplo, para dados resultantes de contagens, a distribuição de Poisson é frequentemente usada.
\[Y_i \sim \text {Poi}(\lambda_i),\, \lambda_i\in \mathbb{R}^+,\]
independentes para \(i=1,\ldots,n\).
Neste caso, a modelação mais usual é \(\log(\lambda_i)=\mathbf{x}_i'\boldsymbol{\beta}\).
Os modelos de regressão fazem parte da família dos modelos lineares (lm no R).
Os modelos de regressão logística e de Poisson fazem parte da família dos modelos lineares generalizados (glm no R).