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.
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.\]
Simule um conjunto de dados do modelo de regressão linear múltipla considerando n=1000 (tamanho amostral) e \(n_p = 10\) (tamanho da previsão), \(\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
no = 1000 # Número de observações
np = 10
n = no+np
# 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
#dividindo a amostra em observacoes e previsao
n = no
xp = x[(n+1):(n+np),]
x = x[1:n,]
yp = y[(n+1):(n+np)]
y = y[1:n]
# Criando um data frame com os dados
dados = data.frame(y, x1=x[,2], x2=x[,3], x3=x[,4], x4=x[,5])Considere \((\beta_0, \beta_1, \ldots, \beta_p)\) desconhecidos mas \(\sigma^2\) conhecido.
Lembrete:
\[\frac{\hat{\beta}_j-\beta_j}{\sqrt{\sigma^2 C_{jj}}} \sim N\left(0 , 1\right),\] sendo \(C_{jj}\) o \(j\)-ésimo elemento da diagonal principal da matriz \((\boldsymbol{X}^T\boldsymbol{X})^{-1}\).
#ajustando o modelo
ajuste = lm(y~x1+x2+x3+x4, data = dados)
resultado = summary(ajuste)
#obtendo as estimativas pontuais
sigma2chapeu = resultado$sigma^2
betachapeu = ajuste$coefficients
#criando os ICs
gama = 0.95 #nivel de confianca
alfa = 1-gama #nivel de significancia
p = ncol(x)-1 #definindo o total de covariaveis
IC1 = matrix(NA, p+1, 2) #criando matriz para guardar os ICs
dC = diag(solve(t(x)%*%x))
erro1 = qnorm(1-alfa/2)*sqrt(sigma2chapeu * dC)
IC1[,1] = betachapeu - erro1
IC1[,2] = betachapeu + erro1
IC1## [,1] [,2]
## [1,] 9.502960 10.0189154
## [2,] 1.922251 2.2852736
## [3,] -1.710952 -1.3577282
## [4,] 2.830168 3.0998082
## [5,] -1.018767 -0.6313613Os valores verdadeiros dos coeficientes estão incluídos no IC criado acima?
for(j in 1:(p+1)){
if((beta[j]>=IC1[j,1]) & (beta[j]<=IC1[j,2])){
print(paste0("Beta_",j," pertence ao IC criado"))}else{
print(paste0("Beta_",j," NAO pertence ao IC criado"))}
}## [1] "Beta_1 pertence ao IC criado"
## [1] "Beta_2 pertence ao IC criado"
## [1] "Beta_3 pertence ao IC criado"
## [1] "Beta_4 pertence ao IC criado"
## [1] "Beta_5 pertence ao IC criado"Considere \((\beta_0, \beta_1, \ldots, \beta_p, \sigma^2)\) desconhecidos.
Crie um intervalo de confiança bilateral para cada parâmetro \(\beta_j\), \(j=0, 1, \ldots, p\), considerando um nível de confiança de 95%. Mas NÃO USE a função confint.
Crie um intervalo de confiança bilateral para cada parâmetro \(\beta_j\), \(j=0, 1, \ldots, p\), considerando um nível de confiança de 95%. Mas agora USE a função confint.
Os intervalos encontrados coincidiram?
_ Compare os intervalos quando \(\sigma^2\) é conhecido com os intervalos obtidos quando \(\sigma^2\) é desconhecido.
Lembrete:
\[\frac{\hat{\beta}_j-\beta_j}{\sqrt{\hat{\sigma}^2 C_{jj}}} \sim t_{(n-p-1)},\] sendo \(C_{jj}\) o \(j\)-ésimo elemento da diagonal principal da matriz \((\boldsymbol{X}^T\boldsymbol{X})^{-1}\).
#criando os ICs
IC = matrix(NA, p+1, 2) #criando matriz para guardar os ICs
erro = qt(1-alfa/2, n-p-1)*sqrt(sigma2chapeu * dC)
IC[,1] = betachapeu - erro
IC[,2] = betachapeu + erro
#comparando os ICs calculados com os obtidos pela funcao confint
cbind(round(IC,3),"",round(confint(ajuste),3))## 2.5 % 97.5 %
## (Intercept) "9.503" "10.019" "" "9.503" "10.019"
## x1 "1.922" "2.285" "" "1.922" "2.285"
## x2 "-1.711" "-1.358" "" "-1.711" "-1.358"
## x3 "2.83" "3.1" "" "2.83" "3.1"
## x4 "-1.019" "-0.631" "" "-1.019" "-0.631"#comparando os ICs de sigma2 conhecido versus desconhecido
cbind(erro1, erro)## erro1 erro
## 0.2579779 0.2582921
## x1 0.1815113 0.1817324
## x2 0.1766118 0.1768269
## x3 0.1348200 0.1349842
## x4 0.1937028 0.1939388Os valores verdadeiros dos coeficientes estão incluídos no IC criado acima?
for(j in 1:(p+1)){
if((beta[j]>=IC[j,1]) & (beta[j]<=IC[j,2])){
print(paste0("Beta_",j," pertence ao IC criado"))}else{
print(paste0("Beta_",j," NAO pertence ao IC criado"))}
}## [1] "Beta_1 pertence ao IC criado"
## [1] "Beta_2 pertence ao IC criado"
## [1] "Beta_3 pertence ao IC criado"
## [1] "Beta_4 pertence ao IC criado"
## [1] "Beta_5 pertence ao IC criado"Crie um intervalo de confiança para a variância \(\sigma^2\).
Lembrete:
\[\frac{(n-p-1) \hat{\sigma}^2}{\sigma^2} \sim \chi^2_{(n-p-1)}.\]
#criando o IC
q1 = qchisq(alfa/2, n-p-1)
q2 = qchisq(1-alfa/2, n-p-1)
linf_IC = (n-p-1)*sigma2chapeu/q2
lsup_IC = (n-p-1)*sigma2chapeu/q1
c(linf_IC, lsup_IC)## [1] 7.932114 9.456959O valor verdadeiro da variância populacional está incluído no IC criado acima?
if((sigma2>=linf_IC) & (sigma2<=lsup_IC)){
print("Sigma2 pertence ao IC criado")}else{
print("Sigma2 NAO pertence ao IC criado")}## [1] "Sigma2 pertence ao IC criado"Vimos que podemos testar as hipóteses \[H_0: \beta_j = 0 \mbox{ versus } H_0: \beta_j \neq 0, \quad j=0, \ldots, p,\] de 3 formas diferentes.
Quais seriam as suas conclusões com base nos ICs criados anteriormente?
b = 0
for(j in 1:(p+1)){
if((b>=IC[j,1]) & (b<=IC[j,2])){
print(paste0("NAO rejeito H0: beta_",j-1," = ",b, " a favor de H1: beta_",j-1," <> ",b))}else{
print(paste0("Rejeito H0: beta_",j-1," = ",b, " a favor de H1: beta_",j-1," <> ",b))}
}## [1] "Rejeito H0: beta_0 = 0 a favor de H1: beta_0 <> 0"
## [1] "Rejeito H0: beta_1 = 0 a favor de H1: beta_1 <> 0"
## [1] "Rejeito H0: beta_2 = 0 a favor de H1: beta_2 <> 0"
## [1] "Rejeito H0: beta_3 = 0 a favor de H1: beta_3 <> 0"
## [1] "Rejeito H0: beta_4 = 0 a favor de H1: beta_4 <> 0"Quais seriam as suas conclusões com base nas Estatísticas de Teste apropriadas?
Lembrando que para o teste de hipóteses \[H_0: \beta_j = b \mbox{ versus } H_0: \beta_j \neq b, \quad j=0, \ldots, p,\] assumindo \(H_0\) verdadeiro, temos que a Estatística de Teste é dada por:
\(\frac{\hat{\beta}_j - b}{\sqrt{\sigma^2 C_{jj}}} \sim N(0,1)\), quando \(\sigma^2\) for conhecido
\(\frac{\hat{\beta}_j - b}{\sqrt{\hat{\sigma}^2 C_{jj}}} \sim t_{(n-p-1)}\), quando \(\sigma^2\) for desconhecido
sendo \(C_{jj}\) o \(j\)-ésimo elemento da diagonal principal de \((\boldsymbol{X}^T\boldsymbol{X})^{-1}\)
#Supondo sigma2 conhecido
ET_betas = betachapeu / sqrt(sigma2*dC)
z = qnorm(1-alfa/2)
for(j in 1:(p+1)){
if((ET_betas[j]>=(-z)) & (ET_betas[j]<=z)){print(paste0("NAO rejeito H0: beta_",j-1," = ",b, " a favor de H1: beta_",j-1," <> ",b))}else{print(paste0("Rejeito H0: beta_",j-1," = ",b, " a favor de H1: beta_",j-1," <> ",b))}
}## [1] "Rejeito H0: beta_0 = 0 a favor de H1: beta_0 <> 0"
## [1] "Rejeito H0: beta_1 = 0 a favor de H1: beta_1 <> 0"
## [1] "Rejeito H0: beta_2 = 0 a favor de H1: beta_2 <> 0"
## [1] "Rejeito H0: beta_3 = 0 a favor de H1: beta_3 <> 0"
## [1] "Rejeito H0: beta_4 = 0 a favor de H1: beta_4 <> 0"#Supondo sigma2 desconhecido
ET2_betas = betachapeu / sqrt(sigma2chapeu*dC)
t = qt(1-alfa/2, n-p-1)
for(j in 1:(p+1)){
if((ET2_betas[j]>=(-t)) & (ET2_betas[j]<=t)){print(paste0("NAO rejeito H0: beta_",j-1," = ",b, " a favor de H1: beta_",j-1," <> ",b))}else{print(paste0("Rejeito H0: beta_",j-1," = ",b, " a favor de H1: beta_",j-1," <> ",b))}
}## [1] "Rejeito H0: beta_0 = 0 a favor de H1: beta_0 <> 0"
## [1] "Rejeito H0: beta_1 = 0 a favor de H1: beta_1 <> 0"
## [1] "Rejeito H0: beta_2 = 0 a favor de H1: beta_2 <> 0"
## [1] "Rejeito H0: beta_3 = 0 a favor de H1: beta_3 <> 0"
## [1] "Rejeito H0: beta_4 = 0 a favor de H1: beta_4 <> 0" #comparando ET2_betas com a saida do summary do lm
resultado##
## Call:
## lm(formula = y ~ x1 + x2 + x3 + x4, data = dados)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.1159 -1.9968 -0.0684 1.9211 8.6458
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.76094 0.13162 74.158 < 2e-16 ***
## x1 2.10376 0.09261 22.716 < 2e-16 ***
## x2 -1.53434 0.09011 -17.027 < 2e-16 ***
## x3 2.96499 0.06879 43.104 < 2e-16 ***
## x4 -0.82506 0.09883 -8.348 2.3e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.94 on 995 degrees of freedom
## Multiple R-squared: 0.7445, Adjusted R-squared: 0.7434
## F-statistic: 724.7 on 4 and 995 DF, p-value: < 2.2e-16 cbind(ET2_betas, resultado$coefficients[,3])## ET2_betas
## (Intercept) 74.157843 74.157843
## x1 22.716479 22.716479
## x2 -17.027467 -17.027467
## x3 43.103920 43.103920
## x4 -8.348334 -8.348334Quais seriam as suas conclusões com base nos p-valores apropriados?
#Supondo sigma2 conhecido
pvalor_betas = NULL
for(j in 1:(p+1)){
pvalor_betas[j] = 2*pnorm(abs(ET_betas[j]), lower.tail = F)
if(pvalor_betas[j] > alfa){print(paste0("NAO rejeito H0: beta_",j-1," = ",b, " a favor de H1: beta_",j-1," <> ",b))}else{print(paste0("Rejeito H0: beta_",j-1," = ",b, " a favor de H1: beta_",j-1," <> ",b))}
}## [1] "Rejeito H0: beta_0 = 0 a favor de H1: beta_0 <> 0"
## [1] "Rejeito H0: beta_1 = 0 a favor de H1: beta_1 <> 0"
## [1] "Rejeito H0: beta_2 = 0 a favor de H1: beta_2 <> 0"
## [1] "Rejeito H0: beta_3 = 0 a favor de H1: beta_3 <> 0"
## [1] "Rejeito H0: beta_4 = 0 a favor de H1: beta_4 <> 0"#Supondo sigma2 desconhecido
pvalor2_betas = NULL
for(j in 1:(p+1)){
pvalor2_betas[j] = 2*pt(abs(ET2_betas[j]), n-p-1, lower.tail = F)
if(pvalor2_betas[j] > alfa){print(paste0("NAO rejeito H0: beta_",j-1," = ",b, " a favor de H1: beta_",j-1," <> ",b))}else{print(paste0("Rejeito H0: beta_",j-1," = ",b, " a favor de H1: beta_",j-1," <> ",b))}
}## [1] "Rejeito H0: beta_0 = 0 a favor de H1: beta_0 <> 0"
## [1] "Rejeito H0: beta_1 = 0 a favor de H1: beta_1 <> 0"
## [1] "Rejeito H0: beta_2 = 0 a favor de H1: beta_2 <> 0"
## [1] "Rejeito H0: beta_3 = 0 a favor de H1: beta_3 <> 0"
## [1] "Rejeito H0: beta_4 = 0 a favor de H1: beta_4 <> 0" #comparando pvalor com a saida do summary do lm
cbind(pvalor2_betas, resultado$coefficients[,4])## pvalor2_betas
## (Intercept) 0.000000e+00 0.000000e+00
## x1 2.305672e-92 2.305672e-92
## x2 2.974355e-57 2.974355e-57
## x3 8.018825e-230 8.018825e-230
## x4 2.298401e-16 2.298401e-16Teste agora se \[H_0: \beta_0 = 8 \mbox{ versus } H_0: \beta_0 \neq 8.\]
Teste agora se \[H_0: \sigma^2 = 15 \mbox{ versus } H_0: \sigma^2 \neq 15.\]