Conjunto de dados “cars”

O R possui pré-instalado diversos conjuntos de dados (datasets) para utilização direta pelo usuário (sem necessidade de importação). Entre com o seguinte comando para visualizar as opções disponíveis: 

utils::data()

Para o nosso exemplo de regressão linear simples em R usaremos o conjunto de dados “cars”. Ele é composto por 50 observações da velocidade de diferentes carros da década de 20 e da respectiva distância alcançada pelo carros durante o experimento. 

str(cars)  # Retorna a estrutura do conjunto de dados.
## 'data.frame':    50 obs. of  2 variables:
##  $ speed: num  4 4 7 7 8 9 10 10 10 11 ...
##  $ dist : num  2 10 4 22 16 10 18 26 34 17 ...

Antes de ajustar qualquer modelo, vamos observar o comportamento dos dados com base em algumas técnicas exploratórias gráficas e parâmetros estatísticos. Com a função stat.desc do pacote pastecs podemos verificar estimativas para a média, mediana, range, variância, dentre outros. 

library(pastecs)
stat.desc(cars)
##                    speed         dist
## nbr.val       50.0000000   50.0000000
## nbr.null       0.0000000    0.0000000
## nbr.na         0.0000000    0.0000000
## min            4.0000000    2.0000000
## max           25.0000000  120.0000000
## range         21.0000000  118.0000000
## sum          770.0000000 2149.0000000
## median        15.0000000   36.0000000
## mean          15.4000000   42.9800000
## SE.mean        0.7477858    3.6443403
## CI.mean.0.95   1.5027319    7.3235761
## var           27.9591837  664.0608163
## std.dev        5.2876444   25.7693775
## coef.var       0.3433535    0.5995667

Neste momento, os valores tabelados acima não nos agregaram mais informações do que visualizando, propriamente dito, os dados. Se queremos um modelo, é esperado que a variável resposta e as preditores possuam uma boa correlação entre si. Para isso, plotemos os valores de velocidade contra as distâncias: 

par(mar = c(5,4,2,2))

plot(x = cars$speed, 
     y = cars$dist,
     xlab = "Velocidade (mph)", 
     ylab = "distância (ft)",
     pch = 16)

O gráfico acima apresenta uma boa relação entre a velocidade do carro e a distância percorrida por eles no experimento. O coeficiente de correlação de pearson concorda com o gráfico, obtendo um valor igual a 0.8068949.

No código acima usamos o argumento pch para a plotagem dos dados. Este argumento seleciona o formato do simbolo plotado. No gráfico abaixo, tem-se o simbolo associado a cada valor de pch:

Para a construção dos histogramas abaixo usamos o argumento mfrow dentro da função par. Esta última controla todos os argumentos gráficos, enquanto mfrow define a organização de múltiplos gráficos num mesmo painel semelhante a uma matriz. Neste caso nosso painel tem uma linha e duas colunas. 

par(mfrow = c(1, 2)) 

# histograma das velocidades
hist(x = cars$speed,
     col = "gray",
     xlab = "Velocidade (mph)",
     ylab = "Frequência",
     main = "")

# histograma das distâncias
hist(x = cars$dist,
     col = "lightblue",
     xlab = "Distância (ft)",
     ylab = "Frequência",
     main = "")

Com os boxplots abaixo, observamos um possível ponto aberrante (outlier) em nossos dados. Em geral, estamos acostumados com a ideia de um outlier em modelo, ou seja, alguma observação (dados) que o nosso modelo não tem boa aderência e por isso produz um alto resíduo. Não é este o caso. 

par(mar = c(4,4,2,2))

# boxplot dos dados de velocidade e distância
boxplot(cars$speed, cars$dist,
        names = c("Velocidade", "Distância"),
        col = c("lightgray", "lightblue"))

Este possível outlier corresponde à observação 49 no conjunto de dados, na qual um carro com velocidade igual a 24 mph alcança 120 ft. Devemos remover este ponto de nossos dados? Para isso devemos perguntar: Esta observação representa nosso conjunto de dados? Estamos tratando de um mesmo tipo de carro? Ou será que ele alcança uma distância menor porque é na realidade um carro de corrida e não um carro urbano? Infelizmente, não temos essa informação. Então, de forma conservadora, vamos manter essa observação em nosso conjunto de dados.


Construção do modelo linear

Para construção de nosso modelo \(distancia = \alpha.velocidade + \beta\) vamos utilizar a função nativa do R lm (do inglês, linear model). Nela escrevemos a equação anterior com a seguinte sintaxe: Y ~ X. Ou seja, Y é nossa variável resposta (distância) e X nossa preditora (velocidade). o simbolo “~” equivale a dizer “é função de”. A função lm, então estima o coeficiente da variável preditora e a constante da função por padrão. 

modelo.simples = lm(cars$dist ~ cars$speed)
summary(modelo.simples)
## 
## Call:
## lm(formula = cars$dist ~ cars$speed)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -29.069  -9.525  -2.272   9.215  43.201 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -17.5791     6.7584  -2.601   0.0123 *  
## cars$speed    3.9324     0.4155   9.464 1.49e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.38 on 48 degrees of freedom
## Multiple R-squared:  0.6511, Adjusted R-squared:  0.6438 
## F-statistic: 89.57 on 1 and 48 DF,  p-value: 1.49e-12

A função retorna no console (acima) diversas informações sobre o modelo e sua adequação.

Na linha Call temos a função ajustada pela função lm. Na linha Residual são apresentadas algumas informações sobre os resíduos do modelo (valor mínimo, máximo, mediana, 1° e 3° Quartis).

Na linha Coeficients é apresentada uma tabela contendo o valor dos coeficientes (Estimate), o erro padrão do coeficiente (Std. Error), o valor da estatística t (t value) e finalmente o p-valor do teste t para significância do coeficiente (Pr(>|t|)).

Ao lado direito do p-valor do teste dos coeficientes, ainda é apresentado por meio de símbolos (explicados na linha abaixo), para qual nível de significância os coeficientes são significativos.

Por fim, é realizado um teste F para adequação do modelo, apresentado para isso o número de graus de liberdade envolvidos, o erro padrão residual, os valores para o R2 e o R2 ajustado, o consequente valor da estatística F e seu respectivo p-valor. Caso deseje-se obter a tabela de ANOVA separadamente, basta entrar com a função: 

anova(modelo.simples)
## Analysis of Variance Table
## 
## Response: cars$dist
##            Df Sum Sq Mean Sq F value   Pr(>F)    
## cars$speed  1  21186 21185.5  89.567 1.49e-12 ***
## Residuals  48  11354   236.5                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Com base nas informações acima, podemos afirmar que o modelo é estatisticamente significante (p-valor do teste F aproximadamente nulo), da mesma forma que com os coeficientes de regressão (p-valor do teste T). O valor do R2 indica um aderência moderada aos dados (~ 0.6).

Contudo, o que garante que as conclusões anteriores sejam verdadeiras? A estimativa dos p-valores, do R2 e dos coeficientes é em sua essência uma função matemática, independente dos valores utilizados, um resultado será obtido. Então, o que garante os resultados tenham valor estatístico?

Resposta: O cumprimento das considerações do modelo de mínimos quadrados:

  • Resíduos com distribuição normal
  • Homoscedasticidade dos resíduos (variância constante)
  • Aleatoriedade dos resíduos frente ao valor predito e às variáveis preditoras

Estas condições vem do desenvolvimento do métodos dos mínimos quadrados a partir do método da máxima verossimilhança e devem ser avaliados a posteriori da construção do modelo.

Validação do modelo linear

Inicialmente, avaliemos a hipótese de homoscedasticidade dos resíduos de forma gráfica, plotando os valores preditos pelo modelo por aqueles observados (distância). 

par(mar = c(5,4,2,2))

plot(x = cars$dist, 
     y = modelo.simples$fitted.values,
     xlim = c(0, 120), 
     ylim = c(0, 120),
     pch = 16,
     xlab = "Valores observados",
     ylab = "Valores preditos")

abline(a = 0, b = 1, lty = 2, col = "red")

Com base no gráfico anterior, verifica-se que para valores acima de 80 para a distância, a predição do modelo começa a ser prejudicada, estimando valores menores do que o observado. Observe também que o modelo pode estimar valores negativos para a distância quando valores baixos de velocidade são utilizados, o que compromete o domínio de aplicação do modelo.

A mesma conclusão pode ser obtida avaliando o gráfico a direta no painel a seguir. À esquerda, no gráfico de resíduos versus a variável preditora, por outro lado, não se encontra indícios visuais da não-aleatoriedade dos resíduos. 

par(mfrow = c(1, 2), mar = c(5,4,2,2))

plot(x = cars$speed, 
     y = modelo.simples$residuals,
     pch = 16,
     col = "red",
     xlab = "Velocidade (mph)", ylab = "Resíduos")

plot(x = cars$dist, 
     y = modelo.simples$residuals,
     pch = 16,
     col = "lightblue",
     xlab = "Distância (ft)", ylab = "Resíduos")

A distribuição de probabilidade dos resíduos pode ser verificada visualmente a partir de um histograma dos resíduos. Que conforme a figura abaixo indica uma cauda não simétrica a direita da distribuição, o que sugere a não normalidade dos resíduos. 

par(mar = c(5,4,2,2))
hist(x = modelo.simples$residuals,
     xlab = "Resíduos",
     ylab = "Densidade de Probabilidade",
     main = "",
     col = "lightgreen",
     probability = TRUE)

lines(density(modelo.simples$residuals))

Aplicando o teste de Shapiro-Wilk para os resíduos (shapiro.test), concluímos que para um nível de significância de 5%, nossos dados não vêm de uma distribuição normal. 

shapiro.test(modelo.simples$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  modelo.simples$residuals
## W = 0.94509, p-value = 0.02152

Conclusões

A validação do modelo demonstrou até certo ponto que as considerações necessárias para o método dos mínimos quadrados não foram atendidos. Aplicabilidade e significância estatística, são, contudo, aspectos diferentes.

A validação não é finalizada apenas com estas análises, técnicas como a chamada “crossvalidation” ainda podem ser aplicadas.

Outra questão interessante a ser discutida é a robustez dos teste estatísticos (F e t) contra a fuga da normalidade ou aleatoriedade dos resíduos.