Unidade III - Modelo de Regressão Linear

#Para numerar tabelas e figuras:
require(captioner)
fig_nums <- captioner(prefix = "Figura")
table_nums <- captioner(prefix = "Tabela")
quadro_nums<- captioner(prefix="Quadro")

require(DT)
require(tidyverse)

1 Unidade III - Inferência Bayesiana

1.1 Exemplo 3.1: Regressão linear simples

O problema envolve as variáveis \(X\): a dose de um medicamento anti-alérgico em estudo, e \(Y\): tempo de duração do efeito (alívio dos sintomas alérgicos). Abaixo temos a representação gráfica dos dados observados. Pelo gráfico, nós concluímos que uma relação linear (reta) é satisfatória para os dados. Também iremos supor que \(X\) é uma variável controlada pelo pesquisador (sem a presença de erros).

x=c(3,3,4,5,6,6,7,8,8,9)
y=c(9,5,12,9,14,16,22,18,24,22)
a=cbind(x,y)%>%
  as.data.frame()
names(a)=c("x_dose","y_tempo")

tab_nums <- captioner(prefix = "Tabela")
tab_cap = tab_nums("tabela1","Dados do Exemplo")
datatable(as.data.frame(a),caption=tab_cap)

ggplot(a,aes(x=x_dose,y=y_tempo))+
       geom_point() +
      labs(x="Dose",y="Tempo")

Figura 1: Diagrama de dispersão dos dados

Modelo Estatístico O modelo estatístico será o modelo de regressão simples com erros i.i.d. normais: \[ y_i=\beta_0+\beta_1 x_i + \epsilon_i, i=1,\ldots,n, \] onde \(\beta_0\): intercepto da linha de regressão com o eixo \(y\);
\(\beta_1\): coeficiente de inclinação da reta: é o nº de unidades em \(y\) que mudam para cada unidade da variável independente \(x\).
\(\epsilon_i\): erros aleatórios com distribuição normal: \(\epsilon_i \sim N(0,\sigma^2)\).

Estimadores de mínimos quadrados da regressão Encontrar \(\widehat{\beta_0}\) e \(\widehat{\beta_1}\) que minimizam a soma de quadrados dos erros:
\(S(\beta_0,\beta_1)=\sum \limits_{i=1}^{n} \epsilon_i^2 = \sum \limits_{i=1}^{n} (y_i-\beta_0-\beta_1 x_i)^2\).
Então temos as equações normais:

\[ \frac{\partial S} {\partial \beta_0} =0 \text{ e } \frac{\partial S} {\partial \beta_1} =0. \]

A solução é dada por:

\[ \begin{array}{lll} \widehat{\beta_1}&=&\frac{\sum \limits_{i=1}^{n} (x_i-\bar{x})(y_i-\bar{y})}{\sum \limits_{i=1}^{n} (x_i-\bar{x})^2}\\ \widehat{\beta_0}&=&\bar{y}-\widehat{\beta_1}\bar{x}. \end{array} \]

Com respeito a variância \(\sigma_2\): \[ \widehat{\sigma_2}=\frac{SQR}{n-2}, \] ou seja, a estimativa da variância é igual à soma dos quadrados dos resíduos sobre o número graus de liberdade do modelo. Os intervalos de confiança e testes de hipóteses para \(\beta_0\) e \(\beta_1\) são baseados na distribuição t-student.

Voltando ao R:
- Método dos mínimos quadrados: cálculo das estimativas passo a passo

soma_xx=sum((x-mean(x))^2)
soma_xy=sum((x-mean(x))*(y-mean(y)))
beta1_hat=soma_xy/soma_xx
beta0_hat=mean(y)-beta1_hat*mean(x)
print(paste0("beta0_hat=",beta0_hat," & beta1_hat=",beta1_hat))

## [1] "beta0_hat=-1.07090464547677 & beta1_hat=2.74083129584352"

• Método dos mínimos quadrados: utilizando a função lm:linear model
  • Resíduos aproximadamente normais com variância constante.

a=lm(y ~ x)
summary(a)

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.6333 -2.0128 -0.3741  2.0428  3.8851 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -1.0709     2.7509  -0.389 0.707219    
## x             2.7408     0.4411   6.214 0.000255 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.821 on 8 degrees of freedom
## Multiple R-squared:  0.8284, Adjusted R-squared:  0.8069 
## F-statistic: 38.62 on 1 and 8 DF,  p-value: 0.0002555

anova(a)

## Analysis of Variance Table
## 
## Response: y
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## x          1 307.247 307.247  38.615 0.0002555 ***
## Residuals  8  63.653   7.957                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(a)

shapiro.test(residuals(a))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(a)
## W = 0.93227, p-value = 0.4706

2 Inferência Bayesiana no modelo de regressão linear simples

Assumimos as seguintes distribuições a priori:
\(\beta_0 \sim N(0,a_0^2)\) com \(a_0\) conhecido \(\beta_0 \sim N(0,a_1^2)\) com \(a_1\) conhecido \(\sigma_2 \sim IG(b,d)\) com \(b\) e \(d\) conhecidos Iremos assumir independência a priori entre os parâmetros.

• Função de verossimilhança:

\[ \begin{array}{lll} L(\beta_0,\beta_1,\sigma^2)&=&\prod \limits_{i=1}^n \frac{1}{2 \pi \sigma^2} \exp\left[-\frac{1}{2\sigma^2}(y_i-\beta_0-\beta_1x_i)^2\right]\\ &=& (\sigma^2)^{-\frac{n}{2}} \exp\left[-\frac{1}{2\sigma^2}\sum \limits_{i=1}^n (y_i-\beta_0-\beta_1x_i)^2\right] \end{array} \]

• Distribuição a posteriori conjunta para \(\beta_0\),\(\beta_1\) e \(\sigma^2\) é dada por:

\[ \begin{array}{lll} f(\beta_0,\beta_1,\sigma^2 \mid \boldsymbol{x}, \boldsymbol{y}) &\propto& \exp\left[-\frac{\beta_0^2}{2a_0^2}\right] \exp\left[-\frac{\beta_1^2}{2a_1^2}\right] (\sigma^2)^{-(b+1)} \exp\left[-\frac{d}{\sigma^2}\right] (\sigma^2)^{-\frac{n}{2}} \exp\left[-\frac{1}{2\sigma^2}\sum \limits_{i=1}^n (y_i-\beta_0-\beta_1x_i)^2\right]\\ &\propto& \exp\left[-\frac{\beta_0^2}{2a_0^2}\right] \exp\left[-\frac{\beta_1^2}{2a_1^2}\right] (\sigma^2)^{-(b+\frac{n}{2}+1)} \exp\left[-\frac{d}{\sigma^2}\right] \exp\left[-\frac{1}{2\sigma^2}\sum \limits_{i=1}^n (y_i-\beta_0-\beta_1x_i)^2\right] \end{array} \]

• Distribuições a posteriori marginais - é necessário integrar com respeito aos outros parâmetros, respeitando os limites de integração:
  \(\checkmark\) A conjunta de \(\beta_0\) e \(\beta_1\) é obtida da integração com respeito à variância \(\sigma^2\): \[ f(\beta_0,\beta_1 \mid \boldsymbol{x}, \boldsymbol{y}) = \int \limits _0^\infty f(\beta_0,\beta_1,\sigma^2 \mid \boldsymbol{x}, \boldsymbol{y}) d\sigma^2\\ \]

\(\checkmark\) A marginal de \(\sigma^2\) é obtida da integração com respeito à \(\beta_0\) e \(\beta_1\): \[ f(\sigma^2\mid \boldsymbol{x}, \boldsymbol{y}) = \int \limits _{-\infty}^\infty \int \limits _{-\infty}^\infty f(\beta_0,\beta_1,\sigma^2 \mid \boldsymbol{x}, \boldsymbol{y}) d\beta_0 d\beta_1 \]

\(\checkmark\) A marginal de \(\beta_0\) é obtida da integração com respeito à \(\beta_1\): \[ f(\beta_0\mid \boldsymbol{x}, \boldsymbol{y}) = \int \limits _{-\infty}^\infty f(\beta_0,\beta_1 \mid \boldsymbol{x}, \boldsymbol{y}) d\beta_1 \]

\(\checkmark\) A marginal de \(\beta_1\) é obtida da integração com respeito à \(\beta_0\): \[ f(\beta_1\mid \boldsymbol{x}, \boldsymbol{y}) = \int \limits _{-\infty}^\infty f(\beta_0,\beta_1 \mid \boldsymbol{x}, \boldsymbol{y}) d\beta_0 \]

\(\checkmark\) A conjunta de \(\beta_0\) e \(\beta_1\) pode ser obtida analiticamente:

\[ \begin{array}{lll} f(\beta_0,\beta_1 \mid \boldsymbol{x}, \boldsymbol{y}) &\propto& \exp\left[-\frac{\beta_0^2}{2a_0^2}\right] \exp\left[-\frac{\beta_1^2}{2a_1^2}\right] \int \limits _{0}^\infty (\sigma^2)^{-(b+\frac{n}{2}+1)} \exp\left[-\frac{1}{\sigma^2} \left(d + \frac{1}{2}\sum \limits_{i=1}^n (y_i-\beta_0-\beta_1x_i)^2\right)\right]d\sigma^2\\ &\propto& \exp\left[-\frac{\beta_0^2}{2a_0^2}\right] \exp\left[-\frac{\beta_1^2}{2a_1^2}\right] \left[d + \frac{1}{2}\sum \limits_{i=1}^n (y_i-\beta_0-\beta_1x_i)^2\right]^{-\left(b+\frac{n}{2}\right)}. \end{array} \]

Este resultado envolve a integral da distribuição
Gamma Invertida \(\left(b+\frac{n}{2},d+\frac{1}{2}\sum \limits_{i=1}^n (y_i - \beta_0 -\beta_1 x_i)^2\right)\) e o fato de que a integral de uma função de densidade sempre é igual a 1.
Observe que mesmo tendo obtido a integral, ela não tem forma conhecida - não conseguimos identificar esta na tábua de distribuições com suporte de \(-\infty\) a \(\infty\).

Já as outras marginais não têm forma fechada - não são obtidas analiticamente - integrais analíticas não são possíveis.

Devido à este inconveniente com respeito às integrais, nós recorremos às distribuições a posteriori condicionais. Este tópico está relacionado com métodos MCMC (método de Monte Carlo com cadeias de Markov) e o amostrador de GIBBS que veremos mais adiante. As distribuições As *posterioris} condicionais são facilmente obtidas:

\(\checkmark\) A condicional de \(\sigma^2\) dado \(\beta_0\), \(\beta_1\) e os dados:

\[ \begin{array}{lll} &&f(\sigma^2 \mid \beta_0, \beta_1, \boldsymbol{x}, \boldsymbol{y}) \propto (\sigma^2)^ {-(b+\frac{n}{2}+1)} \exp \left[-\frac{1}{\sigma^2}\left(d+\frac{1}{2}\sum \limits_{i=1}^n (y_i - \beta_0 -\beta_1 x_i)^2\right) \right],\\ \text{ou seja, } &&\sigma^2 \mid \beta_0, \beta_1, \boldsymbol{x}, \boldsymbol{y} \sim IG \left(b+\frac{n}{2},d+\frac{1}{2}\sum \limits_{i=1}^n (y_i - \beta_0 -\beta_1 x_i)^2\right). \end{array} \]

a idéia de condicional nos diz que \(\beta_0\) e \(\beta_1\) são tratadas como constantes com respeito a \(\sigma^2\).

\(\checkmark\) A condicional de \(\beta_0\) dado \(\beta_1\), \(\sigma^2\) e os dados:

\[ \begin{array}{lll} &&f(\beta_0 \mid \beta_1, \sigma^2, \boldsymbol{x}, \boldsymbol{y}) \propto \exp\left[-\frac{\beta_0^2}{2a_0^2}\right] \exp\left[-\frac{1}{2\sigma^2}\sum \limits_{i=1}^n (y_i-\beta_0-\beta_1x_i)^2\right]\\ && \textbf{ Ideia: }\text{expandir os termos e identificar uma distribuição normal! Tome } \mu_i^{(1)} = y_i - \beta_1 x_i\\ &&\propto \exp\left[-\frac{\beta_0^2}{2a_0^2}\right] \exp\left[-\frac{1}{2\sigma^2}\sum \limits_{i=1}^n (\beta_0-\mu_i^{(1)})^2\right] \\ &&\propto \exp\left[-\frac{\beta_0^2}{2a_0^2}\right] \exp\left[-\frac{1}{2\sigma^2} \left(n \beta_0^2 + 2\beta_0 \sum \limits_{i=1}^n \mu_i^{(1)} + \sum \limits_{i=1}^n \mu_i^{{(1)}^2}\right)\right] \\ &&\propto \exp\left[-\frac{\beta_0^2}{2} \left(\frac{1}{a_0^2}+\frac{n}{\sigma^2}\right)+\frac{\beta_0 \sum \limits_{i=1}^n \mu_i^{(1)} }{\sigma^2}\right], \text{ o termo } \sum \limits_{i=1}^n \mu_i^{{(1)}^2} \text{ "caiu" pois é constante com respeito a } \beta_0\\ && \propto \exp\left[-\frac{1}{2\left(\frac{1}{a_0^2}+\frac{n}{\sigma^2}\right)^{-1}} \left(\beta_0^2-\frac{2\beta_0 \sum \limits_{i=1}^n \mu_i^{(1)} }{\sigma^2 \left(\frac{1}{a_0^2}+\frac{n}{\sigma^2}\right)}\right)\right]. \text{ Mas } \left(\frac{1}{a_0^2}+\frac{n}{\sigma^2}\right)^{-1} = \left(\frac{\sigma^2+a_0^2n}{a_0^2 \sigma^2}\right)^{-1}= \left(\frac{a_0^2 \sigma^2}{\sigma^2+a_0^2n}\right)\\ && \propto \exp\left[-\frac{1}{2\left(\frac{a_0^2 \sigma^2}{\sigma^2+a_0^2n}\right)} \left(\beta_0^2-\frac{2\beta_0 \sum \limits_{i=1}^n \mu_i^{(1)} }{\sigma^2 \left(\frac{a_0^2 \sigma^2}{\sigma^2+a_0^2n}\right)}\right)\right], \text{ e simplificando um pouco mais: }\\ && \propto \exp\left[-\frac{1}{2\left(\frac{a_0^2 \sigma^2}{\sigma^2+a_0^2n}\right)} \left(\beta_0^2-\frac{2\beta_0 a_0^2 \sum \limits_{i=1}^n \mu_i^{(1)} }{ \sigma^2+a_0^2n}\right)\right] \text{ e completando quadrados: }\\ && \propto \exp\left[-\frac{1}{2\left(\frac{a_0^2 \sigma^2}{\sigma^2+a_0^2n}\right)} \left(\beta_0^2-\frac{2\beta_0 a_0^2 \sum \limits_{i=1}^n \mu_i^{(1)} }{ \sigma^2+a_0^2n} + \left(\frac{ a_0^2 \sum \limits_{i=1}^n \mu_i^{(1)} }{ \sigma^2+a_0^2n}\right)^2 \right)\right]\\ && \propto \exp\left[-\frac{1}{2\left(\frac{a_0^2 \sigma^2}{\sigma^2+a_0^2n}\right)} \left(\beta_0-\frac{ a_0^2 \sum \limits_{i=1}^n \mu_i^{(1)} }{ \sigma^2+a_0^2n}\right)^2\right]\\ && \text{ou seja, } \beta_0 \mid \beta_1, \sigma^2, \boldsymbol{x}, \boldsymbol{y} \sim N\left(\frac{ a_0^2 \sum \limits_{i=1}^n \mu_i^{(1)}}{ \sigma^2+a_0^2n},\frac{a_0^2 \sigma^2}{\sigma^2+a_0^2n}\right) \end{array} \]

\(\checkmark\) A condicional de \(\beta_1\) dado \(\beta_0\), \(\sigma^2\) e os dados: Note que o parâmetro \(\beta_1\) acompanha o termo \(x_i\):

\[ \begin{array}{lll} &&f(\beta_1 \mid \beta_0, \sigma^2, \boldsymbol{x}, \boldsymbol{y}) \propto \exp\left[-\frac{\beta_1^2}{2a_1^2}\right] \exp\left[-\frac{1}{2\sigma^2}\sum \limits_{i=1}^n (y_i-\beta_0-\beta_1x_i)^2\right]. \text{ Tome } \mu_i^{(2)} = y_i - \beta_0 \\ && \propto \exp\left[-\frac{\beta_1^2}{2a_1^2}\right] \exp\left[-\frac{1}{2\sigma^2}\sum \limits_{i=1}^n (\beta_1 x_i -\mu_i^{(2)})^2\right]\\ && \propto \exp\left[-\frac{\beta_1^2}{2a_1^2}\right] \exp\left[-\frac{1}{2\sigma^2} \left(\sum \limits_{i=1}^n \beta_1^2 x_i^2 -2\beta_1\sum \limits_{i=1}^n \mu_i^{(2)} x_i + \sum \limits_{i=1}^n \mu_i^{(2)^2} \right) \right]\\ && \propto \exp\left[-\frac{\beta_1^2}{2a_1^2}\right] \exp\left[-\frac{1}{2\sigma^2} \left(\sum \limits_{i=1}^n \beta_1^2 x_i^2 -2\beta_1\sum \limits_{i=1}^n \mu_i^{(2)} x_i \right) \right]\\ && \propto \exp\left[-\frac{\beta_1^2}{2}\left(\frac{1}{a_1^2}+\frac{\sum \limits_{i=1}^n x_i^2}{\sigma^2}\right)+\frac{\beta_1\sum \limits_{i=1}^n \mu_i^{(2)} x_i}{\sigma^2}\right]\\ && \propto \exp\left[-\frac{1}{2 \left(\frac{1}{a_1^2}+\frac{\sum \limits_{i=1}^n x_i^2}{\sigma^2}\right)^{-1}}\left(\beta_1^2-\frac{2\beta_1\sum \limits_{i=1}^n \mu_i^{(2)} x_i}{\sigma^2 \left(\frac{1}{a_1^2}+\frac{\sum \limits_{i=1}^n x_i^2}{\sigma^2}\right)}\right)\right] \text{ Mas } \left(\frac{1}{a_1^2}+\frac{\sum \limits_{i=1}^n x_i^2}{\sigma^2}\right)^{-1}=\frac{a_1^2 \sigma^2}{\sigma^2+a_1^2 \sum \limits_{i=1}^n x_i^2}\\ &&\propto \exp\left[-\frac{1}{2 \left(\frac{a_1^2 \sigma^2}{\sigma^2+a_1^2 \sum \limits_{i=1}^n x_i^2}\right)}\left(\beta_1^2-\frac{2\beta_1\sum \limits_{i=1}^n \mu_i^{(2)} x_i}{\sigma^2 \left(\frac{a_1^2 \sigma^2}{\sigma^2+a_1^2 \sum \limits_{i=1}^n x_i^2}\right)}\right)\right]\\ &&\propto \exp\left[-\frac{1}{2 \left(\frac{a_1^2 \sigma^2}{\sigma^2+a_1^2 \sum \limits_{i=1}^n x_i^2}\right)}\left(\beta_1^2-\frac{2 a_1^2 \beta_1\sum \limits_{i=1}^n \mu_i^{(2)} x_i}{\sigma^2+a_1^2 \sum \limits_{i=1}^n x_i^2}\right)\right]\\ &&\propto \exp\left[-\frac{1}{2 \left(\frac{a_1^2 \sigma^2}{\sigma^2+a_1^2 \sum \limits_{i=1}^n x_i^2}\right)}\left(\beta_1^2-\frac{2 a_1^2 \beta_1\sum \limits_{i=1}^n \mu_i^{(2)} x_i}{\sigma^2+a_1^2 \sum \limits_{i=1}^n x_i^2}+\left( \frac{a_1^2 \sum \limits_{i=1}^n \mu_i^{(2)} x_i}{\sigma^2+a_1^2 \sum \limits_{i=1}^n x_i^2}\right)^2\right)\right]\\ &&\propto \exp\left[-\frac{1}{2 \left(\frac{a_1^2 \sigma^2}{\sigma^2+a_1^2 \sum \limits_{i=1}^n x_i^2}\right)}\left(\beta_1-\frac{a_1^2 \sum \limits_{i=1}^n \mu_i^{(2)} x_i}{\sigma^2+a_1^2 \sum \limits_{i=1}^n x_i^2}\right)^2\right]\\ && \text{ou seja, } \beta_1 \mid \beta_0, \sigma^2, \boldsymbol{x}, \boldsymbol{y} \sim N\left(\frac{a_1^2 \sigma^2}{\sigma^2+a_1^2 \sum \limits_{i=1}^n x_i^2},\frac{a_1^2 \sum \limits_{i=1}^n \mu_i^{(2)} x_i}{\sigma^2+a_1^2 \sum \limits_{i=1}^n x_i^2}\right) \end{array} \]

\(\checkmark\) Por fim, a condicional de \(\sigma^2\) dado \(\beta_1\), \(\beta_0\) e os dados:

\[ \begin{array}{lll} f(\sigma^2 \mid \beta_0,\beta_1,\boldsymbol{x}, \boldsymbol{y}) &\propto& (\sigma^2)^{-(b+\frac{n}{2}+1)} \exp\left[-\frac{d}{\sigma^2}\right] \exp\left[-\frac{1}{2\sigma^2}\sum \limits_{i=1}^n (y_i-\beta_0-\beta_1x_i)^2\right]\\ && \text{ou seja, } \sigma^2 \mid \beta_0,\beta_1,\boldsymbol{x}, \boldsymbol{y} \sim \text{ IG } \left(b+\frac{n}{2},d+\frac{1}{2}\sum \limits_{i=1}^n (y_i - \beta_0 -\beta_1 x_i)^2\right) \end{array} \]

• Esta metodologia de encontrar as posterioris condicionais para os coeficientes da regressão se extende ao modelo de regressão linear múltipla de maneira análoga.
• Uma outra alternativa a este problema é utilizar uma priori conjunta conjugada.

Aplicação da Metodologia

• Atribuir prioris não informativas para \(\beta_0\) e \(\beta_1\): Normais com média igual a zero e variância grande, por exemplo \(10^6\).
  
• Atribuir prioris não informativas para \(\sigma^2\). A distribuição \(IG (0.001,0.001)\) é não informativa, e é muito próxima à priori de Jeffreys para o modelo normal com média e variância desconhecidos.