Modelo de Regressão Linear Múltipla

Modelo

Considere o modelo de regressão linear múltipla: \[Y_i = \beta_0 + \beta_1 X_{i1} + \ldots + \beta_p X_{ip} + e_i, \quad e_i \stackrel{iid}{\sim} N(0,\sigma^2), \quad i=1,\ldots,n.\] Podemos reescrever o modelo acima da seguinte forma: \[Y_i = \boldsymbol{X}_i^T \boldsymbol{\beta} + e_i, \quad e_i \stackrel{iid}{\sim} N(0,\sigma^2), \quad i=1,\ldots,n,\] sendo \(\boldsymbol{X}_i^T = (1 \; X_{i1} \; \ldots \; X_{in})\) e \(\boldsymbol{\beta}^T = (\beta_0 \; \beta_1 \;\ldots \; \beta_p)\).

Novamente podemos reescrever o modelo da seguinte forma: \[\boldsymbol{Y} = \boldsymbol{X} \boldsymbol{\beta} + \boldsymbol{e}, \quad \boldsymbol{e} \sim N_n(\boldsymbol{0},\sigma^2 \boldsymbol{I}),\] sendo \(\boldsymbol{X}\) a matriz de delineamento (ou matriz de covariáveis) com \(n\) linhas e \(p+1\) colunas, incluindo o intercepto dada por \[\boldsymbol{X} = \begin{pmatrix} 1 & X_{11} & \cdots & X_{1p} \\ 1 & X_{21} & \cdots & X_{2p} \\ \vdots & \vdots & \ddots & \vdots \\ 1 & X_{n1} & \cdots & X_{np} \end{pmatrix};\]

\(\boldsymbol{\beta}^T = (\beta_0 ; \beta_1 ; \ldots ; \beta_p)\) o vetor de parâmetros desconhecidos;

\(\boldsymbol{Y} = (Y_1 ; Y_2 ; \ldots ; Y_n)^T\) o vetor de respostas;

\(\boldsymbol{e} = (e_1 ; e_2 ; \ldots ; e_n)^T\) o vetor de erros aleatórios, assumidos independentes e identicamente distribuídos, com média zero e variância constante \(\sigma^2\);

\(\boldsymbol{I}\) a matriz identidade de ordem \(n\).

1 Simulando conjunto de dados

Simule um conjunto de dados do modelo de regressão linear múltipla considerando n=1000 (tamanho amostral), \(\beta_0=10\), \(\beta_1=3\), \(\beta_2=-4\), \(\beta_3=0\), \(\sigma^2=10\) e que \(x_{1i} \stackrel{iid}{\sim} N(0,1)\), \(x_{2i} \stackrel{iid}{\sim} Exp(1)\) e \(x_{3i} \stackrel{iid}{\sim} t_5\).

# Configurações iniciais
n  = 1000  # Número de observações
# Simulando as covariáveis
x1 = rnorm(n, mean = 0, sd = 1)   # Covariável 1
x2 = rexp(n, rate = 1)            # Covariável 2
x3 = rt(n, df = 5)                # Covariável 3
x  = cbind(1, x1, x2, x3)

#definindo total de covariaveis
p = ncol(x)-1

# Definindo os coeficientes reais
beta_0 = 10    # Intercepto
beta_1 = 3     # Coeficiente para x1
beta_2 = -4  # Coeficiente para x2
beta_3 = 0     # Coeficiente para x3
beta   = c(beta_0, beta_1, beta_2, beta_3)
rm(beta_0, beta_1, beta_2, beta_3)

# Gerando a variável resposta com erro aleatório
sigma2 = 10
e = rnorm(n, mean = 0, sd = sqrt(sigma2))  # Erro aleatório
y = x %*% beta + e

# Criando um data frame com os dados
dados   = data.frame(y, x1=x1, x2=x2, x3=x3)

Analise o conjunto de dados simulado.

#analisando os dados gerados
par(mfrow=c(3,2))
hist(y, main="", freq=F, bty="n", ylab="densidade")
hist(e, main="", freq=F, bty="n", ylab="densidade")
plot(x1, y, bty="n")
plot(x2, y, bty="n")
plot(x3, y, bty="n")
plot(e, y, bty="n")

round(cor(dados),3)

##         y     x1     x2     x3
## y   1.000  0.508 -0.675  0.027
## x1  0.508  1.000 -0.016  0.002
## x2 -0.675 -0.016  1.000 -0.020
## x3  0.027  0.002 -0.020  1.000

2 Ajustando um MRLM

Ajuste um modelo de regressão linear múltipla e obtenha as estimativas pontual e intervalar para \(\beta_j\), \(j=0, \ldots, p\), considerando um nível de confiança de 90%.

#ajustando o modelo
ajuste      = lm(y~x1+x2+x3, data = dados) 
resultado   = summary(ajuste)

#vendo os resultados do ajuste
#ajuste
resultado

## 
## Call:
## lm(formula = y ~ x1 + x2 + x3, data = dados)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.6521  -2.3238   0.0322   2.1448  10.4868 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 10.03860    0.14215  70.619   <2e-16 ***
## x1           2.92703    0.10169  28.785   <2e-16 ***
## x2          -4.01916    0.10424 -38.556   <2e-16 ***
## x3           0.05441    0.07391   0.736    0.462    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.193 on 996 degrees of freedom
## Multiple R-squared:  0.7026, Adjusted R-squared:  0.7017 
## F-statistic: 784.3 on 3 and 996 DF,  p-value: < 2.2e-16

cbind(ajuste$coefficients, confint(ajuste, level=0.9), beta)

##                                 5 %       95 % beta
## (Intercept) 10.03859994  9.80456273 10.2726372   10
## x1           2.92703456  2.75962064  3.0944485    3
## x2          -4.01916195 -4.19078518 -3.8475387   -4
## x3           0.05441062 -0.06726939  0.1760906    0

Com base nas estimativas encontradas, você acredita que o modelo se ajustou bem aos dados?

3 Resíduos

3.1 Resíduos brutos ou ordinários

\[\hat{e}_i = Y_i - \hat{Y}_i = Y_i - \boldsymbol{X}_i^T\boldsymbol{\hat{\beta}}.\]

3.2 Resíduos padronizados

\[d_i = \frac{\hat{e}_i}{\sqrt{\hat{\sigma}^2}}.\]

Em amostras grandes considera-se que o ponto \(i\) é um potencial outlier quando \(|d_i|>3\), já que a distribuição de \(d_i\) será aproximadamente uma normal padrão, quando o modelo está bem ajustado.

3.3 Resíduos estudentizados internamente

\[r_i = \frac{\hat{e}_i}{\sqrt{\hat{\sigma}^2 (1-h_{ii})}},\] sendo \(h_{ii}\) o elemento da \(i\)-ésima linha e \(i\)-ésima coluna da matriz chapéu \(\boldsymbol{H} = \boldsymbol{X}(\boldsymbol{X}^T\boldsymbol{X})^{-1}\boldsymbol{X}^T\).

Este resíduo ajusta os resíduos padronizados levando em conta a alavancagem da observação de modo a corrigir o fato de que observações com alta alavancagem tendem a distorcer o ajuste do modelo.

Geralmente considera-se que o ponto \(i\) é um potencial outlier quando \(|r_i|>3\).

3.4 Resíduos estudentizados externamente

\[r_i^{*} = \frac{\hat{e}_i}{\sqrt{\hat{\sigma}^2_{-i} (1-h_{ii})}}, \mbox{ sendo } \hat{\sigma}^2_{-i} = \frac{(n-p-1)\hat{\sigma}^2 - \frac{\hat{e}_i^2}{1-h_{ii}}}{n-p-2}.\]

Este resíduo ajusta os resíduos estudentizados internamente removendo a influência da observação na estimativa da variância residual. Isso fornece uma avaliação mais precisa do impacto do ponto \(i\) no ajuste do modelo.

Como \(r_i^* \stackrel{\cdot}{\sim} t_{n-p-2}\), considera-se que o ponto \(i\) é um potencial outlier quando \(|r_i^*|>t_{\alpha/2, n-p-2}\), sendo \(t_{\alpha/2, n-p-2}\) o quantil \((1-\alpha/2) \times 100\%\) da distribuição t-Student com \(n-p-2\) graus de liberdade.

3.5 Tarefa

Obtenha os resíduos brutos (ou também conhecidos como ordinários), os resíduos padronizados, os resíduos estudentizados internamente e os resíduos externamente estudentizados. Analise o ajuste do modelo com base nestes 4 resíduos. Há alguma diferença na análise?

#residuos brutos ou ordinarios
residuo_ordinario = dados$y - ajuste$fitted.values
  #comparando com funcoes do R
  all.equal(residuo_ordinario, ajuste$residuals)

## [1] TRUE

  all.equal(residuo_ordinario, residuals(ajuste))

## [1] TRUE

#residuos padronizados
residuo_padronizado = residuo_ordinario/sqrt(resultado$sigma^2)

#residuos estudentizados internamente
H     = x%*%solve(t(x)%*%x)%*%t(x)
h_ii  = diag(H)
residuo_estudentizado = residuo_padronizado / sqrt(1-h_ii)
  #comparando com funcoes do R
  all.equal(residuo_estudentizado, rstandard(ajuste))

## [1] TRUE

#residuos estudentizados externamente
  # Calcular sigma_i_sq para cada observação
  MQE_i = ((n-p-1)*resultado$sigma^2 - residuo_ordinario^2/(1-h_ii))/(n-p-2)
residuo_est_ext = residuo_ordinario / sqrt(MQE_i * (1-h_ii) )
  #comparando com funcoes do R
  all.equal(residuo_est_ext, rstudent(ajuste))

## [1] TRUE

#comparando as variancias  
cbind(var(residuo_ordinario),
var(residuo_padronizado),
var(residuo_estudentizado),
var(residuo_est_ext))

##          [,1]     [,2]     [,3]   [,4]
## [1,] 10.16462 0.996997 1.001125 1.0031

par(mfrow = c(2, 2))  # Dividir a janela de gráficos
plot(residuo_ordinario, main = "Resíduos Ordinários")
plot(residuo_padronizado, main = "Resíduos Padronizados")
plot(residuo_estudentizado, main = "Res. Estudentizados Internamente")
plot(residuo_est_ext, main = "Res. Estudentizados Externamente")

par(mfrow = c(2, 2))  # Dividir a janela de gráficos
qqnorm(residuo_ordinario, xlab="quantis teóricos", ylab="quantis amostrais", bty="n", main="Resíduos ordinários")
qqline(residuo_ordinario, col = "red", lwd=2)
qqnorm(residuo_padronizado, xlab="quantis teóricos", ylab="quantis amostrais", bty="n", main = "Resíduos Padronizados")
qqline(residuo_padronizado, col = "red", lwd=2)
qqnorm(residuo_estudentizado, xlab="quantis teóricos", ylab="quantis amostrais", bty="n", main = "Res. Estudentizados Internamente")
qqline(residuo_estudentizado, col = "red", lwd=2)
qqnorm(residuo_est_ext, xlab="quantis teóricos", ylab="quantis amostrais", bty="n", main = "Res. Estudentizados Externamente")
qqline(residuo_est_ext, col = "red", lwd=2)

par(mfrow=c(2,2))
hist(residuo_ordinario, main="", freq=F, ylab="densidade", xlab="resíduos ordinários")
curve(dnorm(x, 0, sqrt(resultado$sigma^2)), add=T, col="red", lwd=2)
hist(residuo_padronizado, main="", freq=F, ylab="densidade", xlab="resíduos padronizados")
curve(dnorm(x, 0, 1), add=T, col="red", lwd=2)
hist(residuo_estudentizado, main="", freq=F, ylab="densidade", xlab="res. estudentizados internamente")
curve(dnorm(x, 0, 1), add=T, col="red", lwd=2)
hist(residuo_est_ext, main="", freq=F, ylab="densidade", xlab="res. estudentizados externamente")
curve(dnorm(x, 0, 1), add=T, col="red", lwd=2)

shapiro.test(residuo_ordinario)

## 
##  Shapiro-Wilk normality test
## 
## data:  residuo_ordinario
## W = 0.99924, p-value = 0.9657

shapiro.test(residuo_padronizado)

## 
##  Shapiro-Wilk normality test
## 
## data:  residuo_padronizado
## W = 0.99924, p-value = 0.9657

shapiro.test(residuo_estudentizado)

## 
##  Shapiro-Wilk normality test
## 
## data:  residuo_estudentizado
## W = 0.99924, p-value = 0.9652

shapiro.test(residuo_est_ext)

## 
##  Shapiro-Wilk normality test
## 
## data:  residuo_est_ext
## W = 0.99925, p-value = 0.9674

ks.test(residuo_ordinario,"pnorm",0,sd(residuo_ordinario))

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  residuo_ordinario
## D = 0.021184, p-value = 0.7606
## alternative hypothesis: two-sided

ks.test(residuo_padronizado,"pnorm",0,sd(residuo_padronizado))

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  residuo_padronizado
## D = 0.021184, p-value = 0.7606
## alternative hypothesis: two-sided

ks.test(residuo_estudentizado,"pnorm",0,sd(residuo_estudentizado))

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  residuo_estudentizado
## D = 0.021567, p-value = 0.7409
## alternative hypothesis: two-sided

ks.test(residuo_est_ext,"pnorm",0,sd(residuo_est_ext))

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  residuo_est_ext
## D = 0.021295, p-value = 0.7549
## alternative hypothesis: two-sided

par(mfrow=c(3,2))
plot(ajuste$fitted.values, residuo_ordinario)
abline(h=0, lwd=2, col="red", lty=3)
plot(x1, residuo_ordinario)
abline(h=0, lwd=2, col="red", lty=3)
plot(x2, residuo_ordinario)
abline(h=0, lwd=2, col="red", lty=3)
plot(x3, residuo_ordinario)
abline(h=0, lwd=2, col="red", lty=3)
plot(y, residuo_ordinario)
abline(h=0, lwd=2, col="red", lty=3)
plot(residuo_ordinario)
abline(h=0, lwd=2, col="red", lty=3)

# Testando Homocedasticidade (Teste de Breusch-Pagan)
# Instalar pacote necessário (se ainda não instalado)
# install.packages("lmtest")
library(lmtest)
bptest(ajuste)

## 
##  studentized Breusch-Pagan test
## 
## data:  ajuste
## BP = 3.2396, df = 3, p-value = 0.3561

#Interpretação: Um p-valor pequeno indica a presença de heterocedasticidade, o que significa que a variância dos resíduos não é constante.

#Avaliando a Independência dos Resíduos (Durbin-Watson)
dwtest(ajuste)

## 
##  Durbin-Watson test
## 
## data:  ajuste
## DW = 2.0495, p-value = 0.7837
## alternative hypothesis: true autocorrelation is greater than 0

#Interpretação: Um valor de estatística Durbin-Watson próximo de 2 sugere que não há autocorrelação dos resíduos.

3.6 Outliers

Com base nos resíduos calculados, há outliers?

#identificando possiveis outliers
  i_rp=which(abs(residuo_padronizado)>3)
  i_rei=which(abs(residuo_estudentizado)>3)
  #i_ree=which(abs(residuo_est_ext)>qt(1-0.05/2, n-p-2))
  i_ree=which(abs(residuo_est_ext)>3)

  for (i in i_rp){
    cat("Índice", i,
    "| Resíduo Padronizado:", residuo_padronizado[i], 
      "| y:", y[i], 
      "| e:", e[i], 
      "| prob:", round(pnorm(e[i],0,sqrt(sigma2)),3), "\n")
  }

## Índice 117 | Resíduo Padronizado: -3.649279 | y: -2.749807 | e: -11.46956 | prob: 0 
## Índice 182 | Resíduo Padronizado: 3.28432 | y: 23.51953 | e: 10.50857 | prob: 1

  for (i in i_rei){
    cat("Índice", i,
    "| Resíduo Est. Int.:", residuo_estudentizado[i], 
      "| y:", y[i], 
      "| e:", e[i], 
      "| prob:", round(pnorm(e[i],0,sqrt(sigma2)),3), "\n")
  }

## Índice 117 | Resíduo Est. Int.: -3.658393 | y: -2.749807 | e: -11.46956 | prob: 0 
## Índice 182 | Resíduo Est. Int.: 3.290389 | y: 23.51953 | e: 10.50857 | prob: 1

  for (i in i_ree){
    cat("Índice", i,
    "| Resíduo Est. Ext.:", residuo_est_ext[i], 
      "| y:", y[i], 
      "| e:", e[i], 
      "| prob:", round(pnorm(e[i],0,sqrt(sigma2)),3), "\n")
  }

## Índice 117 | Resíduo Est. Ext.: -3.681374 | y: -2.749807 | e: -11.46956 | prob: 0 
## Índice 182 | Resíduo Est. Ext.: 3.306758 | y: 23.51953 | e: 10.50857 | prob: 1 
## Índice 902 | Resíduo Est. Ext.: 3.00007 | y: 21.65834 | e: 9.571021 | prob: 0.999

3.7 Resíduos parciais

\[t_{ki} = \hat{e}_i + X_{ik} \hat{\beta}_k.\]

Calcule os resíduos parciais e analise-os.

#residuos parciais
  residuo_parcial_1 = residuo_ordinario + x1*ajuste$coefficients[2]
  residuo_parcial_2 = residuo_ordinario + x2*ajuste$coefficients[3]
  residuo_parcial_3 = residuo_ordinario + x3*ajuste$coefficients[4]

  # Gráfico de resíduos parciais
  par(mfrow=c(1,3))
  plot(dados$x1, residuo_parcial_1, 
     xlab = "x1", ylab = "Resíduos Parciais", 
     main = "Gráfico de Resíduos Parciais", bty="n")
  abline(lm(residuo_parcial_1 ~ dados$x1), col = "red") 
  plot(dados$x2, residuo_parcial_2, 
     xlab = "x2", ylab = "Resíduos Parciais", 
     main = "Gráfico de Resíduos Parciais", bty="n")
  abline(lm(residuo_parcial_2 ~ dados$x2), col = "red") 
  plot(dados$x3, residuo_parcial_3, 
     xlab = "x3", ylab = "Resíduos Parciais", 
     main = "Gráfico de Resíduos Parciais", bty="n")
  abline(lm(residuo_parcial_3 ~ dados$x3), col = "red")

Aula de Modelos Lineares I

Patrícia Lusié Velozo da Costa

2025-06-04