Modelo de Regressão Linear Simples

Nesta aula computacional, analisaremos a função de verossimilhança de um MRLS e faremos análise de resíduos.

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.\]

1ª tarefa - Simulando conjunto de dados

Simule um conjunto de dados do modelo de regressão linear simples considerando n=1000 (tamanho amostral), \(\beta_0=10\), \(\beta_1=2\), \(\sigma^2=0,2\) e que \(x_i \stackrel{iid}{\sim} N(0,1)\). Em seguida, faça o gráfico de dispersão de X e Y.

# Definindo parâmetros reais
n       = 1000    # 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)) 
y       = beta0 + beta1 * x + e # Gerando variável dependente Y com erro aleatório

plot(x,y)

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

2ª tarefa - verossimilhança perfilada para beta0

Desenhe a função log-verossimilhança para \(\beta_0\) considerando \(\beta_1\) e \(\sigma^2\) conhecidos. E plote neste mesmo gráfico uma linha com o valor verdadeiro de \(\beta_0\).

# Função log-verossimilhança para beta_0
log_verossimilhanca = function(beta0, y, x, beta1, sigma2) {
  n     = length(y)
  soma  = -n/2 * log(2 * pi * sigma2) - 1/(2 * sigma2) * sum((y - beta0 - beta1 * x)^2)
  #soma2  = sum(dnorm(y, beta0+beta1*x, sqrt(sigma2), log=T))
  return(soma)
}

#comparando com uma funcao do R
log_verossimilhanca(beta0, y, x, beta1, sigma2)
## [1] -603.4228
sum(dnorm(y, beta0+beta1*x, sqrt(sigma2), log=T))
## [1] -603.4228
# Gerar um vetor de valores de beta_0
beta0_vals = seq(0, 20, by = 0.01)

# Calcular a log-verossimilhança para cada valor de beta0_vals
log_verossimilhanca_vals = sapply(beta0_vals, function(b0) log_verossimilhanca(b0, y, x, beta1, sigma2))

# Plotar a função de verossimilhança
plot(beta0_vals, log_verossimilhanca_vals, type = "l", col = "blue", xlab = expression(beta[0]), ylab = "Log-Verossimilhança",    main = expression("Função de log-verossimilhança para"~beta[0]), bty="n", lwd=2)
abline(v=beta0, lwd=2, col="red")

3ª tarefa - verossimilhança perfilada para beta1

Desenhe a função log-verossimilhança para \(\beta_1\) considerando \(\beta_0\) e \(\sigma^2\) conhecidos. E plote neste mesmo gráfico uma linha com o valor verdadeiro de \(\beta_1\).

# Gerar um vetor de valores de beta_1
beta1_vals = seq(0, 4, by = 0.01)

# Calcular a log-verossimilhança para cada valor de beta1_vals
log_verossimilhanca_vals = sapply(beta1_vals, function(b1) log_verossimilhanca(beta0, y, x, b1, sigma2))

# Plotar a função de verossimilhança
plot(beta1_vals, log_verossimilhanca_vals, type = "l", col = "blue", xlab = expression(beta[1]), ylab = "Log-Verossimilhança",    main = expression("Função de log-verossimilhança para"~beta[1]), bty="n", lwd=2)
abline(v=beta1, lwd=2, col="red")

4ª tarefa - verossimilhança perfilada para sigma2

Desenhe a função log-verossimilhança para \(\sigma^2\) considerando \(\beta_0\) e \(\beta_1\) conhecidos. E plote neste mesmo gráfico uma linha com o valor verdadeiro de \(\sigma^2\).

# Gerar um vetor de valores de beta_1
sigma2_vals = seq(0.1, 0.4, l=500)

# Calcular a log-verossimilhança para cada valor de sigma2_vals
log_verossimilhanca_vals = sapply(sigma2_vals, function(s2) log_verossimilhanca(beta0, y, x, beta1, s2))

# Plotar a função de verossimilhança
plot(sigma2_vals, log_verossimilhanca_vals, type = "l", col = "blue", xlab = expression(sigma^2), ylab = "Log-Verossimilhança",    main = expression("Função de log-verossimilhança para"~sigma^2), bty="n", lwd=2)
abline(v=sigma2, lwd=2, col="red")

5ª tarefa - verossimilhança para beta0 e beta1

Desenhe a função log-verossimilhança para \(\beta_0\) e \(\beta_1\), considerando \(\sigma^2\) conhecido. Plote neste gráfico o valor verdadeiro de \(\beta_0\) e \(\beta_1\).

# Gerar um vetor de valores de beta_1
beta0_vals = seq(0, 20, l=500)
beta1_vals = seq(0, 4, l=500)

# Calcular a log-verossimilhança para cada valor 
log_verossimilhanca_vals = sapply(sigma2_vals, function(s2) log_verossimilhanca(beta0, y, x, beta1, s2))


# Criar uma matriz para armazenar os valores da log-verossimilhança
log_verossimilhanca_matrix = outer(beta0_vals, beta1_vals, 
                           Vectorize(function(b0, b1) log_verossimilhanca(b0, y, x, b1, sigma2)))

# Plotar a lo-verossimilhança conjunta em 3D
library(plotly)
plot_ly(z = ~log_verossimilhanca_matrix, x = ~beta0_vals, y = ~beta1_vals, type = "surface") %>%
  add_markers(x = beta0, y = beta1, z = log_verossimilhanca(beta0, y, x, beta1, sigma2), 
              marker = list(color = 'red', size = 5, symbol = 'circle')) %>%
  layout(scene = list(xaxis = list(title = "beta0"),
                      yaxis = list(title = "beta1"),
                      zaxis = list(title = "Log-Verossimilhança")))

6ª tarefa - curva de contorno

Faça curvas de contorno para \(\beta_0\) e \(\beta_1\), considerando \(\sigma^2\) conhecido e interprete estas curvas.

# Gerar um vetor de valores de beta_1
beta0_vals = seq(0, 20, l=500)
beta1_vals = seq(0, 4, l=500)

# Fazer o gráfico de contorno
contour(beta0_vals, beta1_vals, log_verossimilhanca_matrix, 
        xlab = expression(beta[0]), 
        ylab = expression(beta[1]), 
        main = expression("Curvas de Contorno da Verossimilhança para"~beta[0]~"e"~beta[1]))

# Adicionar o ponto verdadeiro de beta_0 e beta_1
points(beta0, beta1, col = "red", pch = 19, cex = 1.5)

Interpretação das Curvas de Contorno

O gráfico de contorno mostra as curvas de nível da função de log-verossimilhança em função de \(\beta_0\) e \(\beta_1\). As curvas representam diferentes níveis de verossimilhança, indicando combinações de valores de \(\beta_0\) e \(\beta_1\) que resultam em probabilidades similares de observação dos dados sob o modelo de regressão linear. No gráfico, o ponto vermelho representa os valores verdadeiros dos parâmetros \(\beta_0 = 10\) e \(\beta_1 = 2\). Esse ponto está próximo ao centro das curvas de contorno de nível mais elevado, o que sugere que esses valores são os que maximizam a verossimilhança. Ou seja, os valores \(\beta_0 = 10\) e \(\beta_1 = 2\) são os mais prováveis dados os resultados calculados.

  • As curvas de contorno mais próximas ao ponto verdadeiro indicam regiões de maior verossimilhança, onde os valores de \(\beta_0\) e \(\beta_1\) fornecem um bom ajuste ao modelo.
  • À medida que nos afastamos do ponto central, a log-verossimilhança diminui, o que indica que os valores de \(\beta_0\) e \(\beta_1\) nessas regiões são menos prováveis de explicar os dados observados.
  • As curvas de contorno têm uma forma elíptica, o que reflete a incerteza associada às estimativas. Se as curvas são mais alongadas em uma direção, isso significa que há maior incerteza em relação a esse parâmetro.
  • A densidade das curvas também nos dá uma ideia de quão sensível a verossimilhança é às mudanças em \(\beta_0\) e \(\beta_1\). Se as curvas são muito próximas, pequenas variações nesses parâmetros resultam em grandes mudanças na verossimilhança.

7ª tarefa: ANOVA

  • Ajuste um modelo de regressão linear simples.
  • Calcula as medidas abaixo: \[SST = \sum_{i=1}^n{(Y_i - \bar{Y})^2}, \quad SSR = \sum_{i=1}^n{(\hat{Y}_i - \bar{Y})^2}, \quad SSE = \sum_{i=1}^n{(Y_i - \hat{Y}_i)^2},\] \[MSR = \frac{SSR}{1}, \quad MSE = \frac{SSE}{n-2}, \quad F_0 = \frac{MSR}{MSE} \mbox{ e } \text{p-valor} = P(W>F_0), \quad W \sim F-Snedecor(1, n-2).\]
  • Em seguida, utilize a função ANOVA do R para conferir os valores encontrados.
#Fazendo o ajuste do MRLS
ajuste = lm(y~x)
#Criando a tabela da ANOVA
TabelaANOVA = anova(ajuste)
#conferindo os valores encontrados
SSR   = ajuste$coefficients[2]*sum((x-mean(x))*y) #SSR
SSR2  = sum((ajuste$fitted.values-mean(y))^2) #SSR
SSRes = sum((y-ajuste$fitted.values)^2) #SSRes
SST   = sum((y-mean(y))^2) #SST
n     = length(y)
F0 = SSR/(SSRes/(n-2))
pf(F0, 1, n-2, lower.tail = F)
## x 
## 0
TabelaANOVA
## Analysis of Variance Table
## 
## Response: y
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## x           1 4139.4  4139.4   21191 < 2.2e-16 ***
## Residuals 998  194.9     0.2                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
c(SSR, SSR2, SSRes, SST, F0)
##          x                                           x 
##  4139.4194  4139.4194   194.9436  4334.3630 21191.4639
SSRes+SSR
##        x 
## 4334.363

8ª tarefa: Análise de Resíduos

residuos = y - ajuste$fitted.values #residuals(ajuste)
plot(round(residuos-ajuste$residuals,9))

sigma2chapeu = summary(ajuste)$sigma^2
hist(residuos, freq=F, main="Histograma", ylab="densidade")
curve(dnorm(x, 0, sqrt(sigma2chapeu)), add=T, lwd=2, col="red")

qqnorm(residuos, xlab="quatis teóricos", 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.99861, p-value = 0.6311
ks.test(residuos,"pnorm",0, sqrt(sigma2chapeu))
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  residuos
## D = 0.022381, p-value = 0.6983
## alternative hypothesis: two-sided
ks.test(residuos,"pnorm",mean(residuos), sd(residuos))
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  residuos
## D = 0.022452, p-value = 0.6946
## alternative hypothesis: two-sided
#Um p-valor pequeno indica que os resíduos não seguem uma distribuição normal.

par(mfrow=c(2,2))
plot(x, residuos)
plot(y, residuos)
plot(residuos)

# Testando Homocedasticidade (Teste de Breusch-Pagan)
# Instalar pacote necessário (se ainda não instalado)
# install.packages("lmtest")
library(lmtest)
bptest(ajuste)
## 
##  studentized Breusch-Pagan test
## 
## data:  ajuste
## BP = 0.029704, df = 1, p-value = 0.8632
#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(ajuste)
## 
##  Durbin-Watson test
## 
## data:  ajuste
## DW = 2.0752, p-value = 0.8833
## 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.

9ª tarefa: Dados sem linearidade

Altere a variável resposta \(Y_i\) usando \(X_i^2\) ao invés de \(X_i\) na geração dos dados, para \(i=1, \ldots, n\). Depois ajuste um MRLS e analise os resíduos.

#Simulando os dados
ySL       = beta0 + beta1 * x^2 + e # Gerando variável dependente Y com erro aleatório
plot(x, ySL)

plot(x^2, ySL)

#Fazendo o ajuste
ajuste = lm(ySL~x)

#Analisando os residuos
residuos = ajuste$residuals
sigma2chapeu = summary(ajuste)$sigma^2
ajuste$coefficients
## (Intercept)           x 
##  12.0254599   0.2321536
sigma2chapeu
## [1] 7.852317
hist(residuos, freq=F, main="Histograma", ylab="densidade")
curve(dnorm(x, 0, sqrt(sigma2chapeu)), add=T, lwd=2, col="red")

qqnorm(residuos, xlab="quatis teóricos", 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.77333, p-value < 2.2e-16
ks.test(residuos,"pnorm",0, sqrt(sigma2chapeu))
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  residuos
## D = 0.1714, p-value < 2.2e-16
## alternative hypothesis: two-sided
ks.test(residuos,"pnorm",mean(residuos), sd(residuos))
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  residuos
## D = 0.17143, p-value < 2.2e-16
## alternative hypothesis: two-sided
#Um p-valor pequeno indica que os resíduos não seguem uma distribuição normal.

par(mfrow=c(2,2))
plot(x, residuos)
plot(ySL, residuos)
plot(residuos)

# Testando Homocedasticidade (Teste de Breusch-Pagan)
# Instalar pacote necessário (se ainda não instalado)
# install.packages("lmtest")
library(lmtest)
bptest(ajuste)
## 
##  studentized Breusch-Pagan test
## 
## data:  ajuste
## BP = 1.846, df = 1, p-value = 0.1743
#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(ajuste)
## 
##  Durbin-Watson test
## 
## data:  ajuste
## DW = 2.1179, p-value = 0.969
## 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.

10ª tarefa: Dados sem normalidade

Altere a variável resposta de forma que \(Y_i = \beta_0+\beta_1 X_i + \epsilon_i\), com \(\epsilon_i \sim Unif(-2, 2)\). Depois ajuste um MRLS e analise os resíduos.

#Simulando os dados
epsilon       = runif(n, -2, 2) 
ySN       = beta0 + beta1 * x + epsilon # Gerando variável dependente Y com erro aleatório
plot(x, ySN)

#Fazendo o ajuste
ajuste = lm(ySN~x)

#Analisando os residuos
residuos = ajuste$residuals
sigma2chapeu = summary(ajuste)$sigma^2
ajuste$coefficients
## (Intercept)           x 
##   10.069156    2.000282
sigma2chapeu
## [1] 1.352928
hist(residuos, freq=F, main="Histograma", ylab="densidade")
curve(dnorm(x, 0, sqrt(sigma2chapeu)), add=T, lwd=2, col="red")

qqnorm(residuos, xlab="quatis teóricos", 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.95118, p-value < 2.2e-16
ks.test(residuos,"pnorm",0, sqrt(sigma2chapeu))
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  residuos
## D = 0.067102, p-value = 0.0002455
## alternative hypothesis: two-sided
ks.test(residuos,"pnorm",mean(residuos), sd(residuos))
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  residuos
## D = 0.067213, p-value = 0.0002383
## alternative hypothesis: two-sided
#Um p-valor pequeno indica que os resíduos não seguem uma distribuição normal.

par(mfrow=c(2,2))
plot(x, residuos)
plot(ySN, residuos)
plot(residuos)

# Testando Homocedasticidade (Teste de Breusch-Pagan)
# Instalar pacote necessário (se ainda não instalado)
# install.packages("lmtest")
library(lmtest)
bptest(ajuste)
## 
##  studentized Breusch-Pagan test
## 
## data:  ajuste
## BP = 6.7718, df = 1, p-value = 0.009261
#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(ajuste)
## 
##  Durbin-Watson test
## 
## data:  ajuste
## DW = 2.0837, p-value = 0.9074
## 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.

11ª tarefa: Dados sem homocedasticidade

Altere a variável resposta de forma que \(Y_i = \beta_0+\beta_1 X_i + \nu_i\), com \(\nu_i \sim N(0, \sigma^2 X_i^2)\), para \(i=1, \ldots, n\). Depois ajuste um MRLS e analise os resíduos.

#Simulando os dados
nu       = NULL
for(i in 1:n){nu[i] = rnorm(1, 0, sqrt(sigma2*x[i]*x[i]))}
ySH       = beta0 + beta1 * x + nu # Gerando variável dependente Y com erro aleatório
plot(x, ySH)

#Fazendo o ajuste
ajuste = lm(ySH~x)

#Analisando os residuos
residuos = ajuste$residuals
sigma2chapeu = summary(ajuste)$sigma^2
ajuste$coefficients
## (Intercept)           x 
##   10.013153    1.997708
sigma2chapeu
## [1] 0.1933804
hist(residuos, freq=F, main="Histograma", ylab="densidade")
curve(dnorm(x, 0, sqrt(sigma2chapeu)), add=T, lwd=2, col="red")

qqnorm(residuos, xlab="quatis teóricos", 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.93946, p-value < 2.2e-16
ks.test(residuos,"pnorm",0, sqrt(sigma2chapeu))
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  residuos
## D = 0.099707, p-value = 4.634e-09
## alternative hypothesis: two-sided
ks.test(residuos,"pnorm",mean(residuos), sd(residuos))
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  residuos
## D = 0.099667, p-value = 4.709e-09
## alternative hypothesis: two-sided
#Um p-valor pequeno indica que os resíduos não seguem uma distribuição normal.

par(mfrow=c(2,2))
plot(x, residuos)
plot(ySH, residuos)
plot(residuos)

# Testando Homocedasticidade (Teste de Breusch-Pagan)
# Instalar pacote necessário (se ainda não instalado)
# install.packages("lmtest")
library(lmtest)
bptest(ajuste)
## 
##  studentized Breusch-Pagan test
## 
## data:  ajuste
## BP = 0.88798, df = 1, p-value = 0.346
#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(ajuste)
## 
##  Durbin-Watson test
## 
## data:  ajuste
## DW = 1.9318, p-value = 0.1405
## 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.