Aula de Modelos Lineares I
Patrícia Lusié Velozo da Costa
2025-06-18
Modelo de Regressão Linear Múltipla - Medidas Corretivas para Linearidade
1 Simulando um conjunto de dados
- Simule um conjunto de dados considerando n=1000 (tamanho amostral), \(\beta_0=2\), \(\beta_1=5\), \(\beta_2=3\), \(\sigma^2=1\) e que \(x_{i1} \stackrel{iid}{\sim} N(0,1)\), \(x_{i2} \stackrel{iid}{\sim} N(0, 1)\) e o seguinte modelo \[Y_i \stackrel{indep}{\sim} N(\beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2}^2, \sigma^2).\]
# Configurações iniciais
n = 1000 # Número de observações
# Simulando as covariáveis
x1 = rnorm(n, mean = 0, sd = sqrt(1)) # Covariável 1
x2 = rnorm(n, mean = 0, sd = sqrt(1)) # Covariável 2
x2m = x2^2 # Covariável 2
x = cbind(1, x1=x1, x2=x2)
xm = cbind(1, x1=x1, x2=x2m)
#definindo total de covariaveis
p = ncol(x)-1
# Definindo os coeficientes reais
beta_0 = 2 # Intercepto
beta_1 = 5 # Coeficiente para x1
beta_2 = 3 # Coeficiente para x2
beta = c(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 = xm %*% beta + e
# Criando um data frame com os dados
dados = data.frame(y, x1=x1, x2=x2)
dadosm = data.frame(y, x1=x1, x2=x2m)
rm(beta_0, beta_1, beta_2)
#cor(dados)
#cor(dadosm)2 Analisando os dados
- Analise o conjunto de dados simulado.
#analisando os dados gerados
par(mfrow=c(1,2))
hist(y, main="", freq=F, bty="n", ylab="densidade")
hist(e, main="", freq=F, bty="n", ylab="densidade")
par(mfrow=c(1,3))
plot(x1, y, bty="n")
plot(x2, y, bty="n")
plot(x2m, y, bty="n")
par(mfrow=c(1,2))
plot(x1, x2, bty="n")
plot(x1, x2m, bty="n")
round(cor(dados),3)## y x1 x2
## y 1.000 0.776 0.013
## x1 0.776 1.000 0.050
## x2 0.013 0.050 1.000round(cor(dadosm),3)## y x1 x2
## y 1.000 0.776 0.667
## x1 0.776 1.000 0.074
## x2 0.667 0.074 1.0003 Ajustando o MRL aos dados
- Ajuste os seguintes modelos:
\[\mbox{Ajuste 1: }Y_i \stackrel{indep}{\sim} N(\beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2}, \sigma^2).\] \[\mbox{Ajuste 2: }Y_i \stackrel{indep}{\sim} N(\beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2}^2, \sigma^2).\] \[\mbox{Ajuste 3: }Y_i \stackrel{indep}{\sim} N(\beta_0 + \beta_2 X_{i2}, \sigma^2).\]
- Ajuste um modelo de regressão linear múltipla usando as covariáveis \(X_1\) e \(X_2\).
#ajustando o modelo
ajuste1 = lm(y ~ x1 + x2, data = dados)
summary(ajuste1)##
## Call:
## lm(formula = y ~ x1 + x2, data = dados)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.196 -2.689 -1.313 1.238 27.286
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.1164 0.1346 38.012 <2e-16 ***
## x1 5.2857 0.1357 38.957 <2e-16 ***
## x2 -0.1735 0.1315 -1.319 0.188
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.253 on 997 degrees of freedom
## Multiple R-squared: 0.6036, Adjusted R-squared: 0.6028
## F-statistic: 759 on 2 and 997 DF, p-value: < 2.2e-16#ajustando o modelo correto
ajuste2 = lm(y ~ x1 + x2, data = dadosm)
summary(ajuste2)##
## Call:
## lm(formula = y ~ x1 + x2, data = dadosm)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.1479 -0.6738 -0.0166 0.6879 2.9146
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.00742 0.04081 49.18 <2e-16 ***
## x1 4.96880 0.03295 150.82 <2e-16 ***
## x2 2.96191 0.02342 126.46 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.031 on 997 degrees of freedom
## Multiple R-squared: 0.9767, Adjusted R-squared: 0.9767
## F-statistic: 2.089e+04 on 2 and 997 DF, p-value: < 2.2e-16#ajustando o modelo com somente x2
ajuste3 = lm(y ~ x2, data = dados)
summary(ajuste3)##
## Call:
## lm(formula = y ~ x2, data = dados)
##
## Residuals:
## Min 1Q Median 3Q Max
## -16.462 -4.507 -0.591 3.757 34.413
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.19617 0.21363 24.323 <2e-16 ***
## x2 0.08387 0.20855 0.402 0.688
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.751 on 998 degrees of freedom
## Multiple R-squared: 0.000162, Adjusted R-squared: -0.0008398
## F-statistic: 0.1617 on 1 and 998 DF, p-value: 0.68774 Seleção de Modelos
#através de R2
c(summary(ajuste1)$r.squared, summary(ajuste2)$r.squared, summary(ajuste3)$r.squared)## [1] 0.6035921913 0.9766969658 0.0001620301#através de R2 ajustado
c(summary(ajuste1)$adj.r.squared, summary(ajuste2)$adj.r.squared, summary(ajuste3)$adj.r.squared)## [1] 0.6027969901 0.9766502194 -0.0008398116#através de AIC
c(AIC(ajuste1), AIC(ajuste2), AIC(ajuste3))## [1] 5738.250 2904.390 6661.399#através de BIC
c(BIC(ajuste1), BIC(ajuste2), BIC(ajuste3))## [1] 5757.881 2924.021 6676.1235 Analisando os resíduos do modelo 3
hist(ajuste3$residuals, main="", freq=F, ylab="densidade", xlab="resíduos")
curve(dnorm(x, 0, sqrt(summary(ajuste3)$sigma^2)), add=T, col="red", lwd=2)
qqnorm(ajuste3$residuals, xlab="quantis teóricos", ylab="quantis amostrais", bty="n")
qqline(ajuste3$residuals, col="red", lwd=2)
shapiro.test(ajuste3$residuals)##
## Shapiro-Wilk normality test
##
## data: ajuste3$residuals
## W = 0.97486, p-value = 3.867e-12ks.test(ajuste3$residuals,"pnorm",0, sqrt(summary(ajuste3)$sigma^2))##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: ajuste3$residuals
## D = 0.053412, p-value = 0.006655
## alternative hypothesis: two-sidedpar(mfrow=c(1,3))
plot(x2, ajuste3$residuals)
abline(h=0, lwd=2, col="red", lty=3)
plot(y, ajuste3$residuals)
abline(h=0, lwd=2, col="red", lty=3)
plot(ajuste3$residuals)
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(ajuste3)##
## studentized Breusch-Pagan test
##
## data: ajuste3
## BP = 4.9319, df = 1, p-value = 0.02636#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.9609, p-value = 0.2686
## 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.6 Analisando os resíduos do modelo 2
hist(ajuste2$residuals, main="", freq=F, ylab="densidade", xlab="resíduos")
curve(dnorm(x, 0, sqrt(summary(ajuste2)$sigma^2)), add=T, col="red", lwd=2)
qqnorm(ajuste2$residuals, xlab="quantis teóricos", ylab="quantis amostrais", bty="n")
qqline(ajuste2$residuals, col="red", lwd=2)
shapiro.test(ajuste2$residuals)##
## Shapiro-Wilk normality test
##
## data: ajuste2$residuals
## W = 0.99855, p-value = 0.5845ks.test(ajuste2$residuals,"pnorm",0, sqrt(summary(ajuste2)$sigma^2))##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: ajuste2$residuals
## D = 0.015639, p-value = 0.9673
## alternative hypothesis: two-sidedpar(mfrow=c(1,4))
plot(x1, ajuste2$residuals)
abline(h=0, lwd=2, col="red", lty=3)
plot(x2, ajuste2$residuals)
abline(h=0, lwd=2, col="red", lty=3)
plot(y, ajuste2$residuals)
abline(h=0, lwd=2, col="red", lty=3)
plot(ajuste2$residuals)
abline(h=0, lwd=2, col="red", lty=3)
# Testando Homocedasticidade (Teste de Breusch-Pagan)
bptest(ajuste2)##
## studentized Breusch-Pagan test
##
## data: ajuste2
## BP = 1.6183, df = 2, p-value = 0.4452#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 = 2.0192, p-value = 0.6193
## 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.