Aula de Modelos Lineares I
Patrícia Lusié Velozo da Costa
2025-05-30
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=2\), \(\beta_2=-1,5\), \(\beta_3=3\), \(\beta_4=0\), \(\sigma^2=1\) 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)
# 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 = 0 # 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 = 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=x[,2], x2=x[,3], x3=x[,4], x4=x[,5])2 Analisando o conjunto de dados
Calcule a matriz de correlação entre a variável resposta e as covariáveis do conjunto de dados.
Faça os histogramas da variável resposta \(Y\) e de cada uma das covariáveis \(X_j\).
Faça o gráfico de dispersão de cada covariável \(X_j\) versus a variável resposta \(Y\). E faça o gráfico de dispersão da média \(\mu\) versus a variável resposta \(Y\).
# Visualizando as primeiras linhas do conjunto de dados
head(dados)## y x1 x2 x3 x4
## 1 4.978343 -1.2071812 1.68976842 0.2987788 0.33008955
## 2 3.383324 0.6516967 0.07218059 -1.9912121 0.11480027
## 3 5.635815 0.1625368 -0.81525718 -1.9683002 2.44663539
## 4 7.927642 -0.6671540 -1.08085016 -0.4377765 0.24808677
## 5 14.561758 -0.5280685 -1.17726044 1.7177771 0.07799962
## 6 11.201325 -1.1076671 -1.19475207 -0.1005727 0.56333544# Análise descritiva dos dados
# Estatísticas resumo
summary(dados)## y x1 x2 x3
## Min. :-11.516 Min. :-4.19090 Min. :-1.72937 Min. :-7.28617
## 1st Qu.: 7.112 1st Qu.:-0.70009 1st Qu.:-0.89873 1st Qu.:-0.68643
## Median : 10.197 Median :-0.05635 Median :-0.04506 Median : 0.06827
## Mean : 9.991 Mean :-0.05647 Mean :-0.02305 Mean : 0.03708
## 3rd Qu.: 13.165 3rd Qu.: 0.60968 3rd Qu.: 0.84307 3rd Qu.: 0.79114
## Max. : 24.544 Max. : 3.14792 Max. : 1.73194 Max. : 4.87593
## x4
## Min. :0.001057
## 1st Qu.:0.291541
## Median :0.685929
## Mean :0.999521
## 3rd Qu.:1.320853
## Max. :7.094201# Correlação entre as variáveis
round(cor(dados),3)## y x1 x2 x3 x4
## y 1.000 0.454 -0.269 0.827 0.016
## x1 0.454 1.000 0.015 0.056 0.014
## x2 -0.269 0.015 1.000 0.056 0.012
## x3 0.827 0.056 0.056 1.000 0.021
## x4 0.016 0.014 0.012 0.021 1.000# Visualização gráfica
# Histograma para cada variável
par(mfrow = c(2, 3)) # Organizando o layout de gráficos
hist(dados$y, main = "Histograma de Y", xlab = "y", col = "skyblue", border = "white", freq=F, ylab="densidade")
hist(dados$x1, main = "Histograma de X1", xlab = "x1", col = "skyblue", border = "white", freq=F, ylab="densidade")
hist(dados$x2, main = "Histograma de X2", xlab = "x2", col = "skyblue", border = "white", freq=F, ylab="densidade")
hist(dados$x3, main = "Histograma de X3", xlab = "x3", col = "skyblue", border = "white", freq=F, ylab="densidade")
hist(dados$x4, main = "Histograma de X4", xlab = "x4", col = "skyblue", border = "white", freq=F, ylab="densidade")
# Gráficos de dispersão para relação entre y e cada covariável
par(mfrow = c(2, 2)) # Ajustando layout para scatter plots
plot(dados$x1, dados$y, main = "y vs x1", xlab = "x1", ylab = "y", col = "blue", pch = 19)
plot(dados$x2, dados$y, main = "y vs x2", xlab = "x2", ylab = "y", col = "blue", pch = 19)
plot(dados$x3, dados$y, main = "y vs x3", xlab = "x3", ylab = "y", col = "blue", pch = 19)
plot(dados$x4, dados$y, main = "y vs x4", xlab = "x4", ylab = "y", col = "blue", pch = 19)
#plot(x%*%beta, dados$y, main = "y vs mu", xlab = expression(mu), ylab = "y", col = "blue", pch = 19)
# Gráficos de dispersão para relação entre mu=E[Y] e cada covariável
par(mfrow = c(2, 2))
plot(dados$x1, x%*%beta, main = "E[Y] vs x1", xlab = "x1", ylab = "E[Y]", col = "blue", pch = 19)
plot(dados$x2, x%*%beta, main = "E[Y] vs x2", xlab = "x2", ylab = "E[Y]", col = "blue", pch = 19)
plot(dados$x3, x%*%beta, main = "E[Y] vs x3", xlab = "x3", ylab = "E[Y]", col = "blue", pch = 19)
plot(dados$x4, x%*%beta, main = "E[Y] vs x4", xlab = "x4", ylab = "E[Y]", col = "blue", pch = 19)
3 Ajustando um MRLM
Ajuste um modelo de regressão linear múltipla e obtenha as estimativas pontuais e intervalares para \((\beta_0, \ldots, \beta_4, \sigma^2)\) considerando um nível de confiança de 90%.
#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
#definindo nivel de significancia
alfa = 0.1
#visualizando a estimativa pontual e a intervalar de beta
round(cbind(betachapeu, confint(ajuste, level = 1-alfa), beta),3) #obtendo as estimativas intervalares dos betas## betachapeu 5 % 95 % beta
## (Intercept) 9.974 9.900 10.049 10.0
## x1 2.004 1.949 2.058 2.0
## x2 -1.525 -1.578 -1.472 -1.5
## x3 3.031 2.990 3.073 3.0
## x4 -0.018 -0.070 0.035 0.0#definindo total de covariaveis
p = ncol(x)-1
#criando o IC de sigma2
q1 = qchisq(alfa/2, n-p-1)
q2 = qchisq(1-alfa/2, n-p-1)
ICs = c((n-p-1)*sigma2chapeu/q2, (n-p-1)*sigma2chapeu/q1)
#visualizando a estimativa pontual e a intervalar de sigma2
c(round(sigma2chapeu,3), paste0("[",round(ICs[1],3)," , ",round(ICs[2],3),"]"), sigma2)## [1] "1.019" "[0.948 , 1.099]" "1"4 IC para \(\mu_i\)
4.1 Considerando a variância populacional conhecida
Considerando a variância populacional conhecida, crie um intervalo de confiança para cada média populacional \(\mu_i = E[Y_i]\), \(i=1, \ldots, n\) e faça um gráfico para representar esses intervalos.
Lembrete:
\[\frac{\hat{Y}_i - \mu_i}{ \sqrt{\sigma^2 \boldsymbol{x}_i^T(\boldsymbol{X}^T\boldsymbol{X})^{-1}\boldsymbol{x}_i}} \sim N(0,1)\]
#criando os ICs
q = qnorm(1-alfa/2)
linf1 = NULL
lsup1 = NULL
ychapeu = ajuste$fitted.values
M = solve(t(x)%*%x)
for(i in 1:n){
erro = q* sqrt(sigma2 * x[i,]%*%M%*%x[i,])
linf1[i] = ychapeu[i] - erro
lsup1[i] = ychapeu[i] + erro
}
# Criando um data frame com as médias verdadeiras e os ICs
df_pred <- data.frame(
y = y,
Media_Verdadeira = x %*% beta, # Valores verdadeiros de E(Yi)
Predicao = ychapeu,
Limite_Inferior = linf1,
Limite_Superior = lsup1
)
# Gráfico das médias verdadeiras com os ICs
library(ggplot2)
inds = seq(1,10,by=1) #para facilitar a visualizacao do grafico
ggplot(df_pred[inds,], aes(x = inds)) +
geom_point(aes(y = y, color = "Observações"), size = 2, shape = 4, stroke = 1) +
geom_point(aes(y = Media_Verdadeira, color = "Média Verdadeira"), size = 2) +
geom_point(aes(y = Predicao, color = "Predição"), size = 2) +
geom_ribbon(aes(ymin = Limite_Inferior, ymax = Limite_Superior, fill = "Intervalo de Confiança"), alpha = 0.4) +
scale_color_manual(values = c("Observações" = "red", "Média Verdadeira" = "blue", "Predição" = "black")) +
scale_fill_manual(values = "gray") +
labs(
title = "Médias Verdadeiras com Intervalos de Confiança para as Médias",
x = "Observações",
y = "Média de Y",
color = "Legenda", # Título da legenda das linhas
fill = "Legenda IC" # Título da legenda da faixa de IC
) +
theme_minimal()
4.2 Considerando a variância populacional desconhecida
Considerando a variância populacional desconhecida, crie um intervalo de confiança para cada média populacional \(\mu_i = E[Y_i]\), \(i=1, \ldots, n\) e faça um gráfico para representar esses intervalos.
Lembrete:
\[\frac{\hat{Y}_i - \mu_i}{\sqrt{\hat{\sigma}^2 \boldsymbol{x}_i^T(\boldsymbol{X}^T\boldsymbol{X})^{-1}\boldsymbol{x}_i}} \sim t_{(n-p-1)}.\]
#criando os ICs
q = qt(1-alfa/2, n-p-1)
linf = NULL
lsup = NULL
for(i in 1:n){
erro = q* sqrt(sigma2chapeu * x[i,]%*%M%*%x[i,])
linf[i] = ychapeu[i] - erro
lsup[i] = ychapeu[i] + erro
}
# Calculando as predições e ICs para as médias
predicoes = predict(ajuste, interval = "confidence", level = 1 - alfa)
#comparando os valores calculados
summary(round(predicoes[,1] - ychapeu,5))## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 0 0 0 0summary(round(predicoes[,2] - linf,5))## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 0 0 0 0summary(round(predicoes[,3] - lsup,5))## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 0 0 0 0# Criando um data frame com as médias verdadeiras e os ICs
df_pred <- data.frame(
y = y,
Media_Verdadeira = x %*% beta, # Valores verdadeiros de E(Yi)
Predicao = predicoes[, "fit"],
Limite_Inferior = predicoes[, "lwr"],
Limite_Superior = predicoes[, "upr"],
Limite_Inferior1 = linf1,
Limite_Superior1 = lsup1
)
# Gráfico das médias verdadeiras com os ICs
library(ggplot2)
inds = seq(1,10,by=1) #para facilitar a visualizacao do grafico
ggplot(df_pred[inds, ], aes(x = inds)) +
geom_point(aes(y = y, color = "Observações"), size = 2, shape = 4, stroke = 1) +
geom_point(aes(y = Media_Verdadeira, color = "Média Verdadeira"), size = 2) +
geom_point(aes(y = Predicao, color = "Predição"), size = 2) +
geom_ribbon(aes(ymin = Limite_Inferior, ymax = Limite_Superior, fill = "Intervalo de Confiança"), alpha = 0.4) +
geom_ribbon(aes(ymin = Limite_Inferior1, ymax = Limite_Superior1, fill = "Intervalo de Confiança 1"), alpha = 0.4) +
# Bordas do intervalo 1
geom_line(aes(y = Limite_Inferior, color = "Limite IC"), linetype = "dotted") +
geom_line(aes(y = Limite_Superior, color = "Limite IC"), linetype = "dotted") +
scale_color_manual(values = c("Observações" = "red", "Média Verdadeira" = "blue", "Predição" = "black")) +
scale_fill_manual(values = c("Intervalo de Confiança" = "gray", "Intervalo de Confiança 1" = "yellow")) +
labs(
title = "Médias Verdadeiras com Intervalos de Confiança para as Médias",
x = "Observações",
y = "Média de Y",
color = "Legenda", # Título da legenda das linhas
fill = "Legenda IC" # Título da legenda da faixa de IC
) +
theme_minimal()
5 IC para \(Y_i\)
Considerando a variância populacional desconhecida, crie um intervalo de predição para cada \(Y_i\), \(i=1, \ldots, n\) e faça um gráfico para representar esses intervalos.
Lembrete:
\[\frac{Y_i - \hat{Y}_i}{ \sqrt{\sigma^2 \left( 1+ \boldsymbol{x}_i^T(\boldsymbol{X}^T\boldsymbol{X})^{-1}\boldsymbol{x}_i\right)}} \sim N(0,1)\] \[\frac{Y_i - \hat{Y}_i}{ \sqrt{\hat{\sigma}^2 \left( 1+ \boldsymbol{x}_i^T(\boldsymbol{X}^T\boldsymbol{X})^{-1}\boldsymbol{x}_i\right)}} \sim t_{(n-p-1)}.\]
#criando os ICs
q = qt(1-alfa/2, n-p-1)
linfy = NULL
lsupy = NULL
ychapeu = ajuste$fitted.values
M = solve(t(x)%*%x)
for(i in 1:n){
erro = q* sqrt(sigma2chapeu * (1+x[i,]%*%M%*%x[i,]))
linfy[i] = ychapeu[i] - erro
lsupy[i] = ychapeu[i] + erro
}
# Calculando as predições e ICs para as médias
predicoes = predict(ajuste, interval = "prediction", level = 1 - alfa)
#comparando os valores calculados
summary(round(predicoes[,1] - ychapeu,5))## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 0 0 0 0summary(round(predicoes[,2] - linfy,5))## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 0 0 0 0summary(round(predicoes[,3] - lsupy,5))## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 0 0 0 0# Criando um data frame com as médias verdadeiras e os ICs
df_pred <- data.frame(
y = y,
Predicao = predicoes[, "fit"],
Limite_Inferior = predicoes[, "lwr"],
Limite_Superior = predicoes[, "upr"]
)
# Gráfico das médias verdadeiras com os ICs
library(ggplot2)
inds = seq(1,30,by=1) #para facilitar a visualizacao do grafico
ggplot(df_pred[inds,], aes(x = inds)) +
geom_point(aes(y = y, color = "Observações"), size = 2, shape = 4, stroke = 1) +
geom_point(aes(y = Predicao, color = "Predição"), size = 2) +
geom_ribbon(aes(ymin = Limite_Inferior, ymax = Limite_Superior, fill = "Intervalo de Confiança"), alpha = 0.2) +
scale_color_manual(values = c("Observações" = "red", "Predição" = "black")) +
scale_fill_manual(values = "gray") +
labs(
title = "Observações com Intervalos de Predição para as Médias",
x = "Observações",
y = "Y",
color = "Legenda", # Título da legenda das linhas
fill = "Legenda IC" # Título da legenda da faixa de IC
) +
theme_minimal()
6 IC para \(Y_{\text{novo}}\) (interpolação versus extrapolação)
A formulação apresentada anteriormente pode ser utilizada tanto para construir intervalos de confiança para as unidades observadas quanto para prever o valor da variável resposta em unidades não observadas. Neste último caso, ao considerar uma nova unidade com vetor de covariáveis \(\boldsymbol{X}_{\text{novo}}\), temos:
\[ \frac{Y_{\text{novo}} - \hat{Y}_{\text{novo}}}{ \sqrt{\sigma^2 \left( 1 + \boldsymbol{X}_{\text{novo}}^T(\boldsymbol{X}^T\boldsymbol{X})^{-1} \boldsymbol{X}_{\text{novo}} \right)}} \sim N(0,1) \]
e
\[ \frac{Y_{\text{novo}} - \hat{Y}_{\text{novo}}}{ \sqrt{\hat{\sigma}^2 \left( 1 + \boldsymbol{X}_{\text{novo}}^T(\boldsymbol{X}^T\boldsymbol{X})^{-1} \boldsymbol{X}_{\text{novo}} \right)}} \sim t_{(n - p - 1)}. \]
Dizemos que estamos interpolando quando a nova unidade possui um valor de alavancagem inferior ou igual ao máximo valor de alavancagem das unidades observadas, ou seja, quando:
\[ h_{\text{novo}} \leq \max(h_i), \quad i = 1, \ldots, n \]
onde \(h_{\text{novo}} = \boldsymbol{X}_{\text{novo}}^T (\boldsymbol{X}^T \boldsymbol{X})^{-1} \boldsymbol{X}_{\text{novo}}\) é o valor de alavancagem da nova unidade e \(h_i\) são os valores de alavancagem das unidades observadas, obtidos a partir da matriz chapéu \(\boldsymbol{H} = \boldsymbol{X}(\boldsymbol{X}^T\boldsymbol{X})^{-1}\boldsymbol{X}^T\).
Por outro lado, chamamos de extrapolação quando:
\[ h_{\text{novo}} > \max(h_i) \]
Essa situação tende a aumentar a incerteza associada à predição e pode comprometer a validade das inferências baseadas no modelo ajustado.
- Preveja o valor da variável resposta quando \[\boldsymbol{x_i}^T=(1 \quad 0 \quad 0 \quad 0 \quad 0 ).\] Este caso é um caso de extrapolação ou interpolação? Justifique.
H = x%*%solve(t(x)%*%x)%*%t(x)
hmax = max(H)
xnovo = c(1 , 0, 0, 0, 0 )
hnovo = xnovo%*%solve(t(x)%*%x)%*%xnovo
if(hnovo > hmax){print("Extrapolação")}else{print("Interpolação")}## [1] "Interpolação"dados_novos <- data.frame(x1 = xnovo[2], x2 = xnovo[3], x3 = xnovo[4], x4 = xnovo[5])
predict(ajuste, newdata = dados_novos, interval = "prediction", level = 1-alfa)## fit lwr upr
## 1 9.974385 8.310434 11.63834- Preveja o valor da variável resposta quando \[\boldsymbol{x_i}^T=(1 \quad -3,4 \quad -1,7 \quad \quad 5,7 \quad \quad 8,9 ).\] Este caso é um caso de extrapolação ou interpolação? Justifique.
H = x%*%solve(t(x)%*%x)%*%t(x)
hmax = max(H)
xnovo = c(1 , -3.4 , -1.7 , 5.7 , 8.9 )
hnovo = xnovo%*%solve(t(x)%*%x)%*%xnovo
if(hnovo > hmax){print("Extrapolação")}else{print("Interpolação")}## [1] "Extrapolação"dados_novos <- data.frame(x1 = xnovo[2], x2 = xnovo[3], x3 = xnovo[4], x4 = xnovo[5])
predict(ajuste, newdata = dados_novos, interval = "prediction", level = 1-alfa)## fit lwr upr
## 1 22.87564 21.1329 24.618397 ANOVA
- Calcule o valor de SST, SSR e SSE usando as fórmulas abaixo:
\[\begin{align} SST &= \sum_{i=1}^n{(Y_i-\bar{Y})^2}\\ SSE &= \sum_{i=1}^n{(Y_i-\hat{Y_i})^2}\\ SSR &= \sum_{i=1}^n{(\hat{Y_i}-\bar{Y})^2} \end{align}\]
- Faça o teste de significância da regressão: \[H_0: \beta_1 = \ldots = \beta_p = 0 \text{ versus } H_1: \text{ pelo menos algum }\beta_j \neq 0 \]
Para isso calcule o valor da estatística de teste \(F_0 = \frac{SSR/p}{SSE/n-p-1}\) e o p-valor assumindo que \(F_0 \sim F_{(p, n-p-1)}\), quando \(H_0\) é verdadeiro. Compare os valores obtidos com os encontrados na saída do summary(ajuste).
SST = sum((y-mean(y))^2)
SSE = sum((y-ajuste$fitted.values)^2)
SSR = sum((ajuste$fitted.values-mean(y))^2)
c(SST, SSE, SSR)## [1] 22033.937 1014.326 21019.610MST = SST/(n-1)
MSE = SSE/(n-p-1)
MSR = SSR/p
c(MST, MSE, MSR)## [1] 22.055993 1.019423 5254.902625#calculando a estatistica de teste
F0 = MSR/MSE
F0## [1] 5154.779#calculando o pvalor
p_valor = 1-pf(F0, p, n-p-1)
p_valor## [1] 0if(p_valor <= alfa){print("Há indícios de que alguma covariavel seja significativa!")}else{print("Não rejeito a hipótese de que todas as covariáveis não sejam significativas!")}## [1] "Há indícios de que alguma covariavel seja significativa!"#comparando os valores obtidos
resultado##
## Call:
## lm(formula = y ~ x1 + x2 + x3 + x4, data = dados)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.1870 -0.6869 -0.0179 0.6622 2.9227
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.97439 0.04503 221.510 <2e-16 ***
## x1 2.00370 0.03301 60.697 <2e-16 ***
## x2 -1.52470 0.03231 -47.194 <2e-16 ***
## x3 3.03120 0.02516 120.461 <2e-16 ***
## x4 -0.01753 0.03167 -0.553 0.58
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.01 on 995 degrees of freedom
## Multiple R-squared: 0.954, Adjusted R-squared: 0.9538
## F-statistic: 5155 on 4 and 995 DF, p-value: < 2.2e-16c(resultado$fstatistic, F0)## value numdf dendf
## 5154.779 4.000 995.000 5154.779#o p-valor apresentado no summary nao da para ser extraido7.1 Função ANOVA do R
Considere as seguintes definições:
- \(SSR(x_1)\): representa a variação explicada da variável resposta quando apenas a covariável \(x_1\) é considerada no modelo;
- \(SSR(x_1, x_2)\): representa a variação explicada da variável resposta quando as covariáveis \(x_1\) e \(x_2\) são consideradas no modelo;
- \(SSR(x_1, x_2, x_3)\): representa a variação explicada da variável resposta quando as covariáveis \(x_1\), \(x_2\) e \(x_3\) são consideradas no modelo;
- \(SSR(x_1, x_2, x_3, x_4)\): representa a variação explicada da variável resposta quando todas as covariáveis \(x_1\), \(x_2\), \(x_3\) e \(x_4\) são consideradas no modelo;
- \(SSR(x_2 \mid x_1) = SSR(x_1, x_2) - SSR(x_1)\): representa a variação explicada exclusivamente pela covariável \(x_2\), após considerar o efeito da covariável \(x_1\) no modelo;
- \(SSR(x_3 \mid x_1, x_2) = SSR(x_1, x_2, x_3) - SSR(x_1, x_2)\): representa a variação explicada exclusivamente pela covariável \(x_3\), após considerar o efeito das covariáveis \(x_1, x_2\) no modelo;
- \(SSR(x_4 \mid x_1, x_2, x_3) = SSR(x_1, x_2, x_3, x_4) - SSR(x_1, x_2, x_3)\): representa a variação explicada exclusivamente pela covariável \(x_4\), após considerar o efeito das covariáveis \(x_1, x_2, x_3\) no modelo;
Calcule os valores de \(SSR(x_1), SSR(x_1,x_2), SSR(x_2|x_1), SSR(x_1, x_2, x_3), SSR(x_3|x_1, x_2), SSR(x_1, x_2, x_3, x_4), SSR(x_4|x_1, x_2, x_3)\).
Use a função ANOVA do R e compare com os valores encontrados acima.
Consigo obter os valores de SST, SSR e SSE através desta função?
#################################################
# Gerando a tabela ANOVA usando a função do R
#################################################
anova_tabela <- anova(ajuste)
anova_tabela## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## x1 1 4546.6 4546.6 4459.9829 <2e-16 ***
## x2 1 1676.9 1676.9 1644.9141 <2e-16 ***
## x3 1 14795.8 14795.8 14513.9118 <2e-16 ***
## x4 1 0.3 0.3 0.3063 0.5801
## Residuals 995 1014.3 1.0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ajuste_M1 = lm(y ~ x1, data = dados)
SSR_x1 = sum((ajuste_M1$fitted.values-mean(y))^2)
ajuste_M2 = lm(y ~ x1+x2, data = dados)
SSR_x1x2 = sum((ajuste_M2$fitted.values-mean(y))^2)
SSR_x2_dado_x1 = SSR_x1x2 - SSR_x1
ajuste_M3 = lm(y ~ x1+x2+x3, data = dados)
SSR_x1x2x3 = sum((ajuste_M3$fitted.values-mean(y))^2)
SSR_x3_dado_x1x2 = SSR_x1x2x3 - SSR_x1x2
SSR_x4_dado_x1x2x3 = SSR - SSR_x1x2x3
c(SSR_x1, anova_tabela$`Sum Sq`[1])## [1] 4546.611 4546.611c(SSR_x2_dado_x1, anova_tabela$`Sum Sq`[2])## [1] 1676.864 1676.864c(SSR_x3_dado_x1x2, anova_tabela$`Sum Sq`[3])## [1] 14795.82 14795.82c(SSR_x4_dado_x1x2x3, anova_tabela$`Sum Sq`[4])## [1] 0.3122792 0.3122792c(SSE, anova_tabela$`Sum Sq`[5], SSE/(n-p-1), anova_tabela$`Mean Sq`[5])## [1] 1014.326381 1014.326381 1.019423 1.019423anova_tabela## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## x1 1 4546.6 4546.6 4459.9829 <2e-16 ***
## x2 1 1676.9 1676.9 1644.9141 <2e-16 ***
## x3 1 14795.8 14795.8 14513.9118 <2e-16 ***
## x4 1 0.3 0.3 0.3063 0.5801
## Residuals 995 1014.3 1.0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# calculando SSR, SSE e SST pela função do R
# Soma dos Quadrados de Regressão (SSR)
SSR_2 <- sum(anova_tabela$`Sum Sq`[-length(anova_tabela$`Sum Sq`)])
# Soma dos Quadrados dos Erros (SSE)
SSE_2 <- anova_tabela$`Sum Sq`[length(anova_tabela$`Sum Sq`)]
# Soma dos Quadrados Total (SST)
SST_2 <- SSR_2 + SSE_2
cbind(c(SSR, SSE, SST), c(SSR_2, SSE_2, SST_2))## [,1] [,2]
## [1,] 21019.610 21019.610
## [2,] 1014.326 1014.326
## [3,] 22033.937 22033.9377.1.1 Valores F e p-valores
Calcule manualmente os valores encontrados nas 2 últimas colunas da tabela obtida pela função ANOVA do R.
# calculando F value sequencial
cbind( c(SSR_x1/MSE, SSR_x2_dado_x1/MSE,
SSR_x3_dado_x1x2/MSE, SSR_x4_dado_x1x2x3/MSE),
anova_tabela$`F value`[1:p])## [,1] [,2]
## [1,] 4.459983e+03 4.459983e+03
## [2,] 1.644914e+03 1.644914e+03
## [3,] 1.451391e+04 1.451391e+04
## [4,] 3.063292e-01 3.063292e-01#calculando p-valor
round( cbind( c(1-pf(SSR_x1/MSE, 1, n-p-1), 1-pf(SSR_x2_dado_x1/MSE, 1, n-p-1),
1-pf(SSR_x3_dado_x1x2/MSE, 1, n-p-1),
1-pf(SSR_x4_dado_x1x2x3/MSE, 1, n-p-1)),
anova_tabela$`Pr(>F)`[1:4]), 6)## [,1] [,2]
## [1,] 0.000000 0.000000
## [2,] 0.000000 0.000000
## [3,] 0.000000 0.000000
## [4,] 0.580066 0.5800668 O que obtenho com a função AOV do R?
#Funcao aov do R
aov(ajuste)## Call:
## aov(formula = ajuste)
##
## Terms:
## x1 x2 x3 x4 Residuals
## Sum of Squares 4546.611 1676.864 14795.823 0.312 1014.326
## Deg. of Freedom 1 1 1 1 995
##
## Residual standard error: 1.009665
## Estimated effects may be unbalancedrbind(round(c(SSR_x1, SSR_x2_dado_x1, SSR_x3_dado_x1x2, SSR_x4_dado_x1x2x3, SSE),3),
c(1,1,1,1, n-p-1))## [,1] [,2] [,3] [,4] [,5]
## [1,] 4546.611 1676.864 14795.82 0.312 1014.326
## [2,] 1.000 1.000 1.00 1.000 995.000sd(ajuste$residuals)## [1] 1.0076429 Coeficiente de determinação
Lembre que o coeficiente de determinação é dado por:
\[R^2 = 1 - \frac{SSE}{SST}\]
e o coeficiente de determinação ajustado é:
\[R^2_{ajustado} = 1 - \left( \frac{SSE / (n - p - 1)}{SST / (n - 1)} \right)\] Calcule estas medidas usando as fórmulas acima e compare com os valores encontrados no summary(ajuste).
R2 = 1 - SSE/SST
R2a = 1 - (SSE/(n-p-1))/(SST/(n-1))
c(R2, resultado$r.squared)## [1] 0.9539653 0.9539653c(R2a, resultado$adj.r.squared)## [1] 0.9537802 0.953780210 Exercício
Mude os valores de betas de forma que apenas a 3ª covariável seja significativa e analise a tabela da ANOVA.