Como saber se a função de regressão é adequada para explicar a relação entre as variáveis observadas?

• Gráficos

• Testes

Caso a função de regressão não seja adequada, o que fazer?

• Transformação de variáveis.

• Outro tipo de modelo.

Gráficos de diagnósticos para \(Y\) não são úteis em análise de regressão, pois a variável resposta é uma função da preditora.

Os gráficos de diagnósticos são feitos para os resíduos.

• Resíduos:

\[e_i=Y_i-\hat{Y}_i\]

• Erros:

\[\varepsilon_i=Y_i-E(Y_i)\]

Se assumimos \(\varepsilon_i\underset{iid}{\sim}\mathcal{N}(0,\sigma^2)\), esperamos que \(e_i\)'s, os resíduos observados, reflitam as propriedades assumidas para \(\varepsilon_i\)'s.

• Média dos resíduos:

\[\bar{e}=\frac{\sum_{i=1}^ne_i}{n}=0\]

Como \(\bar{e}\) é sempre nula, não fornece informação sobre \(E(\varepsilon_i)\).

• Variância dos resíduos:

\[s^2=\frac{\sum_{i=1}^n(e_i-\bar{e})^2}{n-2}=\frac{SQErro}{n-2}=QMErro\]

Se o modelo é apropriado, \(s^2\) é um estimador não viesado para \(\sigma^2\).

• Independência: \(e_i\)'s não são variáveis aleatórias independentes, pois dependem de \(\hat{Y}_i\), que por sua vez dependem das mesmas estimativas para os parâmetros \(\beta_0\) e \(\beta_1\).

Quanto maior o número de observações com relação ao número de parâmetros estimados, menor é a dependência dos resíduos (pode ser ignorada).

Para detectar outliers, muitas vezes é mais útil observar resíduos padronizados.

Resíduo semi-studentizado:

\[e_i^*=\frac{e_i}{\sqrt{QME}}\]

Podemos utilizar análise de resíduos para verificar se:

1. a função de regressão não é linear

2. erros não possuem variância constante (heterocedasticidade)

3. erros não são independentes

4. o modelo ajustado é adequado, com exceção de algumas observações outliers

5. erros não seguem distribuição normal

6. uma ou mais variáveis preditoras foram omitidas no modelo

• Resíduos x Variável preditora

• Resíduos absolutos x Variável preditora

• Resíduos x Valores ajustados (\(\hat{Y}\))

• Resíduos ao longo do tempo/em sequência

• Resíduos X Variável preditora não incluída no modelo

• Boxplot dos resíduos

• Gráfico de probabilidade normal dos resíduos

Calcular coeficiente de correlação entre \(e_i\) e seus valores esperados, segundo a distribuição normal

Quanto maior a correlação, maior o indício de normalidade dos resíduos.

Detalhes sobre o teste: Looney & Gulledge (1985) - Use of the Correlation Coefficient with Normal Probability Plots

Calcular a correlação entre os resíduos ordenados e

\[r=\{r_1,r_2,\ldots,r_k\ldots,r_n\}\]

em que \(r_k\) é o valor esperado do \(k\)-ésimo resíduos, segundo a suposição de normalidades.

\[r_k\approx \sqrt{QME}\left[z\left(\frac{k-0.375}{n+0.25} \right) \right]\]

\(z(A)\) é o \(A\times 100\) percentil da distribuição normal padrão.

Obs: Outro teste conhecido: Teste de Shapiro-Wilks.

• Não depende da suposição de normalidade

• Aplicável no caso de regressão linear simples quando a variância do erro aumenta ou decresce com a variável preditora \(X\).

• Tamanho amostral grande o suficiente para que possamos assumir que os resíduos são independentes.

• Idéia parecida com um teste-t para duas amostras.

• \(X\) é dividida em \(X_1\) (valores mais baixos de \(X\)) e \(X_2\) (valores mais altos de \(X\)). Isto é, temos dois grupos.

• \(e_{i1}\) é o \(i\)-ésimo resíduo para o grupo 1 e \(e_{i2}\) é o \(i\)-ésimo resíduo para o grupo 2.

• \(n_j\) é o número de observações no grupo \(j\)

• \(n=n_1+n_2\).

• \(\tilde{e}_j\) é a mediana dos resíduos do grupo \(j\).

• \(d_{i1}=| e_{i1}-\tilde{e}_1|\quad\quad d_{i2}=| e_{i2}-\tilde{e}_2|\)

Estatística do teste:

\[t^*_{BF}=\frac{\bar{d}_1-\bar{d}_2}{s\sqrt{\frac{1}{n_1}+\frac{1}{n_2}}}\]

\[s^2=\frac{\sum(d_{i1}-\bar{d}_1)^2+\sum(d_{i2}-\bar{d}_2)^2}{n-2}\]

Se a variância é constante e \(n_1\) e \(n_2\) grandes o suficiente, então \(t^*_{BF}\) segue distribuição \(t\) com \(n-2\) graus de liberdade.

Valores altos de \(t^*_{BF}\) indicam heterocedasticidade.

Assume que \(\varepsilon_i\sim\mathcal{N}(0,\sigma_i^2)\).

O seguinte modelo é assumido:

\[\log_e\sigma_i^2=\gamma_0+\gamma_1X_i\]

Se \(\gamma_1=0\) temos homocedasticidade.

\(H_0\): \(\gamma_1=0\) vs \(H_1\): \(\gamma_1\neq0\).

Estatística do teste:

\[X^2_{BP}=\frac{\frac{SQRq}{2}}{\left(\frac{SQE}{n}\right)^2}\]

em que \(SQRq\) é a soma de quadrados da regressão de \(e_i^2\) em \(X\) e \(SQE\) é a soma de quadrados do erro da regressão de \(Y\) em \(X\).

Sob \(H_0\), \(X^2_{BP}\sim\chi^2_1\).

Diretamente pelo R:

library('car')
ncvTest(fit)

## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 0.8209192    Df = 1     p = 0.3649116

Queremos testar se a função de regressão proposta é adequada para os dados.

Suposições: as observações de \(Y|X\) são

• independentes
• seguem distribuição normal
• variância constante \(\sigma^2\)

É preciso que tenhamos mais de uma observação para cada nível de \(X\) (replicações).

Suponha que tenhamos \(m\) valores distintos para \(X\)

Para cada valor diferente de \(X\), \(x_j\), temos \(n_j\) observações.

• \(Y_{11},Y_{12},\ldots, Y_{1n_1}\) são \(n_1\) repetições para \(X=x_1\).
• \(Y_{21},Y_{22},\ldots, Y_{2n_2}\) são \(n_2\) repetições para \(X=x_2\).
• \(Y_{m1},Y_{m2},\ldots, Y_{mn_m}\) são \(n_m\) repetições para \(X=x_m\).

\[n=\sum_{j=1}^m\sum_{u=1}^{n_j}1=\sum_{j=1}^mn_j\]

Com replicações para cada \(X=x_k\), temos um estimador para \(Var(Y|X=x_k)=\sigma^2\) para todo \(k\).

Por exemplo, quando \(X=x_1\), podemos usar \(Y_{11},Y_{12},\ldots, Y_{1n_1}\) e o seguinte estimador para \(Var(Y|X=x_1)\):

\[\frac{1}{n_1-1}\sum_{j=1}^{n_1}(Y_{1j}-\bar{Y}_1)^2\]

Combinando todas as observações (não apenas as replicações para um dado nível de \(X\)), temos o seguinte estimador para \(\sigma^2\):

\[s_p^2=\frac{\sum_{j=1}^m\sum_{u=1}^{n_j}(Y_{ju}-\bar{Y}_j)^2}{\sum_{j=1}^m(n_j-1)}\]

\[SQE_p=\sum_{j=1}^m\sum_{u=1}^{n_j}=(Y_{ju}-\bar{Y}_j)^2\]

com graus de liberdade \(n_p=\sum_{j=1}^m(n_j-1)=\sum_{j=1}^mn_j-m\).

\[QME_{puro}=s_p^2=\frac{\sum_{j=1}^m\sum_{u=1}^{n_j}(Y_{ju}-\bar{Y}_j)^2}{\sum_{j=1}^mn_j-m}\]

Relembrando:

\[SQE=\sum_{j=1}^m\sum_{u=1}^{n_j}(Y_{ju}-\hat{Y}_j)^2\]

Mostre\(^1\) que a \(SQE\) pode ser particionada da seguinte forma:

\[\sum_{j=1}^m\sum_{u=1}^{n_j}(Y_{ju}-\hat{Y}_j)^2 = \sum_{j=1}^m\sum_{u=1}^{n_j}(Y_{ju}-\bar{Y}_j)^2 + \sum_{j=1}^mn_j(\hat{Y}_j-\bar{Y}_j)^2\]

\[\underbrace{\sum_{j=1}^m\sum_{u=1}^{n_j}(Y_{ju}-\hat{Y}_j)^2}_{SQE} = \underbrace{\sum_{j=1}^m\sum_{u=1}^{n_j}(Y_{ju}-\bar{Y}_j)^2}_{SQE_p} + \underbrace{\sum_{j=1}^mn_j(\hat{Y}_j-\bar{Y}_j)^2}_{SQF_a}\]

\(SQE\) é a soma de quadrados dos resíduos, representando a variação em torno da reta.

\(SQE_p\) é a soma de quadrados do erro puro, representando a variação de \(Y\), para \(X\) fixo, independente do modelo (pois utilizamos as replicações).

\(SQF_a\) é a soma de quadrados da falta de ajuste, representando a falta de ajuste do modelo.

\(SQE\) é o numerador de \(s^2\) e \(SQE_p\) é o numerador de \(s^2_p\).

Valores relativamente altos de \(SQF_a\) indicam discrepância entre esses dois estimadores e, consequentemente, indicam evidência contra a adequação do modelo de regressão linear simples.

Quanto mais alto o valor de \(SQF_a\), pior é o ajuste do modelo, pois, se o modelo fosse perfeito, teríamos que \(SQF_a=0\).

\(H_0\): o MRLS é adequado

\(H_1\): o MRLS não é adequado

Estatística do teste:

\[F^*=\frac{SQF_a/(m-2)}{SQE_p/(n-m)}\overset{\mbox{sob }H_0}{\sim} F_{m-2,n-m}\]

11 agências de um certo banco ofereceram brindes para novos clientes.

Um depósito inicial mínimo era necessário para qualificar para o brinde.

Cada agência especificou o valor do depósito mínimo.

Interesse: valor do depósito mínimo inicial e número de contas que foram abertas na agência.

Modelo completo:

\[Y_{ij}=\mu_{j}+\varepsilon_{ij}\]

Temos que \(E(Y_{ij})=\mu_j\), isto é, temos esperança diferente para cada \(X_j\).

13 ligas de cobre e níquel, mas quantidades diferentes de ferro.

As ligas foram submersas em água do mar por 60 dias e a perda de peso devido à corrosão foi anotada em miligramas por decímetro quadrado por dia.

## 
## Call:
## lm(formula = loss ~ Fe, data = corrosion)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7980 -1.9464  0.2971  0.9924  5.7429 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  129.787      1.403   92.52  < 2e-16 ***
## Fe           -24.020      1.280  -18.77 1.06e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.058 on 11 degrees of freedom
## Multiple R-squared:  0.9697, Adjusted R-squared:  0.967 
## F-statistic: 352.3 on 1 and 11 DF,  p-value: 1.055e-09

anova(modeloreduzido,modelocompleto)

## Analysis of Variance Table
## 
## Model 1: loss ~ Fe
## Model 2: loss ~ factor(Fe) - 1
##   Res.Df     RSS Df Sum of Sq      F   Pr(>F)   
## 1     11 102.850                                
## 2      6  11.782  5    91.069 9.2756 0.008623 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Estatística do teste:

\[F^*=\frac{SQF_a/(m-2)}{SQE_p/(n-m)}=\frac{91.069/5}{11.782/6}=9.28\]

Se nos testes e/ou gráficos de diagnóstico do modelo determinarmos que o mesmo não é adequado, o que fazer?

• Abandonar o modelo e procurar outro mais adequado: regressão logística, regressão não-paramétrica, regressão robusta…

• Aplicar alguma transformação nos dados, de maneira que a regressão linear simples fique adequada.