Modelo de Regressão Linear Simples

Nesta 4ª aula computacional, faremos intervalos de confiança para as medidas desconhecidas.

Modelo

Considere o modelo de regressão linear simples: \[Y_i = \beta_0 + \beta_1 X_i + e_i, \quad e_i \stackrel{iid}{\sim} N(0,\sigma^2), \quad i=1,\ldots,n.\] Suponha \(\beta_0\), \(\beta_1\) e \(\sigma^2\) desconhecidos. Considere os seguintes estimadores não viesados para estes parâmetros: \[\hat{\beta}_0 = \bar{y}-\bar{x}\hat{\beta}_1, \quad \hat{\beta}_1 = \frac{\sum_{i=1}^n{y_ix_i - n\bar{x}\bar{y}}}{\sum_{i=1}^n{(x_i-\bar{x})^2}}, \quad \hat{\sigma}^2 = \frac{\sum_{i=1}^n{(y_i-\hat{y}_i)^2}}{n-2}, \mbox{ sendo } \hat{y}_i = \hat{\beta}_0 + x_i \hat{\beta}_1.\]

Distribuições importantes

\[\frac{\hat{\beta}_0 - \beta_0}{\sqrt{\hat{\sigma}^2\left(\frac{1}{n}+\frac{\bar{x}^2}{\sum_{i=1}^n{(x_i-\bar{x})^2}}\right)}} \sim t_{(n-2)} , \quad \frac{\hat{\beta}_1 - \beta_1}{\sqrt{\frac{\hat{\sigma}^2}{\sum_{i=1}^n{(x_i-\bar{x})^2}}}} \sim t_{(n-2)}, \quad \frac{(n-2)\hat{\sigma}^2}{\sigma^2} \sim \chi^2_{(n-2)}.\]

\[\frac{\hat{y}_i - \mu_i}{\sqrt{\hat{\sigma}^2\left(\frac{1}{n}+\frac{(x_i-\bar{x})^2}{\sum_{i=1}^n{(x_i-\bar{x})^2}}\right)}} \sim t_{(n-2)}, \quad \frac{\hat{y}_{novo} - y_{novo}}{\sqrt{\hat{\sigma}^2\left(1+\frac{1}{n}+\frac{(x_i-\bar{x})^2}{\sum_{i=1}^n{(x_i-\bar{x})^2}}\right)}} \sim t_{(n-2)}.\]

1ª tarefa

Simule um conjunto de dados do modelo de regressão linear simples considerando n=100 (tamanho amostral), \(\beta_0=10\), \(\beta_1=2\), \(\sigma^2=0,2\) e que \(x_i \stackrel{iid}{\sim} N(0,1)\).

# Definindo parâmetros reais
n       = 100    # Tamanho da amostra
beta0   = 10  # Intercepto
beta1   = 2  # Coeficiente de regressão
sigma2  = 0.2

#Simulando os dados
x       = rnorm(n, mean = 0, sd = 1) # Gerando variável independente X
e       = rnorm(n, mean = 0, sd = sqrt(sigma2)) 
mu      = beta0 + beta1 * x
y       = mu + e # Gerando variável dependente Y com erro aleatório

Considere agora que os parâmetros são desconhecidos.

2ª tarefa

Encontre intervalos de confiança bilaterais de 95% para \(\beta_0\), \(\beta_1\) e \(\sigma^2\).

# Ajustando o modelo de regressão linear simples
ajuste = lm(y ~ x)
# Armazenando as estimativas dos parametros
beta0chapeu  = ajuste$coefficients[1]  # Estimativa de beta0
beta1chapeu  = ajuste$coefficients[2]  # Estimativa de beta1
sigma2chapeu = summary(ajuste)$sigma^2 # Estimativa de sigma2

#Definindo algumas variaveis
gama    = 0.95
alfa    = 1 - gama
xbarra  = mean(x)
Sxx     = sum((x-xbarra)^2)

# Intervalo de confiança para beta0
c(beta0chapeu - qt(alfa/2,n-2, lower.tail = F)*sqrt(sigma2chapeu*(1/n + xbarra^2/Sxx)),
beta0chapeu + qt(alfa/2,n-2, lower.tail = F)*sqrt(sigma2chapeu*(1/n + xbarra^2/Sxx)))

## (Intercept) (Intercept) 
##    9.958899   10.110269

c(beta1chapeu - qt(alfa/2,n-2, lower.tail = F)*sqrt(sigma2chapeu/Sxx),
beta1chapeu + qt(alfa/2,n-2, lower.tail = F)*sqrt(sigma2chapeu/Sxx))

##        x        x 
## 1.921925 2.074996

#IC para beta0 e beta1 usando funcao do R
confint(ajuste, level = gama)

##                2.5 %    97.5 %
## (Intercept) 9.958899 10.110269
## x           1.921925  2.074996

# Intervalo de confiança para sigma^2
c((n-2) * sigma2chapeu / qchisq(1-alfa/2, n-2), 
  (n-2) * sigma2chapeu / qchisq(alfa/2, n-2))

## [1] 0.1114108 0.1955920

3ª tarefa

Encontre intervalos de confiança bilaterais de 95% para \(\mu_i = E[Y_i]\), \(i=1, \ldots, n\).

Em seguida, faça um gráfico de dispersão e inclua a reta ajustada e segmentos com estes intervalos.

# Intervalo de confiança para mu_i
quantil    = qt(alfa/2,n-2, lower.tail = F)
DP_ychapeu = sqrt(sigma2chapeu*(1/n + (x-xbarra)^2/Sxx))
muchapeu   = beta0chapeu + beta1chapeu*x
linf_mui   = muchapeu - quantil*DP_ychapeu
lsup_mui   = muchapeu + quantil*DP_ychapeu

#usando funcao do R
pred = predict(ajuste, interval = "confidence", level = gama)

#comparando muchapeu calculado manualmente com o obtido pela funcao lm
round(summary(muchapeu - ajuste$fitted.values),6)

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

#verificando se os limites inferiores obtidos manualmente sao iguais aos obtidos pela funcao predict
round(summary(pred[,2] - linf_mui),6)

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

#verificando se os limites superiores obtidos manualmente sao iguais aos obtidos pela funcao predict
round(summary(pred[,3] - lsup_mui),6)

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

# Grafico com os intervalos de confianca para mu
plot(1:n, mu, main = "", pch=1,
     xlab = "unidade amostral", ylab = expression(mu), bty="n", 
     ylim=c(min(mu, pred), max(mu, pred)+10))
#Adicionando segmentos dos intervalos de confiança
#for(i in 1:n){  segments(i, pred[i,2], i, pred[i,3], col = "blue", lwd=2)}
segments(1:n, pred[,2], 1:n, pred[,3], col = "blue", lwd=2)
points(1:n, ajuste$fitted.values,pch=4, col="blue")
points(1:n,y,pch=19, col="green")
legend("topleft",
       legend = c("Intervalos de confiança", "Valores esperados E(Y|X=x)","Valores observados","Y"),
       col = c("blue", "blue","black","green"),
       lty = c(1, 1, 1, 1),
       pch = c(NA, 4, 1, 19),
       lwd = c(2, NA, NA, NA),
       bty = "n")

# Grafico de dispersao com reta ajustada e intervalos de confianca
plot(x, mu, main = "Gráfico de Dispersão com Intervalos de Confiança para as médias", col="black", pch=1,
     xlab = "x", ylab = "", bty="n", ylim=c(min(y, pred), max(y, pred)))
abline(ajuste, col = "red", lwd = 2)
# Adicionando segmentos dos intervalos de confiança
#for(i in 1:n){  segments(x[i], pred[i,2], x[i], pred[i,3], col = "blue",lwd=2)}
segments(x, pred[,2], x, pred[,3], col = "blue",lwd=2)
points(x, y, pch=19, col="green")
points(x,ajuste$fitted.values,pch=4, col="blue")
legend("topleft",
       legend = c("Reta ajustada", "Intervalos de confiança", "Valores esperados mu=E(Y|X=x)","Valores observados para mu", "Y"),
       col = c("red", "blue", "blue","black","green"),
       lty = c(1, 1, 1, 1,1),
       pch = c(NA, NA, 4, 1, 19),
       lwd = c(2, 2, NA,NA,NA),
       bty = "n")

4ª tarefa

Repita o gráfico anterior, trocando os segmentos por uma região sombreada.

# Grafico de dispersao com reta ajustada e intervalos de confianca
plot(x, y, main = "Gráfico de Dispersão com Intervalos de Confiança para as médias", col="green",
     xlab = "x", ylab = "y", bty="n", ylim=c(min(y, pred), max(y, pred)))
# Ordenando os valores de x para garantir que o sombreado seja desenhado corretamente
ordem = order(x)
# Adicionando o sombreado dos intervalos de confiança
polygon(c(x[ordem], rev(x[ordem])), 
        c(pred[ordem,2], rev(pred[ordem,3])), 
        col = "blue", border = NA)
#adicionando a reta ajustada
abline(ajuste, col = "red", lwd = 2)
#adicionando os demais pontos
points(x,y,col="green",pch=19)
points(x,ajuste$fitted.values, col="blue",pch=4)
points(x, mu, col="black",pch=1)
legend("topleft",
       legend = c("Reta ajustada", "Intervalos de confiança", "Valores esperados mu=E(Y|X=x)", "mu", "Y"),
       col = c("red", "blue", "blue", "black", "green"),
       lty = c(1, 1, NA,NA, NA),
       lwd = c(2, 2, NA,NA, NA),
       pch = c(NA, NA, 4, 1, 19),
       #pt.cex = c(1, 1, 2, 2, 2),
       bty = "n")

5ª tarefa

Faça um gráfico de dispersão e inclua uma região sombreada com o intervalo de predição para \(Y_{novo}\) para todo \(x_{novo}\) maior do que \(x_{min}\) e menor do que \(x_{max}\).

# Intervalo de predicao para y_novo
quantil = qt(alfa/2,n-2, lower.tail = F)
DP_yi   = sqrt(sigma2chapeu*(1 + 1/n + (x-xbarra)^2/Sxx))
linf_y  = muchapeu - quantil*DP_yi
lsup_y  = muchapeu + quantil*DP_yi

#usando uma funcao do R
xnovo = x
pred_y = predict(ajuste, newdata=data.frame(x), interval = "prediction", level = gama)

#verificando se os limites inferiores obtidos manualmente sao iguais aos obtidos pela funcao predict
round(summary(pred_y[,2] - linf_y),6)

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

#verificando se os limites superiores obtidos manualmente sao iguais aos obtidos pela funcao predict
round(summary(pred_y[,3] - lsup_y),6)

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

# Gráfico de dispersão com reta ajustada e intervalos de confiança
plot(x, y, main = "Gráfico de Dispersão com ICs e IPs",
     xlab = "x", ylab = "y", bty="n", ylim=c(min(y, pred_y), max(y, pred_y)), col="green", pch=19)
# Ordenando os valores de x para garantir que o sombreado seja desenhado corretamente
ordem = order(x)
# Adicionando o sombreado dos intervalos de confianca de y
polygon(c(x[ordem], rev(x[ordem])), 
        c(pred_y[ordem,2], rev(pred_y[ordem,3])), 
        col = "darkgray", border = NA)
# Adicionando o sombreado dos intervalos de confianca de mu
polygon(c(x[ordem], rev(x[ordem])), 
        c(pred[ordem,2], rev(pred[ordem,3])), 
        col = "blue", border = NA)
#adicionando a reta ajustada
abline(ajuste, col = "red", lwd = 2)
#adicionando os valores observados
points(x,y, col="green", pch=19)
#legenda
legend("topleft",
       legend = c("Reta ajustada", "IC para mu", "IP para Y", "Y"),
       col = c("red", "blue", "darkgray", "green"),
       lty = c(1, 1, 1, NA),
       lwd = c(2, 2, 2, NA),
       pch = c(NA, NA, NA, 19),
       bty = "n")

6ª tarefa

Repita da 2ª a 5ª tarefa para os dados abaixo.

x = c(-.62, -1.61, 1.3, -0.95, -0.93, 0.19, -0.24, -0.79, 1.05, -0.95)
y = c(6.91, 2.72, 13.37, 5.62, 6.62, 11.52, 11.39, 8.29, 10.14, 6.9)

DESAFIO

Calcular para os 2 conjuntos de dados usados acima:

Quantas observações pertencem aos intervalos de confiança para as médias?
Quantas observações pertencem aos intervalos previstos para Y?

Aula de Modelos Lineares I

Patrícia Lusié Velozo da Costa

2025-04-25