MLGExemplo06

Exemplo

Um joalheiro precifica os diamantes com base na qualidade, com valores de 0 a 8, sendo 8 perfeito e 0 contendo numerosas imperfeições e cor, com valores de 1 a 10, sendo 10 branco puro e 1 amarelo. Com base no preço por quilate, em centenas de dólares, dos 11 diamantes a seguir, pesando entre 1.0 e 1.5 quilates.

Determinar a relação entre qualidade (quality), cor (color) e preço (price).

Observe que o preço do diamante é a resposta; isto deve-se à que a cor e qualidade são características naturais.

Color = c(7,3,5,8,9,5,4,2,8,6,9)                 # color do diamante, entre 1 e 10
Quality = c(5,7,8,1,3,4,0,6,7,4,2)               # qualidade do diamante, entre 0 e 8
Price = c(65,38,51,38,55,43,25,33,71,51,49)      # preço por quilate, em centenas de dólares

dados = data.frame(Color,Quality,Price)

Primeiro uma análise descritiva dos dados.

library(ggplot2)
g1 = ggplot(data = dados, aes(x = Quality, y = Price)) + geom_point() 
g2 = ggplot(data = dados, aes(x = Color, y = Price)) + geom_point()

library(patchwork)
g1 + g2

Não visualizamos muito nas figuras acima. Fazeremos transformações nas variáveis explicativas.

g3 = ggplot(data = dados, aes(x = exp(Quality), y = Price)) + geom_point() 
g4 = ggplot(data = dados, aes(x = I(Quality*Color), y = Price)) + geom_point()

(g1 + g2) / (g3 + g4)

Modelo

Observamos que a interação Quality X Color induz uma excelente explicação do Price (preço). Nosso objetivo é mostrar diferentes métodos de diagnóstico; por isso não incliremos numa primeira análise a descoberta anteriormente mencionada.

modelo0 = lm(Price ~ Quality + Color, data = dados)
summary(modelo0)

## 
## Call:
## lm(formula = Price ~ Quality + Color, data = dados)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -6.672 -4.536 -1.093  3.722 10.190 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.7514     6.9602   0.252 0.807670    
## Quality       3.7584     0.7565   4.968 0.001096 ** 
## Color         4.8953     0.8202   5.968 0.000335 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.888 on 8 degrees of freedom
## Multiple R-squared:  0.8507, Adjusted R-squared:  0.8134 
## F-statistic: 22.79 on 2 and 8 DF,  p-value: 0.0004969

Resíduos

library(car)
influenceIndexPlot(modelo0)

Lembremos que, num modelo linear, Uma medida da influência de cada observação nas estimativas dos parâmetros da regressão é a distância de Cook, dada por \[\begin{equation*} C_i = \frac{1}{p}(\widehat{\beta}-\widehat{\beta}_{(i)})^\top X^\top X (\widehat{\beta}-\widehat{\beta}_{(i)}), \end{equation*}\] onde \(\widehat{\beta}_{(i)}\) são as estimativas dos parâmetros retirando da amostra a \(i\)-ésima observação.

residualPlots(modelo0)

##            Test stat Pr(>|Test stat|)
## Quality      -1.0163           0.3433
## Color        -1.2601           0.2480
## Tukey test    0.9662           0.3340

marginalModelPlots(modelo0)

library(ggfortify)
autoplot(modelo0)

A observação 1 apresentam forte influência no resíduo e a observação 7 é de alta alavancagem (hat-value). Nesta situação não é uma opção a retirada de amostras, a final temos somente 11 observações.

residuos = residuals(modelo0, type="pearson")
summary(residuos)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  -6.672  -4.536  -1.092   0.000   3.722  10.190

sd(residuos)

## [1] 5.266463

plot(density(residuos), ylab = "Densidade estimada", main ="", xlim = c(-15,15), ylim = c(0,0.08))
grid()
rug(residuos, lwd = 2)
curve(dnorm(x, mean = mean(residuos), sd = sd(residuos)), add = TRUE, col = "red", lwd = 2)

Embora parecidas estas curvas, deveriamos observar uma amostra normal padrão. A curva em vermelho é a densidade normal de média \(\overline{x}=0.000\) e desvio padrão \(s=5.266463\). Isso significaria que os resíduos deveriam estar no intervalo (-3,3), não apresentar valos tão extremos e o desvio padrão deveria ser aproximadamente 1. Observemos as mesmas curvas agora com os resíduos estudentizados.

residuos = rstandard(modelo0)
summary(residuos)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -1.323630 -0.908229 -0.233398 -0.006416  0.831506  1.852916

sd(residuos)

## [1] 1.029432

plot(density(residuos), ylab = "Densidade estimada", main ="", xlim = c(-3,3), ylim = c(0,0.4))
grid()
rug(residuos, lwd = 2)
curve(dnorm(x, mean = mean(residuos), sd = sd(residuos)), add = TRUE, col = "red", lwd = 2)

É claro que alguma coisa está errada no modelo. Devemos lembrar que, Como qualquer modelo estatístico devemos fazer certas suposições para poder obtermos estimadores, realizar testes de hipóteses de interesse e verificarmos a adequação deste modelo aos dados, a final, a utilidade de qualquer modelo é sua qualidade em representar os dados.

Um modelo de regressão linear deve satisfazer as seguintes suposições:

• Suposição de Linearidade: esta suposição é implícita na definição do modelo e significa que cada valor observado \(y_i\) da variável resposta \(Y_i\) pode ser escrito como uma função linear da \(i\)-ésima linha da matriz \(X\), \(x_i\), isto é \[\begin{equation*} y_i=x_i\beta+\epsilon_1, \quad i=1,\cdots,n\cdot \end{equation*}\]
• Suposição Computacional: para encontrar estimador único do vetor de parâmetros é necessário que \((X^\top X)^{-1}\) exista.
• Suposição Distribucional: a estimação dos parâmetros é por mínimos quadrados ou por máxima verossimilhança e para isso assume-se que: (a) \(X\) é mensurada sem erros, (b) \(\epsilon_i\) não depende de \(x_i\), \(i = 1,\cdots,n\) e (c) \(\epsilon\sim N_n(0,\sigma^2\mbox{I})\), isto é, a distribuição de probabilidade do erro é normal multivariada de ordem \(n\), com média zero e matriz de variâncias e covariâncias \(\sigma^2\mbox{I}\).
• Suposição Implícita: todas as observações são igualmente possíveis de serem observadas e tem a mesma importância na determinação dos resultados de mínimos quadrados ou de máxima verossimilhança. Significa que todas as observações têm a mesma influência nas conclusões.

dados1 = cbind(dados, Prod = I(Quality*Color))
modelo1 = update(modelo0, . ~ . + Prod, data = dados1)
summary(modelo1)

## 
## Call:
## lm(formula = Price ~ Quality + Color + Prod, data = dados1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.0747 -2.4493 -0.4417  1.2411  7.9992 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  20.2071     8.1325   2.485   0.0419 *
## Quality      -0.3012     1.5086  -0.200   0.8474  
## Color         1.8313     1.2154   1.507   0.1756  
## Prod          0.7280     0.2522   2.887   0.0234 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.253 on 7 degrees of freedom
## Multiple R-squared:  0.9318, Adjusted R-squared:  0.9026 
## F-statistic:  31.9 on 3 and 7 DF,  p-value: 0.0001872

influenceIndexPlot(modelo1)

residuos = rstudent(modelo1)
summary(residuos)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -1.57610 -0.82120 -0.11225 -0.05373  0.36452  3.00850

sd(residuos)

## [1] 1.287887

plot(density(residuos), ylab = "Densidade estimada", main ="", xlim = c(-5,5), ylim = c(0,0.4))
grid()
rug(residuos, lwd = 2)
curve(dnorm(x, mean = mean(residuos), sd = sd(residuos)), add = TRUE, col = "red", lwd = 2)

library(qqplotr); library(ggplot2)
# p-p plot
gg = ggplot(data = data.frame(residuos),
mapping = aes(sample = residuos)) +
stat_pp_band() +
stat_pp_line() +
stat_pp_point() +
labs(x = "Função de distribuição empírica",
y = "Função de distribuição Normal padrão")

# q-q plot
gg1 = ggplot(data = data.frame(residuos),
mapping = aes(sample = residuos)) +
stat_qq_band() +
stat_qq_line() +
stat_qq_point() +
labs(x = "Quantis teóricos", y = "Quantis amostrais")

gg + gg1

Continuamos observando a influência das observações 1 e 7, mas o comportamento dos resíduos Pearson melhorou bastante.

# Teste Anderson-Darling de normalidade
library(DescTools)
AndersonDarlingTest(residuos, null = "pnorm", mean=0, sd=1)

## 
##  Anderson-Darling test of goodness-of-fit
##  Null hypothesis: Normal distribution
##  with parameters mean = 0.000, sd = 1.000
## 
## data:  residuos
## An = 0.62322, p-value = 0.6224

O teste de normalidade aceita esta hipóteses assim como nós aceitamos que todas as suposições do nomelo são cumpridas. São estas as afirmações que podemos tomar com o estudo dos resíduos.

Qualidade ou performance do modelo

Para avaliarmos a qualidade do ajuste utilizamos o coeficiente de determinação \(R^2\), o qual é amplamente favorável ao modelo reservado no objeto R modelo1.

AIC(modelo0,modelo1)

##         df      AIC
## modelo0  4 74.71813
## modelo1  5 68.09161

Também escolhemos o modelo em modelo1 utilizando o Critériio AIC. Mostramos agora os dados observados e os valores preditos.

plot(modelo1$model$Price,modelo1$fitted.values, xlab = "Price", ylab = "Valores preditos", 
     ylim = c(20,80), pch = 19)
grid()

Uma apresentação mais aprimorada seria a seguinte:

pred.int = predict(modelo1, interval = "prediction")  # adicionando predições aos dados
mydados = cbind(dados1, pred.int)                     # reta de regressão + intervalos confidenciais 

p <- ggplot(mydados, aes(Price,fit)) + ylab("Valores preditos") + ylim(c(20,80)) + 
  geom_point() +
  stat_smooth(method = lm)
# 3. Add prediction intervals
p + geom_line(aes(y = lwr), color = "red", linetype = "dashed")+
    geom_line(aes(y = upr), color = "red", linetype = "dashed")