1. Introducción
En esta publicación se muestra como se estiman los parámetros de un modelo de Regresión Lineal Normal Bayesiana mediante el algoritmo de Metropolis-Hastings, también conocido como uno de los Métodos de Cadenas de Markov de Montecarlo (MCMC).
2. Regresión lineal bayesiana
Considerando \(\boldsymbol{\theta}=\left(\theta_{0}, \ldots, \theta_{k}\right)^{\top}\) como el vector de parámetros a estimar, \(\mathbf{y}=\left(y_{1}, \ldots, y_{n}\right)^{\top}\) un vector \(n \times 1\) donde es una realización de la variable aleatoria cuya distribución es \(p\left(y_{i} \mid \boldsymbol{\theta}\right)\). La función de verosimilitud es dado por \(y_{i}\):
\[\begin{equation} \mathcal{L}(\boldsymbol{\theta} \mid \mathbf{y})=\prod_{i=1}^{n} p\left(y_{i} \mid \boldsymbol{\theta}\right). \end{equation}\]
Toda la información de las observaciones \(y_{i}\) sobre \(\boldsymbol{\theta}\) se incluye en la ecuación anterior. La dificultad para estimar \(\boldsymbol{\theta}\) se vuelve un problema de optimización de maximizar la función de probabilidad (o su logaritmo). En contraste, bajo la metodología bayesiana la estimación de \(\boldsymbol{\theta}\) viene dada por la distribución posterior conjunta, que está definida por la Teorema de Bayes:
\[\begin{equation} p(\boldsymbol{\theta} \mid \mathbf{y})=\frac{\mathcal{L}(\theta \mid \mathbf{y}) p(\boldsymbol{\theta})}{\int_{\Theta} \mathcal{L}(\theta \mid \mathbf{y}) p(\boldsymbol{\theta}) d \boldsymbol{\theta}} \end{equation}\]
donde \(\Theta\) representa el espacio paramétrico de \(\boldsymbol{\theta}\) y \(p(\boldsymbol{\theta})\) es la distribución a priori. La ecuación anterior se puede escribir como:
\[\begin{equation} p(\boldsymbol{\theta} \mid \mathbf{y}) \propto \mathcal{L}(\theta \mid \mathbf{y}) p(\boldsymbol{\theta}) \end{equation}\]
donde \(\int_{\Theta} \mathcal{L}(\theta \mid \mathbf{y}) p(\boldsymbol{\theta}) d \boldsymbol{\theta}\) es la distribución marginal de \(\mathbf{y}\) que no depende de \(\boldsymbol{\theta}.\) La distribución a posteriori \(p(\boldsymbol{\theta} \mid \mathbf{y})\) proporciona toda la información que uno puede tener sobre \(\boldsymbol{\theta}\). Por ejemplo, es posible evaluar \(p(\boldsymbol{\theta} \mid \mathbf{y})\) y su media, mediana, varianza y algunas otras cantidades como cuantiles para tener estimaciones puntuales e intervalos. Además, la distribución a posteriori con frecuencia no tiene forma cerrada, por lo que va depender de los métodos computacionales para poder simular la distribución.
3. Modelo de regresión lineal
La regresión lineal es el modelo básico para el modelamiento estadístico que consiste en determinar relaciones de dependencia de tipo lineal entre una variable dependiente o respuesta, respecto de una o varias variables explicativas o independiente.
\[\begin{equation} y=\beta_1 + \beta_2 X_{2i} +\beta_3 X_{3i} +...+\beta_k X_{ki} + \epsilon_i \end{equation}\]
Donde, asumimos que el error aleatorio \(\epsilon_i\) es IID (Independiente e Idénticamente Distribuida) con media cero y varianza contante \(\sigma^2\), es decir,
\[\begin{equation} \epsilon_i \mathrel{\mathop{\sim}\limits^{\rm iid}}\textsf{Normal}(0, \sigma^2). \end{equation}\]
Para la estimación de los parámetro se hace uso del método de Mínimos Cuadrados Ordinarios (MCO), el método de Máxima Verosimilitud, entre otros. Para el ajuste del modelo se puede hacer uso de la función lmdel R.
Para nuestra aplicación consideremos el siguiente modelo de regresión lineal simple \(y_i=\beta_{0}+ \beta_{1}\,x_{i1}+\epsilon_i\), para \(i=1,\,2,\,3,\,\ldots\,n\), \(\epsilon_i\sim N(0,\sigma_\epsilon^2)\). Para fines de aplicación consideremos \(\beta_{0}=10\), \(\beta_{1}=5\) y \(\sigma_{\epsilon}=5\).
4. Simulando la data
#Fijando los parámetros
b1 <- 5
b0 <- 10
desv <- 5
sampleSize <- 31
# Creando la variable independiente
x <- (-(sampleSize-1)/2):((sampleSize-1)/2)
# Creando la variable dependiente considerando y=b0+b1*N(0,sd)
y <- b0 + b1 * x + rnorm(n=sampleSize,mean=0,sd=desv)
#Plot (x,y)
plot(x,y, main="Diagrama de dispersión de la data simulada",col="blue")
5. Función de Verosimilitud del modelo
Para estimar los parámetros en un modelo, necesitamos derivar la función de verosimilitud del modelo que queremos ajustar.
La función de verosimilitud del modelo de regresión lineal simple está dada por:
\[\begin{align} L(\beta_{0},\,\beta_{1},\,\sigma_\epsilon|Datos)&=\prod^{n}_{i=1}\frac{1}{\sigma_\epsilon\sqrt{2\pi}}\exp\left(-\frac{\left(y_i-\beta_{0}-\beta_{1}\,x_i\right)^2}{2\sigma^2_\epsilon}\right)\\ &=\frac{1}{{(2\pi\sigma^2_\epsilon)}^{n/2}}\exp\left(-\frac{\sum^n_{i=1}\left(y_i-\beta_{0}-\beta_{1}\,x_i\right)^2}{2\sigma^2_\epsilon}\right) \end{align}\]
y su logaritmo es:
\[\begin{align} \log\left(L(\beta_{0},\,\beta_{1},\,\sigma_\epsilon|Datos)\right)=-\frac{n}{2}\log\left(2\pi\sigma^2_\epsilon\right)-\frac{\sum^n_{i=1}\left(y_i-\beta_{0}-\beta_{1}\,x_i\right)^2}{2\sigma^2_\epsilon} \end{align}\]
#Definiendo la función de verosimilitud
likelihood <- function(param){
a = param[1]
b = param[2]
sd = param[3]
pred = a*x + b
singlelikelihoods = dnorm(y, mean = pred, sd = sd, log = T)
sumll = sum(singlelikelihoods)
return(sumll)
}
# Example: plot the likelihood profile of the slope a
slopevalues <- function(x){return(likelihood(c(x, b1, desv)))}
slopelikelihoods <- lapply(seq(3, 7, by=.05), slopevalues )
plot (seq(3, 7, by=.05), slopelikelihoods , type="l", xlab = "values of slope parameter a",
ylab = "Log likelihood",col="blue")
6. Definiendo la distribución a priori
En estadística Bayesiana tenemos que especificar una distribución a priori para cada parámetro. Para hacerlo más fácil, usaremos distribuciones normales para los dos parámetros de la ecucación, mientras que para la varianza usaremos la distribución gamma. Para la implementación Bayesiana del modelo se asumirán las siguientes distribuciones a priori para \(\beta_{0}\), \(\beta_{1}\) y \(\sigma^2_\epsilon\)
\[\begin{align} a\sim &\, N(\mu=5,\,\sigma^2=0.1^2)\\ b\sim &\, N(\mu=10,\,\sigma^2=0.1^2)\\ \sigma^2_\epsilon\sim &\, Gamma(\alpha=4,\,\beta=1.25) \end{align}\]
Denotando \(p(\beta_{0})\), \(p(\beta_{1})\) y \(p(\sigma^2_\epsilon)\) como las distribuciones a priori para \(\beta_{0}\), \(\beta_{1}\) y \(\sigma^2_\epsilon\), respectivamente. Se tiene que la distribución conjunta a priori está dada por:
\[\begin{align} p(\beta_{0},\,\beta_{1},\,\sigma^2_\epsilon)=p(\beta_{0})p(\beta_{1})p(\sigma_\epsilon^2) \end{align}\]
y su logaritmo es:
\[\begin{align} \log\left(p(\beta_{0},\,\beta_{1},\,\sigma^2_\epsilon)\right)=\log(p(\beta_{0}))+\log(p(\beta_{1}))+\log(p(\sigma_\epsilon^2)) \end{align}\]
#Distribución a priori
prior <- function(param){
a <- param[1]
b <- param[2]
sd <- param[3]
aprior <- dnorm(a,mean=a,sd=0.1,log=T)
bprior <- dnorm(b,mean=b,sd=0.1, log=T)
sdprior <- dgamma(sd, shape=2.5, scale=2, log=T)
return(aprior + bprior + sdprior)
}7. Definiendo la distribución a posteriori
Se sabe que por el teorema de Bayes, la distribución a posteriori tiene la siguiente forma:
\[\begin{align} p(a,\,b,\,\sigma^{2}_\epsilon|Datos)\propto L(a,\,b,\,\sigma_\epsilon|Datos)p(a,\,b,\,\sigma^2_\epsilon) \end{align}\]
y su logaritmo:
\[\begin{align} \log\left(p(a,\,b,\,\sigma^{2}_\epsilon|Datos)\right)\propto \log(L(a,\,b,\,\sigma_\epsilon|Datos))+\log(p(a))+\log(p(b))+\log(p(\sigma_\epsilon^2)) \end{align}\]
#Distribución a posteriori
posterior <- function(param){
return (likelihood(param) + prior(param))
}8. Métodos de Cadenas de Markov de Montecarlo (MCMC)
Ahora explicaremos el algoritmo de Metropolis-Hastings. Una de las aplicaciones más frecuentes de este algoritmo (como en este ejemplo) es la simulación de la distribución a posteriori. En principio, el algoritmo puede usarse para muestrear desde cualquier función integrable. Entonces, el objetivo de este algoritmo es recorrer el espacio paramétrico, pero de una manera que la probabilidad de estar en un punto sea proporcional a la función de la que tomamos muestras (esto generalmente se llama función objetivo).En nuestro caso, este es la distribución a posteriori definido anteriormente.
Iniciamos con un valor de parámetro aleatorio.
Elegimos un nuevo valor de parámetro cercano al valor anterior basado en alguna densidad de probabilidad que se llama función de propuesta.
Pasar a este nuevo punto con una probabilidad \(p(new)/p(old)\), donde \(p\) es la función objetivo.
Tenga en cuenta que tenemos una distribución simétrica propuesta es \(q(x'|x)\).
\[\begin{align} q(x'|x) \sim N(x, c(0.1,0.5,0.3) ) \end{align}\]
La desviación estándar \(\sigma\) con fijas. Las distribuciones son simétricas \(q(x'|x) = q(x|x')\). Por tanto, la probabilidad de aceptación es igual a:
\[\begin{align} \alpha(x'|x)=\min \left(1,\frac{P(x')}{P(x)} \right) \end{align}\]
#Metropolis algorithm
proposalfunction <- function(param){
return(rnorm(3,mean = param, sd= c(0.1,0.5,0.3)))
}
run_metropolis_MCMC <- function(startvalue, iterations){
chain = array(dim = c(iterations+1,3))
chain[1,] = startvalue
for (i in 1:iterations){
proposal = proposalfunction(chain[i,])
probab = exp(posterior(proposal) - posterior(chain[i,]))
if (runif(1) < probab){
chain[i+1,] = proposal
}else{
chain[i+1,] = chain[i,]
}
}
return(chain)
}Trabajar con los logaritmos de la distribución a posteriori puede ser un poco confuso al principio, en particular cuando se observa la línea donde se calcula la probabilidad de aceptación.
9. Implementación
En este punto ejecutaremos la función implementada del algoritmo de Metropolis-Hastings y luego procederemos con la estimación de los parámetros.
startvalue = c(4,0,10)
burnIn = 5000
chain = run_metropolis_MCMC(startvalue, 55000)#Estimación del parámetro b0
mean(chain[-(1:burnIn),2])## [1] 9.832702#Estimación del parámetro b1
mean(chain[-(1:burnIn),1])## [1] 5.070862#Estimación del parámetro desv
mean(chain[-(1:burnIn),3])## [1] 5.919968A continuación obtendremos el valor aceptación e interpretaremos dicho valor.
acceptance = 1-mean(duplicated(chain[-(1:burnIn),]))
acceptance## [1] 0.6399072Los primeros pasos del algoritmo pueden estar sesgados por el valor inicial y, por lo tanto, generalmente se descartan para el análisis de la distribución a posteriori (burnin). Un resultado interesante a considerar es la tasa de aceptación: ¿con qué frecuencia se rechazó una propuesta por el criterio de aceptación de Metropolis-Hastings?. La tasa de aceptación puede verse influenciada por la función de propuesta: generalmente, cuanto más cercanas están las propuestas, mayor es la tasa de aceptación. Sin embargo, las tasas de aceptación muy altas no suelen ser beneficiosas, pues significa que los algoritmos se “quedan” en el mismo punto, lo que da como resultado un muestreo no óptimo del espacio paramétrico.
Finalmente, podemos graficar los resultados y analizaremos la convergencia de las cadenas.
par(mfrow = c(2,3))
hist(chain[-(1:burnIn),1],nclass=30, main="Posterior of a", xlab="True value = red line" )
abline(v = mean(chain[-(1:burnIn),1]))
abline(v = b1, col="red" )
hist(chain[-(1:burnIn),2],nclass=30, main="Posterior of b", xlab="True value = red line")
abline(v = mean(chain[-(1:burnIn),2]))
abline(v = b0, col="red" )
hist(chain[-(1:burnIn),3],nclass=30, main="Posterior of sd", xlab="True value = red line")
abline(v = mean(chain[-(1:burnIn),3]))
abline(v = desv, col="red")
plot(chain[-(1:burnIn),1], type = "l", xlab="True value = red line" , main = "Chain values of a",)
abline(h = b1, col="red")
plot(chain[-(1:burnIn),2], type = "l", xlab="True value = red line" , main = "Chain values of b",)
abline(h = b0, col="red")
plot(chain[-(1:burnIn),3], type = "l", xlab="True value = red line" , main = "Chain values of sd",)
abline(h = desv, col="red")
Observando los gráficos anteriores podemos concluir que las cadenas son estacionarias, lo cual es un requisito indispensable en el análisis bayesiano.
Para comparar los resultados de la estimación de los parámetros usaremos las estimaciones bajo el enfoque clásico.
summary(lm(y~x))##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.707 -3.508 -0.911 4.743 11.159
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.8083 1.0435 9.399 2.64e-10 ***
## x 5.0702 0.1167 43.456 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.81 on 29 degrees of freedom
## Multiple R-squared: 0.9849, Adjusted R-squared: 0.9844
## F-statistic: 1888 on 1 and 29 DF, p-value: < 2.2e-1610. Conclusión
- Según las estimaciones obtenidas pudimos recuperar más o menos los parámetros originales que se utilizaron para simular nuestros datos, y también se verifica que estas estimaciones son muy aproximadas a las obtenidas por el enfoque clásico. Para una mejor aproximación siempre se recomienda incrementar el tamaño de muestra y el número de iteraciones.