1.7 Regressão robusta e estimadores M

Dados

x <- c(5.3,5.2,4.3,3.8,6.1,3.7,6.4,5.9,5.3,4.2,4.5,4.5)
y <- c(28.6,30.0,27.1,27.7,32.0,28.6,32.1,30.9,29.7,28.3,28.5,29.5)

dados <- data.frame(x = x, y = y)

Relação entre as variáveis

plot(dados$y ~ dados$x, xlab = "x", ylab = "y")

Ajuste dos modelos

ajuste_lm <- lm(y ~ x, data = dados)
ajuste_rlm <- MASS::rlm(y ~ x, data = dados, method = "MM")

pred_lm <- predict(ajuste_lm, type = "response")
pred_rlm <- predict(ajuste_rlm, type = "response")

plot(dados$y ~ dados$x, xlab = "x", ylab = "y")
lines(dados$x, pred_lm, col = "steelblue", type = "l", lwd = 2)
lines(dados$x, pred_rlm, col = "#ff5454", type = "l", lwd = 2)
legend("topleft", legend = c("LM", "RLM"),col = c("steelblue", "#ff5454"), lty = 1, lwd = 2)

A princípio, os dois modelos parecem se ajustar de maneira similar aos dados.

Bootstrap

Vamos considerar nossa amostra como a população alvo do estudo e retirar amostras com reposição de tamanho n dela. Esse processo será repetido muitas vezes, e em cada uma delas vamos armazenar a estimativa de interesse em uma posição de um vetor pré definido. Esse processo será feito de duas maneiras: dentro de um looping for, e com o auxílio do pacote boot.

Nesse problema, temos interesse em comparar a distribuição do coeficiente de regressão de uma regressão linear usual e de uma regressão linear robusta. Pra isso, além das estimativas, vamos analisar os erros associados a elas.

Passo a passo

vetor_b1_lm <- numeric(1000)
vetor_b1_rlm <- numeric(1000)

vetor_ep_lm <- numeric(1000)
vetor_ep_rlm <- numeric(1000)

set.seed(123)
for (i in 1:1000) {
  
  indices <- sample(nrow(dados), replace = TRUE)
  
  ajuste_lm <- lm(y ~ x, data = dados[indices, ])
  ajuste_rlm <- MASS::rlm(y ~ x, data = dados[indices, ], method = "MM")
  
  vetor_b1_lm[i] <- coef(ajuste_lm)[2]
  vetor_b1_rlm[i] <- coef(ajuste_rlm)[2]
  
  vetor_ep_lm[i] <- summary(ajuste_lm)$coefficients["x", "Std. Error"]
  vetor_ep_rlm[i] <- summary(ajuste_rlm)$coefficients["x", "Std. Error"]
  
}

Distribuição do coeficiente

plot(density(vetor_b1_lm), type = "l", col = "steelblue", main = "", xlim = c(-2,5))
lines(density(vetor_b1_rlm), col = "#ff5454")
legend("topleft", legend = c("LM", "RLM"),col = c("steelblue", "#ff5454"), lty = 1, lwd = 2)

Apesar das distribuições das estimativas serem muito parecidas e apresentarem pontos centrais muito próximos, o range das estimativas geradas pelo modelo de regressão robusta é significativamente maior.

Distribuição do erro padrão

plot(density(vetor_ep_lm), type = "l", col = "steelblue", main = "", xlim = c(-0.1,0.8), ylim = c(0,5.2))
lines(density(vetor_ep_rlm), col = "#ff5454")
legend("topleft", legend = c("LM", "RLM"),col = c("steelblue", "#ff5454"), lty = 1, lwd = 2)

Apesar dos erros associados as estimativas da regressão robusta apresentarem maior dispersão e um padrão mais instável, é importante notar que muitos deles ficaram mais próximos de zero em relação aos erros da regressão linear usual e o range de valores das duas distribuições é muito parecido.

Intervalos de confiança

# LM
quantile(vetor_b1_lm, c(0.025, 0.975))

##      2.5%     97.5% 
## 0.9725039 2.0745163

# RLM
quantile(vetor_b1_rlm, c(0.025, 0.975))

##      2.5%     97.5% 
## 0.7422977 2.3813168

Apesar dos intervalos serem parecidos, o intervalo para as estimativas geradas pela regressão linear usual tem um range menor.

Pacote boot

Com o auxílio do pacote boot, podemos realizar o processo anterior de uma forma mais fácil.

Regressão usual

set.seed(4)

b1_lm <- function(dados, indices){
  ajuste <- lm(y ~ x, data = dados[indices, ])
  coef(ajuste)[2]
}

boot.b1.lm <- boot(data = dados, statistic = b1_lm, R = 1000)

# Distribuição do coeficiente
plot(boot.b1.lm)

Assim como antes, as estimativas geradas pelo modelo de regressão linear usual seguiram um padrão aproximadamente normal, com um range moderado de valores. Apesar disso, o gráfico quantil-quantil indica uma pequena fuga da normalidade nas duas caldas (incidência de alguns valores mais distantes).

# Intervalo de confiança
boot.ci(boot.b1.lm, type = c("basic", "norm", "perc"))

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = boot.b1.lm, type = c("basic", "norm", "perc"))
## 
## Intervals : 
## Level      Normal              Basic              Percentile     
## 95%   ( 1.040,  2.095 )   ( 1.054,  2.127 )   ( 1.004,  2.077 )  
## Calculations and Intervals on Original Scale

Regressão robusta

set.seed(4)
b1_rlm <- function(dados, indices){
  ajuste <- MASS::rlm(y ~ x, data = dados[indices, ], method = "MM")
  coef(ajuste)[2]
}

boot.b1.rlm <- boot(data = dados, statistic = b1_rlm, R = 1000)

# Distribuição do coeficiente
plot(boot.b1.rlm)

Enquanto os valores das estimativas da regressão linear usual variam entre 0,5 e 3, os valores das estimativas da regressão robusta variam entre -10 e 5. Por isso, o uso de métodos de reamostragem se faz necessário, pois quando trabalhamos com apenas uma amostra e não temos condições de repetir o experimento, não sabemos se ela representa nossa população.

# Intervalo de confiança
boot.ci(boot.b1.rlm, type = c("basic", "norm", "perc"))

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 998 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = boot.b1.rlm, type = c("basic", "norm", "perc"))
## 
## Intervals : 
## Level      Normal              Basic              Percentile     
## 95%   ( 0.289,  2.903 )   ( 0.791,  2.727 )   ( 0.433,  2.369 )  
## Calculations and Intervals on Original Scale