Nesta aula abordaremos Outliers, Pontos de Alavancagem e Pontos Influentes. E também inclusão de variáveis categóricas na modelagem.

Modelo de Regressão Linear Múltipla

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\), \(\sigma^2=1\) e que \(x_{1i} \stackrel{iid}{\sim} N(0,1)\) e \(x_{2i} \stackrel{iid}{\sim} N(0,1)\).

# 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 = rnorm(n, mean = 0, sd = 1)   # Covariável 2            
x  = cbind(1, x1, x2)

#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   = c(beta_0, beta_1, beta_2)
rm(beta_0, beta_1, beta_2)

# Gerando a variável resposta com erro aleatório
sigma2 = 1
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)

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(e, y, bty="n")

round(cor(dados),3)

##         y     x1     x2
## y   1.000  0.576 -0.798
## x1  0.576  1.000 -0.007
## x2 -0.798 -0.007  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, data = dados) 
resultado   = summary(ajuste)

#note que eu encontraria o mesmo valor de p
#p = length(coef(ajuste)) - 1

#vendo os resultados do ajuste
#ajuste
resultado

## 
## Call:
## lm(formula = y ~ x1 + x2, data = dados)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.0832 -0.6858 -0.0583  0.6965  4.2997 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 10.02837    0.03151  318.27   <2e-16 ***
## x1           3.03011    0.03247   93.33   <2e-16 ***
## x2          -4.02167    0.03095 -129.93   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9963 on 997 degrees of freedom
## Multiple R-squared:  0.9628, Adjusted R-squared:  0.9627 
## F-statistic: 1.288e+04 on 2 and 997 DF,  p-value: < 2.2e-16

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

##                             5 %      95 % beta
## (Intercept) 10.028372  9.976496 10.080248   10
## x1           3.030110  2.976655  3.083564    3
## x2          -4.021675 -4.072633 -3.970716   -4

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

residuo_ordinario = ajuste$residuals
residuo_padronizado = residuo_ordinario/sqrt(resultado$sigma^2)
residuo_estudentizado = rstandard(ajuste)
residuo_est_ext = rstudent(ajuste)


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.99777, p-value = 0.1993

shapiro.test(residuo_padronizado)

## 
##  Shapiro-Wilk normality test
## 
## data:  residuo_padronizado
## W = 0.99777, p-value = 0.1993

shapiro.test(residuo_estudentizado)

## 
##  Shapiro-Wilk normality test
## 
## data:  residuo_estudentizado
## W = 0.99778, p-value = 0.203

shapiro.test(residuo_est_ext)

## 
##  Shapiro-Wilk normality test
## 
## data:  residuo_est_ext
## W = 0.99774, p-value = 0.1912

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

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

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

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

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

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  residuo_estudentizado
## D = 0.026597, p-value = 0.479
## 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.026575, p-value = 0.4801
## 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(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 = 1.3184, df = 2, p-value = 0.5173

#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 = 1.9228, p-value = 0.1107
## 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 Pontos a serem analisados

Outliers: observações que não são bem ajustadas pelo modelo;
Observações influentes: observações que afetam alguma propriedade do modelo ajustado de maneira substancial;
Leverage (ponto de alavanca): é um ponto extremo no espaço das variáveis explicativas.

Uma mesma observação pode apresentar duas ou até as três características simultaneamente.

3.1 Identificando outliers

Identifique possíveis outliers — observações com comportamento discrepante em relação ao modelo ajustado — por meio dos resíduos padronizados, dos estudentizados internamente e dos estudentizados externamente. Analise os pontos identificados.

Geralmente considera-se que uma observação é outlier, quando o valor absoluto do resíduo (padronizado, estudentizados internamente e estudentizados externamente) é maior do que 3.

#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 150 | Resíduo Padronizado: 4.315688 | y: 14.05511 | e: 4.315985 | prob: 1 
## Índice 819 | Resíduo Padronizado: -3.094638 | y: 2.335637 | e: -3.070488 | prob: 0.001

  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 150 | Resíduo Est. Int.: 4.319754 | y: 14.05511 | e: 4.315985 | prob: 1 
## Índice 819 | Resíduo Est. Int.: -3.101004 | y: 2.335637 | e: -3.070488 | prob: 0.001

  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 150 | Resíduo Est. Ext.: 4.358568 | y: 14.05511 | e: 4.315985 | prob: 1 
## Índice 819 | Resíduo Est. Ext.: -3.114505 | y: 2.335637 | e: -3.070488 | prob: 0.001

3.2 Identificando pontos de alavancagem

Geralmente, considera-se uma observação \(i\) como sendo ponto de alavancagem quando \(H_{ii}\) é maior que \(2(p+1)/n\). Identifique possíveis pontos de alavancagem e analise os pontos identificados.

#calculando matriz chapeu
H     = x%*%solve(t(x)%*%x)%*%t(x)
h_ii  = diag(H)
h_ii_2 = hatvalues(ajuste)
round(summary(h_ii_2 - h_ii),10)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0       0       0       0       0       0

#traço da matriz chapéu é igual ao número de coeficientes betas
c(sum(diag(H)), p+1)

## [1] 3 3

#Indices das unidades que sao pontos de alavanca
i_pa = which(h_ii>2*(p+1)/n) #limite moderado
#i_pa = which(h_ii>3*(p+1)/n) #limite conservador
length(i_pa) #qtde de pontos de alavanca

## [1] 93

#analisando a variacao das covariaveis
rbind(round(c(min(x[,2]), max(x[,2])),4),
      round(c(min(x[,3]), max(x[,3])),4))

##         [,1]   [,2]
## [1,] -2.8668 3.0032
## [2,] -2.9608 3.3808

rbind(round(quantile(x[,2],c(0.025,0.975)),4),
      round(quantile(x[,3],c(0.025,0.975)),4))

##         2.5%  97.5%
## [1,] -2.0738 1.7946
## [2,] -2.0277 2.0769

#analisando os pontos identificados como alavanca
for (i in i_pa){
  p1 <- pnorm(x[i,2], mean(x[,2]), sd(x[,2]))
  p2 <- pnorm(x[i,3], mean(x[,3]), sd(x[,3]))
  p_conj <- p1 * p2

    cat("Índice", i,
    "| y:", y[i], 
    "| x1:", x[i,2], 
    "| x2:", x[i,3], 
      "| e:", e[i], 
      "| prob:", round(p_conj,3), "\n")
}

## Índice 3 | y: 22.88738 | x1: 0.822308 | x2: -2.392648 | e: 0.8498614 | prob: 0.007 
## Índice 26 | y: 4.372439 | x1: 0.8245212 | x2: 2.559005 | e: 2.134894 | prob: 0.801 
## Índice 33 | y: -0.3030749 | x1: 0.1894809 | x2: 2.730007 | e: 0.04850867 | prob: 0.58 
## Índice 39 | y: -1.843123 | x1: -1.457581 | x2: 1.739871 | e: -0.5108963 | prob: 0.065 
## Índice 40 | y: -4.646651 | x1: -2.240705 | x2: 2.162852 | e: 0.7268695 | prob: 0.011 
## Índice 46 | y: -4.380257 | x1: -1.865102 | x2: 2.124268 | e: -0.2878798 | prob: 0.028 
## Índice 48 | y: 20.78487 | x1: 0.400483 | x2: -2.915826 | e: -2.079881 | prob: 0.001 
## Índice 61 | y: -1.422379 | x1: 0.3082229 | x2: 3.38082 | e: 1.176234 | prob: 0.629 
## Índice 89 | y: 0.8064964 | x1: -2.491961 | x2: 0.7959731 | e: 1.466273 | prob: 0.004 
## Índice 90 | y: 4.839448 | x1: -2.247318 | x2: 0.1761963 | e: 2.286188 | prob: 0.006 
## Índice 100 | y: 19.44489 | x1: 1.798639 | x2: -1.312162 | e: -1.199672 | prob: 0.095 
## Índice 103 | y: 12.93055 | x1: 2.289858 | x2: 0.943771 | e: -0.1639431 | prob: 0.814 
## Índice 120 | y: 21.98769 | x1: 2.018755 | x2: -1.601605 | e: -0.4749965 | prob: 0.056 
## Índice 138 | y: -3.81883 | x1: -1.783498 | x2: 1.823389 | e: -1.174778 | prob: 0.033 
## Índice 156 | y: 6.312673 | x1: -2.306204 | x2: -0.5224653 | e: 1.141424 | prob: 0.003 
## Índice 157 | y: 6.155317 | x1: 1.295011 | x2: 1.91415 | e: -0.07311671 | prob: 0.883 
## Índice 174 | y: 13.62979 | x1: 2.138703 | x2: 0.9686366 | e: 1.088225 | prob: 0.816 
## Índice 176 | y: 12.7589 | x1: -1.408311 | x2: -2.114393 | e: -1.473743 | prob: 0.001 
## Índice 179 | y: 3.249744 | x1: -2.220952 | x2: 0.1824405 | e: 0.6423621 | prob: 0.007 
## Índice 182 | y: 7.563814 | x1: -2.201099 | x2: -0.4954539 | e: 2.185296 | prob: 0.004 
## Índice 204 | y: -3.744922 | x1: -1.319631 | x2: 1.951517 | e: -1.979962 | prob: 0.087 
## Índice 209 | y: 5.338963 | x1: 1.232997 | x2: 2.026899 | e: -0.2524345 | prob: 0.879 
## Índice 211 | y: -5.111897 | x1: -2.685922 | x2: 1.468242 | e: -1.181162 | prob: 0.003 
## Índice 216 | y: -0.9300773 | x1: -1.948533 | x2: 1.174112 | e: -0.3880308 | prob: 0.02 
## Índice 271 | y: 1.832247 | x1: 1.633218 | x2: 2.884802 | e: -1.528198 | prob: 0.953 
## Índice 300 | y: 0.3328915 | x1: -2.325755 | x2: 0.9034867 | e: 0.9241045 | prob: 0.007 
## Índice 313 | y: 15.93894 | x1: -0.9085907 | x2: -2.241725 | e: -0.3021872 | prob: 0.002 
## Índice 331 | y: 21.37837 | x1: 1.752023 | x2: -1.519982 | e: 0.04236898 | prob: 0.065 
## Índice 338 | y: -1.074099 | x1: -1.841141 | x2: 1.534321 | e: 0.5866082 | prob: 0.028 
## Índice 341 | y: 12.28461 | x1: 2.900913 | x2: 1.614645 | e: 0.04044677 | prob: 0.941 
## Índice 363 | y: 7.852052 | x1: 1.675807 | x2: 1.620747 | e: -0.6923808 | prob: 0.905 
## Índice 366 | y: 6.37656 | x1: -2.518388 | x2: -0.805974 | e: 0.7078277 | prob: 0.001 
## Índice 372 | y: 8.637359 | x1: -2.180022 | x2: -1.537632 | e: -0.9731036 | prob: 0.001 
## Índice 377 | y: 13.56357 | x1: 1.960508 | x2: 1.156053 | e: 2.306261 | prob: 0.852 
## Índice 396 | y: 20.06076 | x1: 0.7921559 | x2: -2.259573 | e: -1.354 | prob: 0.01 
## Índice 410 | y: 21.88905 | x1: 1.213485 | x2: -2.227375 | e: -0.6609099 | prob: 0.013 
## Índice 413 | y: 13.64113 | x1: -0.9605667 | x2: -2.094473 | e: -1.855062 | prob: 0.003 
## Índice 414 | y: 1.674838 | x1: 0.4760483 | x2: 2.396855 | e: -0.165889 | prob: 0.686 
## Índice 434 | y: 3.61704 | x1: -2.399741 | x2: -0.3362954 | e: -0.5289195 | prob: 0.003 
## Índice 452 | y: -2.667834 | x1: -1.465803 | x2: 1.763637 | e: -1.215876 | prob: 0.064 
## Índice 462 | y: -0.3080204 | x1: -2.866833 | x2: 0.3617679 | e: -0.2604502 | prob: 0.001 
## Índice 476 | y: 24.27919 | x1: 1.144516 | x2: -2.506822 | e: 0.8183592 | prob: 0.006 
## Índice 488 | y: 3.246747 | x1: -2.572414 | x2: 0.02820272 | e: 1.076801 | prob: 0.002 
## Índice 489 | y: 16.93999 | x1: -0.5264014 | x2: -2.230705 | e: -0.4036242 | prob: 0.004 
## Índice 498 | y: 9.527074 | x1: 1.751458 | x2: 1.527081 | e: 0.3810254 | prob: 0.9 
## Índice 530 | y: -0.5933575 | x1: -0.3113079 | x2: 2.763298 | e: 1.393759 | prob: 0.378 
## Índice 537 | y: 12.47597 | x1: -1.739338 | x2: -1.476201 | e: 1.789186 | prob: 0.003 
## Índice 542 | y: -2.82063 | x1: -2.684862 | x2: 0.859052 | e: -1.329838 | prob: 0.002 
## Índice 549 | y: 7.763972 | x1: -2.073311 | x2: -0.8722457 | e: 0.494924 | prob: 0.003 
## Índice 568 | y: 11.16764 | x1: -1.956702 | x2: -2.103578 | e: -1.376565 | prob: 0 
## Índice 575 | y: -1.590176 | x1: -0.3965574 | x2: 2.740846 | e: 0.5628814 | prob: 0.345 
## Índice 585 | y: 21.87087 | x1: 1.170689 | x2: -1.952641 | e: 0.5482342 | prob: 0.024 
## Índice 586 | y: -2.017808 | x1: -2.538192 | x2: 0.9992906 | e: -0.4060681 | prob: 0.004 
## Índice 593 | y: 5.172499 | x1: 1.206872 | x2: 2.424477 | e: 1.249791 | prob: 0.887 
## Índice 598 | y: 5.096904 | x1: 1.203115 | x2: 2.140013 | e: 0.04760904 | prob: 0.878 
## Índice 601 | y: 8.422015 | x1: -1.972759 | x2: -0.959131 | e: 0.5037684 | prob: 0.004 
## Índice 621 | y: 11.93359 | x1: 2.34043 | x2: 1.548027 | e: 1.10441 | prob: 0.928 
## Índice 646 | y: 19.4845 | x1: 0.08541139 | x2: -2.314082 | e: -0.02806269 | prob: 0.006 
## Índice 657 | y: 21.02987 | x1: -0.05752975 | x2: -2.585506 | e: 0.8604296 | prob: 0.003 
## Índice 662 | y: 21.56296 | x1: 0.1287761 | x2: -2.960768 | e: -0.6664372 | prob: 0.001 
## Índice 663 | y: -0.9963547 | x1: -0.5963034 | x2: 2.232162 | e: -0.2787951 | prob: 0.27 
## Índice 676 | y: 11.35281 | x1: 2.1609 | x2: 1.244462 | e: -0.152044 | prob: 0.877 
## Índice 677 | y: -0.4298555 | x1: -0.3057459 | x2: 2.353785 | e: -0.09747808 | prob: 0.377 
## Índice 685 | y: 3.385071 | x1: -2.428626 | x2: 0.03162263 | e: 0.7974384 | prob: 0.003 
## Índice 703 | y: 15.14957 | x1: -1.160589 | x2: -1.979545 | e: 0.713153 | prob: 0.003 
## Índice 707 | y: 5.153172 | x1: 2.009883 | x2: 2.743067 | e: 0.09579346 | prob: 0.978 
## Índice 752 | y: 7.949978 | x1: -1.923489 | x2: -1.107813 | e: -0.7108077 | prob: 0.003 
## Índice 766 | y: -0.3785491 | x1: -2.263106 | x2: 0.5552888 | e: -1.368077 | prob: 0.007 
## Índice 789 | y: 7.382225 | x1: 1.450287 | x2: 2.171232 | e: 1.716294 | prob: 0.918 
## Índice 794 | y: 21.35409 | x1: 0.05735863 | x2: -2.835653 | e: -0.1606017 | prob: 0.001 
## Índice 796 | y: 21.3349 | x1: -0.07474359 | x2: -2.758794 | e: 0.5239594 | prob: 0.002 
## Índice 813 | y: -1.336483 | x1: -2.138389 | x2: 1.18287 | e: -0.1898367 | prob: 0.012 
## Índice 815 | y: 17.83998 | x1: -0.4100886 | x2: -2.435321 | e: -0.6710386 | prob: 0.003 
## Índice 824 | y: 6.936204 | x1: 1.618409 | x2: 1.598688 | e: -1.524271 | prob: 0.897 
## Índice 825 | y: 11.67566 | x1: 2.204862 | x2: 0.7494605 | e: -1.94108 | prob: 0.758 
## Índice 841 | y: 8.512864 | x1: -2.275328 | x2: -1.502439 | e: -0.6709092 | prob: 0.001 
## Índice 843 | y: 22.21702 | x1: 3.003199 | x2: -1.087249 | e: -1.14157 | prob: 0.141 
## Índice 845 | y: -0.3807564 | x1: -2.468563 | x2: 0.8083814 | e: 0.2584572 | prob: 0.004 
## Índice 847 | y: 11.85217 | x1: -1.567716 | x2: -1.780042 | e: -0.5648487 | prob: 0.002 
## Índice 870 | y: -6.35969 | x1: -1.318053 | x2: 2.907709 | e: -0.774694 | prob: 0.089 
## Índice 871 | y: 9.228393 | x1: -2.366597 | x2: -1.754638 | e: -0.6903687 | prob: 0 
## Índice 873 | y: 18.32846 | x1: -0.116331 | x2: -2.55633 | e: -1.54787 | prob: 0.003 
## Índice 891 | y: 21.92543 | x1: 1.493356 | x2: -1.87722 | e: -0.06351689 | prob: 0.03 
## Índice 892 | y: 20.05597 | x1: -0.07477547 | x2: -2.378043 | e: 0.7681243 | prob: 0.005 
## Índice 899 | y: 16.59777 | x1: 2.497181 | x2: 0.06343914 | e: -0.6400205 | prob: 0.519 
## Índice 930 | y: -1.410353 | x1: -1.12791 | x2: 2.629543 | e: 2.491549 | prob: 0.125 
## Índice 931 | y: 1.24323 | x1: -2.556324 | x2: 0.007637847 | e: -1.057245 | prob: 0.002 
## Índice 946 | y: 21.08615 | x1: 0.7621125 | x2: -2.414801 | e: -0.859395 | prob: 0.007 
## Índice 975 | y: 3.690461 | x1: -2.184054 | x2: -0.1850804 | e: -0.4976997 | prob: 0.005 
## Índice 980 | y: 2.217786 | x1: 0.4519589 | x2: 2.43208 | e: 0.5902281 | prob: 0.678 
## Índice 982 | y: -3.078187 | x1: -1.290897 | x2: 2.212582 | e: -0.3551684 | prob: 0.092 
## Índice 988 | y: -0.08953675 | x1: 0.3023031 | x2: 2.526247 | e: -0.8914595 | prob: 0.623 
## Índice 990 | y: 14.40146 | x1: 2.188525 | x2: 0.4806503 | e: -0.2415134 | prob: 0.671

3.3 Identificando pontos influentes

Identifique possíveis pontos influentes usando os métodos abaixo:

Distância de Cook: \(D_i = \frac{\hat{e}_i^2}{(p+1) \hat{\sigma}^2}\frac{H_{ii}}{(1-H_{ii})^2}\). A regra usual é classificar como influente observações tais que \(D_i > 1\) ou \(D_i >4/n\).
DFBETAS

\[\text{DFBETAS}_{ij} = \frac{\hat{\beta}_j - \hat{\beta}_{j(-i)}}{\hat{\sigma}_{(-i)} \sqrt{C_{jj}}},\] sendo \(\beta_j\) O coeficiente estimado no modelo com todas as observações, \(\beta_{j(-i)}\) O coeficiente estimado no modelo sem a \(i\)-ésima observação, \(\hat{\sigma}_{(-i)}\) O desvio padrão residual estimado sem a \(i\)-ésima observação, \(C_{jj}\) o \(j\)-ésimo elemento diagonal da matriz de covariância dos coeficientes \((\boldsymbol{X}^T\boldsymbol{X})^{-1}\) com todas as observações.

A regra usual é classificar como influente observações tais que \(|\text{DFBETAS}_{ij}|>1\) ou \(|\hat{\beta}_j - \hat{\beta}_{j(-i)}| > \frac{2}{\sqrt{n}}\).

DFFITS

\[\text{DFFITS}_i = \frac{\hat{Y}_i - \hat{Y}_{i(-i)}}{\hat{\sigma}_{(-i)} \sqrt{h_{ii}}},\] sendo \(\hat{Y}_{i(-i)}\) o valor previsto para a \(i\)-ésima observação considerando o modelo ajustado sem essa observação.

A regra usual é classificar como influente observações tais que \(|\text{DFFITS}_{i}|>2\sqrt{\frac{p+1}{n}}\).

#######################
#Distancia de Cook
#######################

dist_cook = residuo_ordinario^2/((p+1)*resultado$sigma^2)*h_ii/((1-h_ii)^2)
all.equal(dist_cook, cooks.distance(ajuste))

## [1] TRUE

#Identificando pontos influentes
i_pi = which(dist_cook > 1)
#i_pi = which(dist_cook > (4/n))
length(i_pi)

## [1] 0

for (i in i_pi){
    cat("Índice", i,
    "| y:", y[i], 
    "| x1:", x[i,2], 
    "| x2:", x[i,3], 
      "| e:", e[i], 
      "| prob:", round(pnorm(e[i],0,sqrt(sigma2)),3), "\n")
}

#######################
# DFBETAS
#######################

# Calculando DFBETA para todas as observações

  # Função para calcular DFBETA para uma observação i
  funcaodfbeta = function(i) {
    #Excluindo a i-ésima observação
    dados_menos_i = data.frame(y = dados$y[-i], 
                               x1 = dados$x1[-i], 
                               x2 = dados$x2[-i])
    
    #Reajustando o modelo sem a observação i
    ajuste_menos_i = lm(y ~ x1 + x2, data=dados_menos_i)
    resultado_menos_i = summary(ajuste_menos_i)
    
    #Calculando a matriz de covariância
    C_jj = solve(t(x)%*%x)
    # Diferença nos coeficientes
    diferenca = ajuste$coefficients - ajuste_menos_i$coefficients
    # Calculando DFBETA
    dfbeta_i = diferenca/
      (resultado_menos_i$sigma * sqrt(diag(C_jj)))
    
    # Retorna uma lista com as duas matrizes
    return(list(diferenca = diferenca, dfbeta = dfbeta_i))
  }

# Calculando diferencas e DFBETAs para todas as observações
vdiferencas  = matrix(NA, n, p + 1)
vdfbetas     = matrix(NA, n, p + 1)
for (i in 1:n){
  aux              = funcaodfbeta(i)
  vdiferencas[i,]  = aux$diferenca
  vdfbetas[i, ]    = aux$dfbeta
}

# Comparando resultados com funções predefinidas
vdiferencas2  = dfbeta(ajuste)
vdfbetas2     = dfbetas(ajuste)

for (j in 1:(p+1)) {
  print(all.equal(as.vector(vdfbetas[,j]), as.vector(vdfbetas2[,j])))
}

## [1] TRUE
## [1] TRUE
## [1] TRUE

for (j in 1:(p + 1)) {
  print(all.equal(as.vector(vdiferencas[, j]), as.vector(vdiferencas2[, j])))
}

## [1] TRUE
## [1] TRUE
## [1] TRUE

#Identificando pontos influentes por coeficiente
i_pi_diferenca = which(abs(vdiferencas)>2/sqrt(n), arr.ind = TRUE)
length(i_pi_diferenca)

## [1] 0

i_pi_dfbetas = which(abs(vdfbetas)>1, arr.ind = TRUE)
length(i_pi_dfbetas)

## [1] 0

if(length(i_pi_dfbetas)>0){
  for (i in 1:nrow(i_pi_dfbetas)){
      cat("Observação", i_pi_dfbetas[i,1],
      "| coeficiente:", i_pi_dfbetas[i,2], 
      "| y:", y[i_pi_dfbetas[i,1]], 
      "| x1:", x[i_pi_dfbetas[i,1],2], 
      "| x2:", x[i_pi_dfbetas[i,1],3], 
        "| e:", e[i_pi_dfbetas[i,1]], 
        "| prob:", round(pnorm(e[i_pi_dfbetas[i,1]],0,sqrt(sigma2)),3), "\n")
  }
}

#######################
# DFFITS
#######################

# Funcao para calcular DFFITS para uma observacao i
funcaodffits = function(i) {
  #Excluindo a i-C)sima observacao
  dados_menos_i = data.frame(y = dados$y[-i], 
                             x1 = dados$x1[-i], 
                             x2 = dados$x2[-i])
  
  # Reajustando o modelo sem a observacao i
  ajuste_menos_i    = lm(y~x1+x2, data=dados_menos_i)
  resultado_menos_i = summary(ajuste_menos_i)
  y_chapeu_i        = predict(ajuste_menos_i, newdata = dados[i,-1])
  
  # Calculando DFFITS
  vdffits = (ajuste$fitted.values[i] - y_chapeu_i)/(resultado_menos_i$sigma * sqrt(h_ii[i]))
  
  return(vdffits)
}

# Calculando DFFITS para todas as observacoes
vdffits = sapply(1:n, funcaodffits)
all.equal(vdffits, dffits(ajuste))

## [1] TRUE

#Identificando pontos influentes por coeficiente
i_pi_dffits = which(abs(vdffits)>2*sqrt((p+1)/n))
if(length(i_pi_dffits)>0){
  for (i in i_pi_dffits){
    cat("Observação", i,
        "| y:", y[i], 
        "| x1:", x[i,2], 
        "| x2:", x[i,3], 
        "| e:", e[i], 
        "| prob:", round(pnorm(e[i],0,sqrt(sigma2)),3), "\n")
  }
}

## Observação 26 | y: 4.372439 | x1: 0.8245212 | x2: 2.559005 | e: 2.134894 | prob: 0.984 
## Observação 37 | y: 7.31991 | x1: -1.572929 | x2: 0.1175274 | e: 2.508807 | prob: 0.994 
## Observação 48 | y: 20.78487 | x1: 0.400483 | x2: -2.915826 | e: -2.079881 | prob: 0.019 
## Observação 61 | y: -1.422379 | x1: 0.3082229 | x2: 3.38082 | e: 1.176234 | prob: 0.88 
## Observação 89 | y: 0.8064964 | x1: -2.491961 | x2: 0.7959731 | e: 1.466273 | prob: 0.929 
## Observação 90 | y: 4.839448 | x1: -2.247318 | x2: 0.1761963 | e: 2.286188 | prob: 0.989 
## Observação 150 | y: 14.05511 | x1: -0.7722187 | x2: -0.5139464 | e: 4.315985 | prob: 1 
## Observação 166 | y: 8.697254 | x1: -1.360454 | x2: -0.1926751 | e: 2.007916 | prob: 0.978 
## Observação 176 | y: 12.7589 | x1: -1.408311 | x2: -2.114393 | e: -1.473743 | prob: 0.07 
## Observação 182 | y: 7.563814 | x1: -2.201099 | x2: -0.4954539 | e: 2.185296 | prob: 0.986 
## Observação 204 | y: -3.744922 | x1: -1.319631 | x2: 1.951517 | e: -1.979962 | prob: 0.024 
## Observação 211 | y: -5.111897 | x1: -2.685922 | x2: 1.468242 | e: -1.181162 | prob: 0.119 
## Observação 219 | y: 6.582804 | x1: 0.5890755 | x2: 0.6399348 | e: -2.624684 | prob: 0.004 
## Observação 232 | y: 8.903186 | x1: -1.360999 | x2: -0.2105709 | e: 2.143898 | prob: 0.984 
## Observação 241 | y: -0.1837041 | x1: -1.622058 | x2: 0.9037209 | e: -1.702646 | prob: 0.044 
## Observação 271 | y: 1.832247 | x1: 1.633218 | x2: 2.884802 | e: -1.528198 | prob: 0.063 
## Observação 275 | y: 2.854965 | x1: -0.4271697 | x2: 0.7739318 | e: -2.767799 | prob: 0.003 
## Observação 280 | y: 8.619424 | x1: -1.81632 | x2: -0.5152574 | e: 2.007354 | prob: 0.978 
## Observação 287 | y: 8.754844 | x1: -1.62435 | x2: -0.457339 | e: 1.798536 | prob: 0.964 
## Observação 377 | y: 13.56357 | x1: 1.960508 | x2: 1.156053 | e: 2.306261 | prob: 0.989 
## Observação 396 | y: 20.06076 | x1: 0.7921559 | x2: -2.259573 | e: -1.354 | prob: 0.088 
## Observação 412 | y: 19.81558 | x1: 1.161328 | x2: -1.068847 | e: 2.056204 | prob: 0.98 
## Observação 413 | y: 13.64113 | x1: -0.9605667 | x2: -2.094473 | e: -1.855062 | prob: 0.032 
## Observação 431 | y: 19.38734 | x1: 0.4989811 | x2: -1.485232 | e: 1.949468 | prob: 0.974 
## Observação 470 | y: 6.474111 | x1: 0.9938529 | x2: 1.077408 | e: -2.197816 | prob: 0.014 
## Observação 492 | y: 7.122174 | x1: 0.2564287 | x2: 1.445259 | e: 2.133922 | prob: 0.984 
## Observação 525 | y: 1.133038 | x1: -1.294539 | x2: 0.7197132 | e: -2.104491 | prob: 0.018 
## Observação 530 | y: -0.5933575 | x1: -0.3113079 | x2: 2.763298 | e: 1.393759 | prob: 0.918 
## Observação 537 | y: 12.47597 | x1: -1.739338 | x2: -1.476201 | e: 1.789186 | prob: 0.963 
## Observação 542 | y: -2.82063 | x1: -2.684862 | x2: 0.859052 | e: -1.329838 | prob: 0.092 
## Observação 568 | y: 11.16764 | x1: -1.956702 | x2: -2.103578 | e: -1.376565 | prob: 0.084 
## Observação 579 | y: 20.13965 | x1: 1.310321 | x2: -1.077901 | e: 1.897089 | prob: 0.971 
## Observação 592 | y: 0.4399612 | x1: -1.827071 | x2: 0.3454315 | e: -2.6971 | prob: 0.003 
## Observação 593 | y: 5.172499 | x1: 1.206872 | x2: 2.424477 | e: 1.249791 | prob: 0.894 
## Observação 609 | y: 13.90946 | x1: -0.8052681 | x2: -1.038422 | e: 2.171575 | prob: 0.985 
## Observação 614 | y: -0.7347587 | x1: -1.505964 | x2: 1.087036 | e: -1.868724 | prob: 0.031 
## Observação 658 | y: 13.85024 | x1: -0.5977614 | x2: -1.921985 | e: -2.044419 | prob: 0.02 
## Observação 670 | y: 19.18354 | x1: 2.119951 | x2: -0.1879731 | e: 2.071797 | prob: 0.981 
## Observação 689 | y: 17.44954 | x1: -0.2463994 | x2: -1.483002 | e: 2.256729 | prob: 0.988 
## Observação 746 | y: 0.2085876 | x1: -0.8899715 | x2: 1.050535 | e: -2.919359 | prob: 0.002 
## Observação 789 | y: 7.382225 | x1: 1.450287 | x2: 2.171232 | e: 1.716294 | prob: 0.957 
## Observação 811 | y: 17.80154 | x1: 1.445334 | x2: -0.3645195 | e: 2.007461 | prob: 0.978 
## Observação 819 | y: 2.335637 | x1: 0.6657867 | x2: 1.647809 | e: -3.070488 | prob: 0.001 
## Observação 824 | y: 6.936204 | x1: 1.618409 | x2: 1.598688 | e: -1.524271 | prob: 0.064 
## Observação 825 | y: 11.67566 | x1: 2.204862 | x2: 0.7494605 | e: -1.94108 | prob: 0.026 
## Observação 843 | y: 22.21702 | x1: 3.003199 | x2: -1.087249 | e: -1.14157 | prob: 0.127 
## Observação 873 | y: 18.32846 | x1: -0.116331 | x2: -2.55633 | e: -1.54787 | prob: 0.061 
## Observação 909 | y: 13.78534 | x1: 0.5652866 | x2: -1.079497 | e: -2.228512 | prob: 0.013 
## Observação 919 | y: 17.27718 | x1: 1.025254 | x2: -0.454347 | e: 2.384031 | prob: 0.991 
## Observação 930 | y: -1.410353 | x1: -1.12791 | x2: 2.629543 | e: 2.491549 | prob: 0.994 
## Observação 962 | y: 2.406069 | x1: -1.446752 | x2: 1.330303 | e: 2.067538 | prob: 0.981

4 Covariáveis qualitativas

Variáveis categóricas representam grupos ou categorias, como:

Sexo: Masculino / Feminino
Tipo de Produto: A / B / C
Região: Norte / Sul / Leste / Oeste

Essas variáveis não possuem uma relação numérica natural entre os seus valores, por isso não podem ser incluídas diretamente em um modelo de regressão linear como covariáveis numéricas.

4.1 Como a regressão linear “interpreta” as covariáveis categóricas?

O modelo de regressão linear trata cada categoria como se representasse uma “reta de regressão diferente”, mas todas essas retas compartilham alguns parâmetros em comum.

Para isso, o modelo cria variáveis indicadoras (dummies), que funcionam como “interruptores” que ligam ou desligam o efeito de cada categoria. Dessa forma, a regressão consegue ajustar diferentes interceptos (ou inclinações, caso haja interação) para cada grupo, mantendo uma estrutura comum para os demais parâmetros.

4.2 Exemplo prático: salário, sexo e idade

Suponha que desejamos avaliar o salário de uma população. Para isso, usamos como covariáveis o sexo e a idade.

A variável sexo assume os valores “Feminino” e “Masculino” e pode ser codificada da seguinte forma:

sexoMasculino = 1, se o indivíduo for do sexo masculino;
sexoMasculino = 0, se for do sexo feminino.

Assim, o modelo de regressão ajusta uma reta para cada grupo de sexo, mas ambas podem compartilhar o mesmo coeficiente da idade e a mesma variância. Ou seja:

Quando o sexo for feminino (sexoMasculino = 0), teremos um modelo de regressão do salário em função da idade, com um determinado intercepto.
Quando o sexo for masculino (sexoMasculino = 1), o modelo ajustará um intercepto diferente, mas a variância e a inclinação (coeficiente da idade) podem ser as mesmas, caso não haja interação entre sexo e idade.

Isso permite ao modelo representar diferenças sistemáticas entre os grupos, sem perder a capacidade de interpretar o efeito geral da idade sobre o salário.

5 Conjunto de dados 2

5.1 Simulação

Simule um conjunto de dados do modelo de regressão linear da seguinte forma:

Para as 500 primeiras observações, considere que \(Y_i = 50 + 10 X_{i1} + e_i, \quad e_i \sim N(0, 4)\);
Para as 500 últimas observações, considere que \(Y_i = 70 + 10 X_{i1} + e_i, \quad e_i \sim N(0, 4)\).

Considere que \(x_{i1} \stackrel{iid}{\sim} N(0,1)\) para todas as observações.

# Configurações iniciais
n  = 1000  # Número de observações
# Simulando as covariáveis
x1 = rnorm(n, mean = 0, sd = 1)   # Covariável 1
x  = cbind(1, x1)

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

# Definindo os coeficientes reais
beta_0 = 70    # Intercepto do grupo 1
beta_1 = 10    # Coeficiente para x1
beta_2 = -20    # parte do Intercepto do grupo 1
beta   = c(beta_0, beta_1, beta_2)

# Gerando a variável resposta com erro aleatório
sigma2 = 4
e = rnorm(n, mean = 0, sd = sqrt(sigma2))  # Erro aleatório
y = NULL
y[1:500] = (beta_0+beta_2) + x1[1:500]*beta_1 + e[1:500]
y[501:n] = beta_0 + x1[501:n]*beta_1 + e[501:n]

rm(beta_0, beta_1, beta_2)

# Criando um data frame com os dados
dados   = data.frame(y, x1=x[,2])

5.2 Analisando os dados

Analise o conjunto de dados simulado.

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

round(cor(dados),3)

##        y    x1
## y  1.000 0.709
## x1 0.709 1.000

5.3 Ajustando um MRL aos dados (ajuste 1)

Ajuste o modelo \(Y_i = \beta_0 + \beta_1 X_{i1} + e_i\), \(e_i \sim N(0, \sigma^2)\). Em seguida, faça o gráfico de dispersão com a reta ajustada incluída e analise os resíduos.

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

plot(x1, y, bty="n")
lines(x1, ajuste$fitted.values, col="red", lwd=2)

#analisando os valores estimados
parametros = c("Est. pontual", "Lim Inf do IC", "Lim Sup do IC", "Valor esperado")# Nomes das colunas
aux   = round(cbind(ajuste$coefficients, confint(ajuste, level = 0.9), 
                    c((beta[1] + beta[1] + beta[3]) / 2, beta[2])),4)
print(rbind(parametros, aux))

##                            5 %             95 %                            
## parametros  "Est. pontual" "Lim Inf do IC" "Lim Sup do IC" "Valor esperado"
## (Intercept) "60.0139"      "59.478"        "60.5499"       "60"            
## x1          "10.1799"      "9.6518"        "10.708"        "10"

alfa = 0.1
aux = round(c(summary(ajuste)$sigma^2, 
  (n-p-1)*summary(ajuste)$sigma^2/qchisq(1-alfa/2, n-p-1), 
  (n-p-1)*summary(ajuste)$sigma^2/qchisq(alfa/2, n-p-1), sigma2),4)
print(rbind(parametros, aux))

##            [,1]           [,2]            [,3]            [,4]            
## parametros "Est. pontual" "Lim Inf do IC" "Lim Sup do IC" "Valor esperado"
## aux        "105.9663"     "98.5957"       "114.2465"      "4"

#analisando os resíduos
ks.test(ajuste$residuals,"pnorm",0,sd(ajuste$residuals))

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  ajuste$residuals
## D = 0.22758, p-value < 2.2e-16
## alternative hypothesis: two-sided

hist(ajuste$residuals, main="", freq=F, ylab="densidade", xlab="resíduos ordinários")
curve(dnorm(x, 0, sqrt(resultado$sigma^2)), add=T, col="red", lwd=2)

qqnorm(ajuste$residuals, xlab="quantis teóricos", ylab="quantis amostrais", bty="n", main="Resíduos ordinários")
qqline(ajuste$residuals, col = "red", lwd=2)

5.4 Ajustando um MRL aos dados (ajuste 2)

Ajuste o modelo \(Y_i = \beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2} + e_i\), tendo \(e_i \sim N(0, \sigma^2)\). Considere que \(X_{i2} = \left\{ \begin{array}{cc} 1, & i=1, \ldots, 500,\\ 0, & i=501, \ldots, 1000. \end{array}\right .\) Em seguida, faça o gráfico de dispersão com a reta ajustada incluída e analise os resíduos.

dados$x2 = c(rep(1,500), rep(0, 500))

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

#analisando os valores estimados
parametros = c("Est. pontual", "Lim Inf do IC", "Lim Sup do IC", "Valor verdadeiro")# Nomes das colunas
aux   = round(cbind(ajuste$coefficients, confint(ajuste, level = 0.9), beta),4)
print(rbind(parametros, aux))

##                            5 %             95 %            beta              
## parametros  "Est. pontual" "Lim Inf do IC" "Lim Sup do IC" "Valor verdadeiro"
## (Intercept) "70.0901"      "69.9382"       "70.2419"       "70"              
## x1          "10.0976"      "9.9918"        "10.2034"       "10"              
## x2          "-20.1516"     "-20.3663"      "-19.9369"      "-20"

alfa = 0.1
aux = round(c(summary(ajuste)$sigma^2, 
  (n-p-1)*summary(ajuste)$sigma^2/qchisq(1-alfa/2, n-p-1), 
  (n-p-1)*summary(ajuste)$sigma^2/qchisq(alfa/2, n-p-1), sigma2),4)
print(rbind(parametros, aux))

##            [,1]           [,2]            [,3]            [,4]              
## parametros "Est. pontual" "Lim Inf do IC" "Lim Sup do IC" "Valor verdadeiro"
## aux        "4.2522"       "3.9565"        "4.5845"        "4"

#analisando os ajustes
reta1 = ajuste$coefficients[1] + ajuste$coefficients[2]*x1
reta2 = reta1 + ajuste$coefficients[3]
plot(x1, y, bty="n")
lines(x1, reta1, col="red", lwd=2)
lines(x1, reta2, col="blue", lwd=2)

#analisando os resíduos
ks.test(ajuste$residuals,"pnorm",0,sd(ajuste$residuals))

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  ajuste$residuals
## D = 0.020329, p-value = 0.803
## alternative hypothesis: two-sided

hist(ajuste$residuals, main="", freq=F, ylab="densidade", xlab="resíduos ordinários")
curve(dnorm(x, 0, sqrt(resultado$sigma^2)), add=T, col="red", lwd=2)

qqnorm(ajuste$residuals, xlab="quantis teóricos", ylab="quantis amostrais", bty="n", main="Resíduos ordinários")
qqline(ajuste$residuals, col = "red", lwd=2)

5.5 Ajustando um MRL aos dados (ajuste 3)

Ajuste o modelo \(Y_i = \beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2} + \beta_3 X_{i1} X_{i2} + e_i\), tendo \(e_i \sim N(0, \sigma^2)\). Considere que \(X_{i2} = \left\{ \begin{array}{cc} 1, & i=1, \ldots, 500,\\ 0, & i=501, \ldots, 1000. \end{array}\right .\) Em seguida, faça o gráfico de dispersão com a reta ajustada incluída e analise os resíduos.

dados$x2 = c(rep(1,500), rep(0, 500))
dados$x1x2 = dados$x1*dados$x2

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

#analisando os valores estimados
parametros = c("Est. pontual", "Lim Inf do IC", "Lim Sup do IC", "Valor esperado")# Nomes das colunas
aux   = round(cbind(ajuste$coefficients, confint(ajuste, level = 0.9), 
                    c(beta,0)),4)
print(rbind(parametros, aux))

##                            5 %             95 %                            
## parametros  "Est. pontual" "Lim Inf do IC" "Lim Sup do IC" "Valor esperado"
## (Intercept) "70.0906"      "69.9387"       "70.2425"       "70"            
## x1          "10.0572"      "9.9029"        "10.2116"       "10"            
## x2          "-20.152"      "-20.3668"      "-19.9372"      "-20"           
## x1x2        "0.0762"       "-0.1359"       "0.2882"        "0"

alfa = 0.1
aux = round(c(summary(ajuste)$sigma^2, 
  (n-p-1)*summary(ajuste)$sigma^2/qchisq(1-alfa/2, n-p-1), 
  (n-p-1)*summary(ajuste)$sigma^2/qchisq(alfa/2, n-p-1), sigma2),4)
print(rbind(parametros, aux))

##            [,1]           [,2]            [,3]            [,4]            
## parametros "Est. pontual" "Lim Inf do IC" "Lim Sup do IC" "Valor esperado"
## aux        "4.255"        "3.959"         "4.5875"        "4"

#analisando os ajustes
reta1 = ajuste$coefficients[1] + ajuste$coefficients[2]*x1
reta2 = reta1 + ajuste$coefficients[3]
plot(x1, y, bty="n")
lines(x1, reta1, col="red", lwd=2)
lines(x1, reta2, col="blue", lwd=2)

#analisando os resíduos
ks.test(ajuste$residuals,"pnorm",0,sd(ajuste$residuals))

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  ajuste$residuals
## D = 0.020314, p-value = 0.8037
## alternative hypothesis: two-sided

hist(ajuste$residuals, main="", freq=F, ylab="densidade", xlab="resíduos ordinários")
curve(dnorm(x, 0, sqrt(resultado$sigma^2)), add=T, col="red", lwd=2)

qqnorm(ajuste$residuals, xlab="quantis teóricos", ylab="quantis amostrais", bty="n", main="Resíduos ordinários")
qqline(ajuste$residuals, col = "red", lwd=2)

6 Conjunto de dados 3

6.1 Simulação

Simule um conjunto de dados do modelo de regressão linear da seguinte forma:

Para as 500 primeiras observações, considere que \(Y_i = 50 + 10 X_{i1} + e_i, \quad e_i \sim N(0, 4)\);
Para as 500 últimas observações, considere que \(Y_i = 70 + 30 X_{i1} + e_i, \quad e_i \sim N(0, 4)\).

Considere que \(X_{i1} \stackrel{iid}{\sim} N(0,1)\) para todas as observações.

# Configurações iniciais
n  = 1000  # Número de observações
# Simulando as covariáveis
x1 = rnorm(n, mean = 0, sd = 1)   # Covariável 1
x  = cbind(1, x1)

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

# Definindo os coeficientes reais
beta_0 = 70    # Intercepto do grupo 0
beta_1 = 30    # Coeficiente para x1 do grupo 0
beta_2 = -20   # parte do intercepto do grupo 1
beta_3 = -20   # parte do Coeficiente para x1 do grupo 1
beta   = c(beta_0, beta_1, beta_2, beta_3)

# Gerando a variável resposta com erro aleatório
sigma2 = 4
e = rnorm(n, mean = 0, sd = sqrt(sigma2))  # Erro aleatório
#y = x %*% beta + e
y[1:500] = beta_0 + beta_2 + x1[1:500]*(beta_1+beta_3) + e[1:500]
y[501:n] = beta_0 + x1[501:n]*beta_1 + e[501:n]

rm(beta_0, beta_1, beta_2, beta_3)

# Criando um data frame com os dados
dados   = data.frame(y, x1=x[,2])

6.2 Analisando os dados

Analise o conjunto de dados simulado.

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

round(cor(dados),3)

##       y   x1
## y  1.00 0.82
## x1 0.82 1.00

6.3 Ajustando um MRL aos dados (ajuste 1)

Ajuste o modelo \(Y_i = \beta_0 + \beta_1 X_{i1} + e_i\), \(e_i \sim N(0, \sigma^2)\). Em seguida, faça o gráfico de dispersão com a reta ajustada incluída e analise os resíduos.

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

#analisando os valores estimados
parametros = c("Est. pontual", "Lim Inf do IC", "Lim Sup do IC", "Valor esperado")# Nomes das colunas
aux   = round(cbind(ajuste$coefficients, confint(ajuste, level = 0.9), 
                    c((beta[1]+beta[1]+beta[3])/2, (beta[2]+beta[2]+beta[4])/2)),4)
print(rbind(parametros, aux))

##                            5 %             95 %                            
## parametros  "Est. pontual" "Lim Inf do IC" "Lim Sup do IC" "Valor esperado"
## (Intercept) "60.3083"      "59.5753"       "61.0413"       "60"            
## x1          "20.0859"      "19.3552"       "20.8167"       "20"

alfa = 0.1
aux = round(c(summary(ajuste)$sigma^2, 
  (n-p-1)*summary(ajuste)$sigma^2/qchisq(1-alfa/2, n-p-1), 
  (n-p-1)*summary(ajuste)$sigma^2/qchisq(alfa/2, n-p-1), sigma2),4)
print(rbind(parametros, aux))

##            [,1]           [,2]            [,3]            [,4]            
## parametros "Est. pontual" "Lim Inf do IC" "Lim Sup do IC" "Valor esperado"
## aux        "197.989"      "184.2178"      "213.46"        "4"

plot(x1, y, bty="n")
lines(x1, ajuste$fitted.values, col="red", lwd=2)

#analisando os resíduos
ks.test(ajuste$residuals,"pnorm",0,sd(ajuste$residuals))

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  ajuste$residuals
## D = 0.025174, p-value = 0.5505
## alternative hypothesis: two-sided

hist(ajuste$residuals, main="", freq=F, ylab="densidade", xlab="resíduos ordinários")
curve(dnorm(x, 0, sqrt(resultado$sigma^2)), add=T, col="red", lwd=2)

qqnorm(ajuste$residuals, xlab="quantis teóricos", ylab="quantis amostrais", bty="n", main="Resíduos ordinários")
qqline(ajuste$residuals, col = "red", lwd=2)

6.4 Ajustando um MRL aos dados (ajuste 2)

Ajuste o modelo \(Y_i = \beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2} + \beta_3 X_{i1} X_{i2} + e_i\), tendo \(e_i \sim N(0, \sigma^2)\). Considere que \(X_{i2} = \left\{ \begin{array}{cc} 1, & i=1, \ldots, 500,\\ 0, & i=501, \ldots, 1000. \end{array}\right .\) Em seguida, faça o gráfico de dispersão com a reta ajustada incluída e analise os resíduos.

dados$x2 = c(rep(1,500), rep(0, 500))

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


#analisando os valores estimados
parametros = c("Est. pontual", "Lim Inf do IC", "Lim Sup do IC", "Valor esperado")# Nomes das colunas
aux   = round(cbind(ajuste$coefficients, confint(ajuste, level = 0.9), beta),4)
print(rbind(parametros, aux))

##                            5 %             95 %            beta            
## parametros  "Est. pontual" "Lim Inf do IC" "Lim Sup do IC" "Valor esperado"
## (Intercept) "69.9176"      "69.7746"       "70.0607"       "70"            
## x1          "30.0016"      "29.8573"       "30.1458"       "30"            
## x2          "-19.9238"     "-20.1263"      "-19.7213"      "-20"           
## x1:x2       "-20.0886"     "-20.2905"      "-19.8866"      "-20"

alfa = 0.1
aux = round(c(summary(ajuste)$sigma^2, 
  (n-p-1)*summary(ajuste)$sigma^2/qchisq(1-alfa/2, n-p-1), 
  (n-p-1)*summary(ajuste)$sigma^2/qchisq(alfa/2, n-p-1), sigma2),4)
print(rbind(parametros, aux))

##            [,1]           [,2]            [,3]            [,4]            
## parametros "Est. pontual" "Lim Inf do IC" "Lim Sup do IC" "Valor esperado"
## aux        "3.7747"       "3.5122"        "4.0697"        "4"

reta1 = ajuste$coefficients[1] + ajuste$coefficients[2]*x1
reta2 = reta1 + ajuste$coefficients[3] + ajuste$coefficients[4]*x1
plot(x1, y, bty="n")
lines(x1, reta1, col="red", lwd=2)
lines(x1, reta2, col="blue", lwd=2)

#analisando os resíduos
ks.test(ajuste$residuals,"pnorm",0,sd(ajuste$residuals))

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  ajuste$residuals
## D = 0.017412, p-value = 0.9222
## alternative hypothesis: two-sided

hist(ajuste$residuals, main="", freq=F, ylab="densidade", xlab="resíduos ordinários")
curve(dnorm(x, 0, sqrt(resultado$sigma^2)), add=T, col="red", lwd=2)

qqnorm(ajuste$residuals, xlab="quantis teóricos", ylab="quantis amostrais", bty="n", main="Resíduos ordinários")
qqline(ajuste$residuals, col = "red", lwd=2)

7 Conjunto de dados do R

Vamos usar o dataset mtcars para ilustrar.

data(mtcars)
mtcars$am <- factor(mtcars$am, labels = c("Automático", "Manual"))
head(mtcars)

##                    mpg cyl disp  hp drat    wt  qsec vs         am gear carb
## Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0     Manual    4    4
## Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0     Manual    4    4
## Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1     Manual    4    1
## Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1 Automático    3    1
## Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0 Automático    3    2
## Valiant           18.1   6  225 105 2.76 3.460 20.22  1 Automático    3    1

7.1 Ajustando o modelo

modelo <- lm(mpg ~ wt + am, data = mtcars)
summary(modelo)

## 
## Call:
## lm(formula = mpg ~ wt + am, data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.5295 -2.3619 -0.1317  1.4025  6.8782 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 37.32155    3.05464  12.218 5.84e-13 ***
## wt          -5.35281    0.78824  -6.791 1.87e-07 ***
## amManual    -0.02362    1.54565  -0.015    0.988    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.098 on 29 degrees of freedom
## Multiple R-squared:  0.7528, Adjusted R-squared:  0.7358 
## F-statistic: 44.17 on 2 and 29 DF,  p-value: 1.579e-09

7.1.1 Interpretação

wt: A cada 1000 libras a mais, espera-se uma redução de β1 mpg.
amManual: Carros manuais, em média, fazem β2 mpg a mais que automáticos, para o mesmo peso.

7.2 Como o R codifica variáveis categóricas?

contrasts(mtcars$am)

##            Manual
## Automático      0
## Manual          1

7.3 Visualização com `ggplot2`

library(ggplot2)

ggplot(mtcars, aes(x = wt, y = mpg, color = am)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "Consumo vs Peso por Tipo de Transmissão",
       x = "Peso (1000 lbs)", y = "Milhas por galão (mpg)") +
  theme_minimal()

8 Exercício

Use o dataset iris para ajustar um modelo de regressão com:

Sepal.Length como variável resposta
Sepal.Width como preditora quantitativa
Species como preditora categórica

Qual é o grupo de referência de Species?
Qual o efeito de Sepal.Width?
Qual a diferença entre as espécies?

Aula de Modelos Lineares I

Patrícia Lusié Velozo da Costa

2025-06-13

Modelo de Regressão Linear Múltipla

1 Simulando conjunto de dados

2 Ajustando um MRLM

3 Pontos a serem analisados

3.1 Identificando outliers

3.2 Identificando pontos de alavancagem

3.3 Identificando pontos influentes

4 Covariáveis qualitativas

4.1 Como a regressão linear “interpreta” as covariáveis categóricas?

4.2 Exemplo prático: salário, sexo e idade

5 Conjunto de dados 2

5.1 Simulação

5.2 Analisando os dados

5.3 Ajustando um MRL aos dados (ajuste 1)

5.4 Ajustando um MRL aos dados (ajuste 2)

5.5 Ajustando um MRL aos dados (ajuste 3)

6 Conjunto de dados 3

6.1 Simulação

6.2 Analisando os dados

6.3 Ajustando um MRL aos dados (ajuste 1)

6.4 Ajustando um MRL aos dados (ajuste 2)

7 Conjunto de dados do R

7.1 Ajustando o modelo

7.1.1 Interpretação

7.2 Como o R codifica variáveis categóricas?

7.3 Visualização com `ggplot2`

8 Exercício

Aula de Modelos Lineares I

Patrícia Lusié Velozo da Costa

2025-06-13

Modelo de Regressão Linear Múltipla

1 Simulando conjunto de dados

2 Ajustando um MRLM

3 Pontos a serem analisados

3.1 Identificando outliers

3.2 Identificando pontos de alavancagem

3.3 Identificando pontos influentes

4 Covariáveis qualitativas

4.1 Como a regressão linear “interpreta” as covariáveis categóricas?

4.2 Exemplo prático: salário, sexo e idade

5 Conjunto de dados 2

5.1 Simulação

5.2 Analisando os dados

5.3 Ajustando um MRL aos dados (ajuste 1)

5.4 Ajustando um MRL aos dados (ajuste 2)

5.5 Ajustando um MRL aos dados (ajuste 3)

6 Conjunto de dados 3

6.1 Simulação

6.2 Analisando os dados

6.3 Ajustando um MRL aos dados (ajuste 1)

6.4 Ajustando um MRL aos dados (ajuste 2)

7 Conjunto de dados do R

7.1 Ajustando o modelo

7.1.1 Interpretação

7.2 Como o R codifica variáveis categóricas?

7.3 Visualização com ggplot2

8 Exercício

7.3 Visualização com `ggplot2`