Aula de Modelos Lineares I
Patrícia Lusié Velozo da Costa
2025-05-16
Modelo de Regressão Linear Múltipla
Nesta aula computacional, analisaremos um conjunto de dados proveniente de um 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.\]
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=2\), \(\beta_2=-1,5\), \(\beta_3=3\), \(\beta_4=-1\), \(\sigma^2=9\) e que \(X_{i1} \stackrel{iid}{\sim} N(0,1)\), \(X_{i2} \stackrel{iid}{\sim} Unif(-\sqrt{3}, \sqrt{3})\), \(X_{i3} \stackrel{iid}{\sim} t_5\) e , \(X_{i4} \stackrel{iid}{\sim} Exp(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 = runif(n, -sqrt(3), sqrt(3)) # Covariável 2
x3 = rt(n, df = 5) # Covariável 3
x4 = rexp(n, rate = 1) # Covariável 4
x = cbind(1, x1, x2, x3, x4)
rm(x1, x2, x3, x4)
# Definindo os coeficientes reais
beta_0 = 10 # Intercepto
beta_1 = 2 # Coeficiente para x1
beta_2 = -1.5 # Coeficiente para x2
beta_3 = 3 # Coeficiente para x3
beta_4 = -1 # Coeficiente para x4
beta = c(beta_0, beta_1, beta_2, beta_3, beta_4)
rm(beta_0, beta_1, beta_2, beta_3, beta_4)
# Gerando a variável resposta com erro aleatório
sigma2 = 9
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=x[,2], x2=x[,3], x3=x[,4], x4=x[,5])2 Analisando os dados gerados
Calcule a correlação entre a variável resposta e cada uma das covariáveis e entre cada par de covariáveis. Faça gráficos de dispersão da variável resposta com cada covariável.
# Correlação entre as variáveis
round(cor(dados),3)## y x1 x2 x3 x4
## y 1.000 0.350 -0.308 0.644 -0.176
## x1 0.350 1.000 0.010 -0.007 -0.008
## x2 -0.308 0.010 1.000 -0.025 -0.015
## x3 0.644 -0.007 -0.025 1.000 0.034
## x4 -0.176 -0.008 -0.015 0.034 1.000#Analisando as variaveis selecionadas
correlacao <- cor(dados)
#correlacao <- round(cor(dados), 3)
#as.data.frame(correlacao)
#representando graficamente
#install.packages("corrplot")
library(corrplot)
par(mfrow=c(1,1))
corrplot(correlacao, method = "color",
type = "upper",
tl.col = "black",
tl.srt = 45,
addCoef.col = "black") # Adiciona os valores
par(mfrow=c(2,2))
plot(dados$x1, dados$y, ylab="y", bty="n", xlab=expression(x[1]))
plot(dados$x2, dados$y, ylab="y", bty="n", xlab=expression(x[2]))
plot(dados$x3, dados$y, ylab="y", bty="n", xlab=expression(x[3]))
plot(dados$x4, dados$y, ylab="y", bty="n", xlab=expression(x[4]))
3 Estimando pontualmente os coeficientes de um MRLM
Suponha que \(\boldsymbol{\beta} = (\beta_0, \beta_1, \ldots, \beta_p)^T\) e \(\sigma^2\) são desconhecidos. Considere que os dados seguem um modelo de regressão linear múltipla com todas as \(p=4\) covariáveis geradas.
Obtenha as estimativas pontuais destes parâmetros usando os estimadores não tendenciosos vistos em sala de aula: \[\hat{\boldsymbol{\beta}} = (\boldsymbol{X}^T\boldsymbol{X})^{-1}\boldsymbol{X}^T\boldsymbol{Y}\] \[\hat{\sigma}^2 = \frac{\hat{\boldsymbol{e}}^T\hat{\boldsymbol{e}}}{n-p-1},\] sendo \(\hat{\boldsymbol{e}} = \boldsymbol{Y} - \boldsymbol{X}\hat{\boldsymbol{\beta}}\).
betachapeu = solve(t(x)%*%x)%*%t(x)%*%y
echapeu = y - x%*%betachapeu
p = ncol(x)-1
sigma2chapeu = t(echapeu)%*%echapeu/(n-p-1)4 Estimando pontualmente os coeficientes usando lm de um MRLM
Obtenha as estimativas pontuais destes parâmetros usando a função lm e compare com os encontrados anteriormente.
ajuste = lm(y~x1+x2+x3+x4, data=dados)
#ajuste = lm(y~x1+x2+x3+x4)
resultado = summary(ajuste)
cbind(betachapeu, ajuste$coefficients)## [,1] [,2]
## 10.346403 10.346403
## x1 1.940737 1.940737
## x2 -1.615854 -1.615854
## x3 2.948995 2.948995
## x4 -1.139961 -1.139961c(sigma2chapeu, resultado$sigma^2)## [1] 9.530079 9.5300795 Encontrando a matriz de covariância do betachapeu
Vimos que \(\hat{\boldsymbol{\beta}} \sim N_{p+1}(\boldsymbol{\beta}, \sigma^2 \left(\boldsymbol{X}^T\boldsymbol{X})^{-1}\right)\).
Encontre essa matriz de covariância substituindo a variância pela estimativa encontrada.
Em seguida, para cada estimador \(\hat{\beta}_j\), calcule seu respectivo desvio padrão e compare com os valores encontrados pelo summary do lm.
#calculando "manualmente"
var_beta0chapeu = resultado$sigma^2 * solve(t(x)%*%x)
round(var_beta0chapeu,3)## x1 x2 x3 x4
## 0.019 -0.001 0.00 0.000 -0.010
## x1 -0.001 0.010 0.00 0.000 0.000
## x2 0.000 0.000 0.01 0.000 0.000
## x3 0.000 0.000 0.00 0.007 0.000
## x4 -0.010 0.000 0.00 0.000 0.011#comparando com o encontrado usando as funções lm e summary
cbind(sqrt(diag(var_beta0chapeu)), resultado$coefficients[,2])## [,1] [,2]
## 0.13746044 0.13746044
## x1 0.09942185 0.09942185
## x2 0.09860356 0.09860356
## x3 0.08339254 0.08339254
## x4 0.10415831 0.10415831#calculando a correlacao
cor_beta0chapeu = matrix(NA, p+1, p+1)
for(i in 1:(p+1)){
for(j in 1:(p+1)){
cor_beta0chapeu[i,j] = var_beta0chapeu[i,j] / sqrt(var_beta0chapeu[i,i]*var_beta0chapeu[j,j])
}
}
round(cor_beta0chapeu,3)## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.000 -0.039 0.019 0.012 -0.703
## [2,] -0.039 1.000 -0.010 0.007 0.008
## [3,] 0.019 -0.010 1.000 0.024 0.015
## [4,] 0.012 0.007 0.024 1.000 -0.033
## [5,] -0.703 0.008 0.015 -0.033 1.0006 Considerando somente a 1ª covariável e ajustando um MRLS
Ajuste um modelo de regressão linear simples considerando apenas a primeira covariável. Analise o modelo ajustado.
OBS: Não altere o dado simulado. Apenas o modelo ajustado.
#Ajustando o modelo
ajuste1 = lm(y~x1, data=dados)
resultado1 = summary(ajuste1)
#IC para beta0 e beta1 usando funcao do R
gama = 0.95
confint(ajuste1, level = gama)## 2.5 % 97.5 %
## (Intercept) 9.104468 9.727500
## x1 1.589684 2.223446# Intervalo de confiança para sigma^2
alfa = 1-gama
p = 1
c((n-p-1) * resultado1$sigma^2 / qchisq(1-alfa/2, n-p-1),
(n-p-1) * resultado1$sigma^2 / qchisq(alfa/2, n-p-1))## [1] 23.07767 27.50676residuos = ajuste1$residuals
hist(residuos, freq=F, main="Histograma", ylab="densidade")
curve(dnorm(x, 0, sqrt(resultado1$sigma^2)), add=T, lwd=2, col="red")
qqnorm(residuos, xlab="quantis teoricos", ylab="quantis amostrais", bty="n")
qqline(residuos, col="red", lwd=2)
#Interpretação: No gráfico Q-Q, os resíduos devem seguir aproximadamente a linha reta se forem normalmente distribuídos.
shapiro.test(residuos)##
## Shapiro-Wilk normality test
##
## data: residuos
## W = 0.99546, p-value = 0.004616ks.test(residuos,"pnorm",mean(residuos), sd(residuos))##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: residuos
## D = 0.024147, p-value = 0.6044
## alternative hypothesis: two-sided#Um p-valor pequeno indica que os resíduos não seguem uma distribuição normal.
# Testando Homocedasticidade (Teste de Breusch-Pagan)
# Instalar pacote necessário (se ainda não instalado)
# install.packages("lmtest")
library(lmtest)
bptest(ajuste1)##
## studentized Breusch-Pagan test
##
## data: ajuste1
## BP = 0.080958, df = 1, p-value = 0.776#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(ajuste1)##
## Durbin-Watson test
##
## data: ajuste1
## DW = 1.9502, p-value = 0.2149
## 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.
#Analisando os residuos
par(mfrow=c(2,2))
plot(dados$x1, residuos)
plot(ajuste1$fitted.values, residuos)
plot(residuos)
plot(dados$x1, dados$y, main = "y ~ x1", xlab = "x1", ylab = "y")
abline(ajuste1, col = "red", lwd = 2)
7 Considerando somente a 2ª covariável e ajustando um MRLS
Ajuste um modelo de regressão linear simples considerando apenas a segunda covariável. Analise o modelo ajustado.
OBS: Não altere o dado simulado. Apenas o modelo ajustado.
# Correlação entre as variáveis
#round(cor(dados),3)
#Ajustando o modelo
ajuste2 = lm(y~x2, data=dados)
resultado2 = summary(ajuste2)
#IC para beta0 e beta1 usando funcao do R
gama = 0.95
confint(ajuste2, level = gama)## 2.5 % 97.5 %
## (Intercept) 9.118697 9.751248
## x2 -1.984305 -1.346234# Intervalo de confiança para sigma^2
alfa = 1-gama
p = 1
c((n-p-1) * resultado2$sigma^2 / qchisq(1-alfa/2, n-p-1),
(n-p-1) * resultado2$sigma^2 / qchisq(alfa/2, n-p-1))## [1] 23.79921 28.36678residuos = ajuste2$residuals
hist(residuos, freq=F, main="Histograma", ylab="densidade")
curve(dnorm(x, 0, sqrt(resultado2$sigma^2)), add=T, lwd=2, col="red")
qqnorm(residuos, xlab="quantis teoricos", ylab="quantis amostrais", bty="n")
qqline(residuos, col="red", lwd=2)
#Interpretação: No gráfico Q-Q, os resíduos devem seguir aproximadamente a linha reta se forem normalmente distribuídos.
shapiro.test(residuos)##
## Shapiro-Wilk normality test
##
## data: residuos
## W = 0.99455, p-value = 0.001115ks.test(residuos,"pnorm",mean(residuos), sd(residuos))##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: residuos
## D = 0.02537, p-value = 0.5404
## alternative hypothesis: two-sided#Um p-valor pequeno indica que os resíduos não seguem uma distribuição normal.
# Testando Homocedasticidade (Teste de Breusch-Pagan)
# Instalar pacote necessário (se ainda não instalado)
# install.packages("lmtest")
library(lmtest)
bptest(ajuste2)##
## studentized Breusch-Pagan test
##
## data: ajuste2
## BP = 0.20556, df = 1, p-value = 0.6503#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(ajuste2)##
## Durbin-Watson test
##
## data: ajuste2
## DW = 1.8946, p-value = 0.04746
## 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.
#Analisando os residuos
par(mfrow=c(2,2))
plot(dados$x2, residuos)
plot(ajuste2$fitted.values, residuos)
plot(residuos)
plot(dados$x2, dados$y, main = "y ~ x2", xlab = "x2", ylab = "y")
abline(ajuste2, col = "red", lwd = 2)
8 Considerando somente a 3ª covariável e ajustando um MRLS
Ajuste um modelo de regressão linear simples considerando apenas a terceira covariável. Analise o modelo ajustado.
OBS: Não altere o dado simulado. Apenas o modelo ajustado.
# Correlação entre as variáveis
#round(cor(dados),3)
#Ajustando o modelo
ajuste3 = lm(y~x3, data=dados)
resultado3 = summary(ajuste3)
#IC para beta0 e beta1 usando funcao do R
gama = 0.95
confint(ajuste3, level = gama)## 2.5 % 97.5 %
## (Intercept) 9.189717 9.698025
## x3 2.723024 3.156793# Intervalo de confiança para sigma^2
alfa = 1-gama
p = 1
c((n-p-1) * resultado3$sigma^2 / qchisq(1-alfa/2, n-p-1),
(n-p-1) * resultado3$sigma^2 / qchisq(alfa/2, n-p-1))## [1] 15.39001 18.34368residuos = ajuste3$residuals
hist(residuos, freq=F, main="Histograma", ylab="densidade")
curve(dnorm(x, 0, sqrt(resultado3$sigma^2)), add=T, lwd=2, col="red")
qqnorm(residuos, xlab="quantis teoricos", ylab="quantis amostrais", bty="n")
qqline(residuos, col="red", lwd=2)
#Interpretação: No gráfico Q-Q, os resíduos devem seguir aproximadamente a linha reta se forem normalmente distribuídos.
shapiro.test(residuos)##
## Shapiro-Wilk normality test
##
## data: residuos
## W = 0.99834, p-value = 0.4535ks.test(residuos,"pnorm",mean(residuos), sd(residuos))##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: residuos
## D = 0.017962, p-value = 0.9036
## alternative hypothesis: two-sided#Um p-valor pequeno indica que os resíduos não seguem uma distribuição normal.
# Testando Homocedasticidade (Teste de Breusch-Pagan)
# Instalar pacote necessário (se ainda não instalado)
# install.packages("lmtest")
library(lmtest)
bptest(ajuste3)##
## studentized Breusch-Pagan test
##
## data: ajuste3
## BP = 0.94616, df = 1, p-value = 0.3307#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(ajuste3)##
## Durbin-Watson test
##
## data: ajuste3
## DW = 1.9329, p-value = 0.1441
## 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.
#Analisando os residuos
par(mfrow=c(2,2))
plot(dados$x3, residuos)
plot(ajuste3$fitted.values, residuos)
plot(residuos)
plot(dados$x3, dados$y, main = "y ~ x3", xlab = "x3", ylab = "y")
abline(ajuste3, col = "red", lwd = 2)
9 Considerando somente a 4ª covariável e ajustando um MRLS
Ajuste um modelo de regressão linear simples considerando apenas a quarta covariável. Analise o modelo ajustado.
OBS: Não altere o dado simulado. Apenas o modelo ajustado.
# Correlação entre as variáveis
#round(cor(dados),3)
#Ajustando o modelo
ajuste4 = lm(y~x4, data=dados)
resultado4 = summary(ajuste4)
#IC para beta0 e beta1 usando funcao do R
gama = 0.95
p = 1
confint(ajuste4, level = gama)## 2.5 % 97.5 %
## (Intercept) 9.976473 10.8964235
## x4 -1.354376 -0.6570941# Intervalo de confiança para sigma^2
alfa = 1-gama
c((n-p-1) * resultado4$sigma^2 / qchisq(1-alfa/2, n-p-1),
(n-p-1) * resultado4$sigma^2 / qchisq(alfa/2, n-p-1))## [1] 25.48289 30.37360residuos = ajuste4$residuals
hist(residuos, freq=F, main="Histograma", ylab="densidade")
curve(dnorm(x, 0, sqrt(resultado4$sigma^2)), add=T, lwd=2, col="red")
qqnorm(residuos, xlab="quantis teoricos", ylab="quantis amostrais", bty="n")
qqline(residuos, col="red", lwd=2)
#Interpretação: No gráfico Q-Q, os resíduos devem seguir aproximadamente a linha reta se forem normalmente distribuídos.
shapiro.test(residuos)##
## Shapiro-Wilk normality test
##
## data: residuos
## W = 0.99682, p-value = 0.04229ks.test(residuos,"pnorm",mean(residuos), sd(residuos))##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: residuos
## D = 0.027609, p-value = 0.431
## alternative hypothesis: two-sided#Um p-valor pequeno indica que os resíduos não seguem uma distribuição normal.
# Testando Homocedasticidade (Teste de Breusch-Pagan)
# Instalar pacote necessário (se ainda não instalado)
# install.packages("lmtest")
library(lmtest)
bptest(ajuste4)##
## studentized Breusch-Pagan test
##
## data: ajuste4
## BP = 1.9395, df = 1, p-value = 0.1637#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(ajuste4)##
## Durbin-Watson test
##
## data: ajuste4
## DW = 1.9426, p-value = 0.1819
## 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.
#Analisando os residuos
par(mfrow=c(2,2))
plot(dados$x4, residuos)
plot(ajuste4$fitted.values, residuos)
plot(residuos)
plot(dados$x4, dados$y, main = "y ~ x4", xlab = "x4", ylab = "y")
abline(ajuste4, col = "red", lwd = 2)
10 Ajuste com parte das covariáveis
Ajuste um modelo de regressão linear múltipla para os dados gerados considerando apenas 2 covariáveis. O que podemos dizer sobre os resíduos destes ajustes? E se fizermos um modelo usando 3 covariáveis destas?
OBS: Não altere o dado simulado. Apenas o modelo ajustado.
11 Ajuste com outras covariáveis
Gere uma quinta covariável da seguinte forma: \[X_{i5} \sim N(0,1), \quad i=1, \ldots, n.\]
Ajuste o seguinte modelo de regressão linear múltipla para os dados gerados na primeira tarefa: \[Y_i = \beta_0 + \beta_1 X_{i1} + \ldots + \beta_5 X_{i5} + e_i, \quad e_i \stackrel{iid}{\sim} N(0,\sigma^2), \quad i=1,\ldots,n.\]
OBS: Não altere o dado simulado. Apenas o modelo ajustado.
12 Ajuste com covariáveis dependentes
Gere uma sexta covariável da seguinte forma: \[X_{i5} = 3 X_{i2} + \nu_i, \quad \nu_i \stackrel{iid}{\sim} N(0, 1/2), \quad i=1, \ldots, n.\]
Ajuste o seguinte modelo de regressão linear múltipla para os dados gerados na primeira tarefa: \[Y_i = \beta_0 + \beta_1 X_{i1} + \ldots + \beta_4 X_{i4} + \beta_6 X_{i6} + e_i, \quad e_i \stackrel{iid}{\sim} N(0,\sigma^2), \quad i=1,\ldots,n.\]
OBS: Não altere o dado simulado. Apenas o modelo ajustado.