1 Modelo
Verosimilitud: \[ y_i\mid\theta_i\stackrel{\text {iid}}{\sim}\textsf{Poison}(\theta_i)\,,\quad i = 1,\ldots,n \] donde \[ \log(\theta_i) = \boldsymbol{\beta}^{\textsf{T}}\boldsymbol{x}_i = \sum_{j=1}^k\beta_j\, x_{i,j} \] con \(\boldsymbol{\beta}=(\beta_1,\ldots,\beta_k)\) y \(\boldsymbol{x} = (x_{i,1},\ldots,x_{i,k})\).
Observe que la función de enlace (link function) en este modelo lineal generalizado es la función \(\log\).
Distribuciones previas conjugadas para modelos lineales generalizados no existen (a excepción del modelo de regresión Normal).
Distribución previa: \[ \boldsymbol{\beta}\sim\textsf{N}(\boldsymbol{\beta}_0,\Sigma_0) \]
Los parámetros del modelo son \((\beta_1,\ldots,\beta_k)\).
Los hiper-parámetros del model son \((\boldsymbol{\beta}_0,\Sigma_0)\).
2 Distibución conjunta posterior
\[ \begin{align*} p(\boldsymbol{\beta}\mid\boldsymbol{y}) &\propto p(\boldsymbol{y}\mid\boldsymbol{\beta})\,p(\boldsymbol{\beta}) \\ &= \prod_{i=1}^n \textsf{Poison}(y_i\mid\theta_i) \times\textsf{N}(\boldsymbol{\beta}_0,\Sigma_0) \\ &= \prod_{i=1}^n \frac{e^{-\theta_i}\,\theta_i^{y_i}}{y_i!} \times \exp\left\{ -\tfrac{1}{2}(\boldsymbol{\beta}-\boldsymbol{\beta}_0)^{\textsf{T}}\Sigma_0^{-1}(\boldsymbol{\beta}-\boldsymbol{\beta}_0) \right\} \\ &\propto \prod_{i=1}^n e^{-\theta_i}\,\theta_i^{y_i} \times \exp\left\{ -\tfrac{1}{2} \left[ \boldsymbol{\beta}^{\textsf{T}}\Sigma_0^{-1}\boldsymbol{\beta} - 2\boldsymbol{\beta}^{\textsf{T}}\Sigma_0^{-1}\boldsymbol{\beta}_0 \right] \right\} \end{align*} \] con \(\log(\theta_i) = \boldsymbol{\beta}^{\textsf{T}}\boldsymbol{x}_i\).
Por lo tanto, en escala logaritmica, se tiene que \[ \log p(\boldsymbol{\beta}\mid\boldsymbol{y}) \propto \sum_{i=1}^n(y_i\log(\theta_i) - \theta_i) -\tfrac{1}{2} \left[ \boldsymbol{\beta}^{\textsf{T}}\Sigma_0^{-1}\boldsymbol{\beta} - 2\boldsymbol{\beta}^{\textsf{T}}\Sigma_0^{-1}\boldsymbol{\beta}_0 \right] \]
Observe que \(p(\boldsymbol{\beta}\mid\boldsymbol{y})\) no corresponde a ninguna familia paramétrica de distribuciones conocida.
3 Algorithmo de Metropolis-Hastings
El algoritmo Metropolis-Hastings es un método genérico de aproximación de cualquier distribución posterior.
El problema radica cuando \(p(\boldsymbol{\theta}\mid \boldsymbol{y})\) no tiene una distribución estándar conocida de la cual sea posible simular fácilmente.
3.1 Metropolis
Algoritmo:
Dado un estado actual \(\theta^{(b)}\) del parámetro \(\theta\), se genera el siguiente estado \(\theta^{(b+1)}\) como sigue:
- Simular \(\theta^*\sim J(\theta\mid\theta^{(b)})\) donde \(J\) es una distribución simétrica: \(J(\theta_1\mid\theta_2) = J(\theta_2\mid\theta_1)\). Esta distribución se denomina distribución de propuestas (proposal distribution). Típicamente, \[ J(\theta\mid\theta^{(b)}) = \textsf{N}(\theta\mid\theta^{(b)},\delta^2)\qquad\text{o}\qquad J(\theta\mid\theta^{(b)}) = \textsf{U}(\theta\mid\theta^{(b)} - \delta, \theta^{(b)} + \delta) \] donde \(\delta\), conocido como parámetro de ajuste, se escoge de forma que el algoritmo funcione eficientemente.
- Calcular la tasa de aceptación: \[ r = \frac{p(\theta^*\mid y)}{p(\theta^{(b)}\mid y)}\,. \] Para brindar estabilidad numérica en el procesamiento, esta tasa usualmente se calcula primero en escala logarítmica: \[ r = \exp\left( \log p(\theta^*\mid y) - \log p(\theta^{(b)}\mid y) \right)\,. \]
- Establecer: \[ \theta^{(b+1)} = \left\{ \begin{array}{ll} \theta^*\, , & \hbox{con probabilidad $\min(1,r)$;} \\ \theta^{(b)}\, , & \hbox{con probabilidad $1-\min(1,r)$.} \end{array} \right. \]
Observaciones:
- Este algoritmo produce una secuencia de valores \(\theta^{(1)},\ldots,\theta^{(B)}\), donde \(\theta^{(b+1)}\) depende de los valores previos únicamente a través de \(\theta^{(b)}\), y por lo tanto la secuencia es una Cadena de Markov.
- La demostración de por qué funciona este algoritmo se base en el Teorema Ergodico para cadenas de Markov.
- La correlación de la cadena depende del valor de \(\delta\) en la distribución de propuestas. Es una práctica compun elegir un valor de \(\delta\) tal que la tasa de aceptación esté entre 30% y 50%. Existen algoritmos adaptativos para elegir automaticamente el valor de \(\delta\).
3.2 Metropolis-Hastings
El algoritmo Metrópolis-Hastings generaliza los algoritmos de Gibbs y Metropolis al permitir una distribución de propuesta arbitraria:
Algoritmo:
Dado un estado actual \(\theta^{(b)}\) del parámetro \(\theta\), se genera el siguiente estado \(\theta^{(b+1)}\) como sigue:
Simular \(\theta^*\sim J(\theta\mid\theta^{(b)})\) donde \(J\) es una distribución de propuestas arbitraria.
Calcular la tasa de aceptación: \[ r = \frac{p(\theta^*\mid y)}{p(\theta^{(b)}\mid y)}\,\frac{J(\theta^{(b)}\mid\theta^*)}{J(\theta^*\mid \theta^{(b)})}\,. \]
Establecer: \[ \theta^{(b+1)} = \left\{ \begin{array}{ll} \theta^*\, , & \hbox{con probabilidad $\min(1,r)$;} \\ \theta^{(b)}\, , & \hbox{con probabilidad $1-\min(1,r)$.} \end{array} \right. \]
Observaciones:
- Cada paso del algoritmo de Gibbs se puede ver como la generación de una propuesta de una distribución condicional completa y luego aceptarla con probabilidad 1.
- El algoritmo de Gibbs y el algoritmo de Metropolis son casos particulares del algoritmo de Metropolis-Hastings.
4 Ejemplo: Modelo Normal
Algoritmo de Metropolis en un modelo Normal con varianza conocida.
# MODELO (sigma^2 conocido)
# 1. y_i | theta ~ Normal(theta, sigma^2)
# 2. theta ~ Normal(mu, tau^2)
# 3. sigma^2 = 1
# simular datos
n <- 5
s2 <- 1
set.seed(1)
y <- round(rnorm(n,10,1),2)
# previa
t2 <- 10
mu <- 5
# posterior
mu.n <- ( mean(y)*n/s2 + mu/t2 )/( n/s2 + 1/t2 )
t2.n <- 1/( n/s2 + 1/t2 )
######## Metropolis
theta <- 0 # valor inicial -> rnorm(1, mu, sqrt(t2))
delta <- 1.75 # parametro de ajuste
S <- 10000 # numero de muestras
THETA <- NULL # almacenamiento
# cadena
set.seed(1)
for (s in 1:S) {
# 1. propuesta
theta.star <- rnorm(1,theta,sqrt(delta))
# 2. tasa de aceptacion
log.r <- ( sum(dnorm(y,theta.star,sqrt(s2),log=TRUE)) + dnorm(theta.star,mu,sqrt(t2),log=TRUE) ) -
( sum(dnorm(y,theta, sqrt(s2),log=TRUE)) + dnorm(theta, mu,sqrt(t2),log=TRUE) )
# 3. actualizar
if (runif(1) < exp(log.r)) {
theta<-theta.star
}
# 4. almacenar
THETA <- c(THETA, theta)
}
######## fin MCMC# grafico
par(mfrow=c(1,2), mar=c(3,3,1,1),mgp=c(1.75,.75,0))
# cadena
# adelgazamiento de la cadena (reducir autocorrelacion)
# skeep<-seq(10,S,by=10)
skeep <- 1:S
plot(skeep,THETA[skeep],type="l",xlab="Iteración",ylab=expression(theta))
# histograma metropolis y posterior analitica
# se omiten las primeras 50 observaciones (periodo de calentamiento)
hist(THETA[-(1:50)],prob=TRUE,main="",xlab=expression(theta),ylab="Densidad",col="gray95",border="gray")
th <- seq(min(THETA),max(THETA),length=100)
lines(th, dnorm(th,mu.n,sqrt(t2.n)), col = 2, lty = 2, lwd=2)
ACR <- NULL # tasa de aceptaciones
ACF <- NULL # autocorrelaciones
THETAA <- NULL # muestras
for(delta2 in 2^c(-5,-1,1,5,7) ) {
# parametros iniciales
THETA <- NULL
S <- 10000
theta <- 0
acs <- 0 # tasa de aceptacion
# cadena
set.seed(1)
for(s in 1:S) {
# 1. propuesta
theta.star<-rnorm(1,theta,sqrt(delta2))
# 2. tasa de aceptacion
log.r <- sum( dnorm(y,theta.star,sqrt(s2),log=TRUE) - dnorm(y,theta,sqrt(s2),log=TRUE) ) +
dnorm(theta.star,mu,sqrt(t2),log=TRUE) - dnorm(theta,mu,sqrt(t2),log=TRUE)
# 3. actualizar
if(runif(1) < exp(log.r)) {
theta <- theta.star
acs <- acs + 1
}
# 4. almacenar
THETA <- c(THETA, theta)
}
# fin MCMC
# almacenar valores de todos los casos (delta2)
ACR <- c(ACR, acs/s)
ACF <- c(ACF,acf(THETA,plot=FALSE)$acf[2])
THETAA <- cbind(THETAA,THETA)
}# graficos
par(mfrow=c(1,3),mar=c(2.75,2.75,.5,.5),mgp=c(1.7,.7,0))
laby <- c(expression(theta),"","","","")
for(k in c(1,3,5)) {
plot(THETAA[1:500,k],type="l",xlab="Iteración",ylab=laby[k], ylim=range(THETAA) )
abline(h=mu.n,lty=2)
}
# tasas de aceptacion
ACR## [1] 0.8713 0.5681 0.3549 0.0980 0.0487# autocorrelaciones
ACF ## [1] 0.9750710 0.7724628 0.6904980 0.8401160 0.8636223# tamaños efectivos de muestra
coda::effectiveSize(THETAA)## THETA THETA THETA THETA THETA
## 116.2963 793.0188 1427.4509 612.6138 307.32005 Ejemplo: Regresión Poisson
Actividades de reproducción de gorriones en función de la edad (Arcese et al., 1992).
Caracterizar el número de crias en función de la edad.
- \(y_i\,\,\,\,\): número de crias, para \(i=1,\ldots,n\).
- \(x_{i,j}\,\): predictor \(j\) observado en el individuo \(i\), para \(i=1,\ldots,n\) y \(j=1,\ldots,k\).
#-------------------------------------------------------------------------------
# Descripcion:
# Actividades de reproduccion de gorriones en funcion de la edad (Arcese et al, 1992).
# n = 52 gorriones hembras.
# "age" : edad.
# "fledged" : numero de crias.
#-------------------------------------------------------------------------------
############
### data ###
############
spfage <- structure(c(3, 1, 1, 2, 0, 0, 6, 3, 4, 2, 1, 6, 2, 3, 3, 4, 7, 2, 2, 1,
1, 3, 5, 5, 0, 2, 1, 2, 6, 6, 2, 2, 0, 2, 4, 1, 2, 5, 1, 2,
1, 0, 0, 2, 4, 2, 2, 2, 2, 0, 3, 2, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2,
2, 2, 2, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 4, 4,
4, 4, 5, 5, 5, 5, 3, 3, 3, 3, 3, 3, 3, 6, 1, 1, 9, 9, 1, 1,
1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 25, 25, 16, 16,
16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 25, 16, 16, 16, 16,
25, 25, 25, 25, 9, 9, 9, 9, 9, 9, 9, 36, 1, 1),
.Dim = c(52L, 4L),
.Dimnames = list(NULL, c("fledged", "intercept", "age", "age2")))
spfage <- as.data.frame(spfage)
spf <- spfage$fledged # y = variable dependiente (respuesta)
age <- spfage$age # x1 = variable independiente 1
age2 <- age^2 # x2 = variable independiente 2############################
### analisis descriptivo ###
############################
par(mar=c(3,3,1,1),mgp=c(1.75,.75,0))
plot(spf~as.factor(age), range=0, xlab="Edad (años)", ylab="No. Crias", col="gray", border="lightgray")
# GLM (frecuentista)
summary(glm(spf~age+age2,family="poisson"))##
## Call:
## glm(formula = spf ~ age + age2, family = "poisson")
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4650 -0.6355 -0.2298 0.4937 2.0429
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.27662 0.44219 0.626 0.5316
## age 0.68174 0.33850 2.014 0.0440 *
## age2 -0.13451 0.05786 -2.325 0.0201 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 76.081 on 51 degrees of freedom
## Residual deviance: 67.837 on 49 degrees of freedom
## AIC: 198.78
##
## Number of Fisher Scoring iterations: 5###################
### Monte Carlo ###
###################
y <- spf # variable respuesta
X <- cbind(rep(1,length(y)),age,age^2) # matriz de diseño
n <- length(y) # tamaño de la muestrsa
p <- dim(X)[2] # numero de predictores
# previa
pmn.beta <- rep(0, p) # beta0 = 0
psd.beta <- rep(10, p) # Sigma0 = 10*I
#-------------------------------------------------------------------------------
# En muchos problemas, la varianza posterior es una elección eficiente para la
# distribución de propuestas. Aunque no se conoce la varianza posterior antes de
# ejecutar el algoritmo de Metrópolis, a menudo basta con utilizar una aproximación.
# Si esto da como resultado una tasa de aceptación demasiado alta o demasiado baja,
# siempre es posible ajustar la variabilidad de la propuesta en consecuencia.
#-------------------------------------------------------------------------------
# log: funcion de enlace
# y + 1/2 evitar problemas en la frontera con 0
var.prop <- var(log(y + 1)) * solve( t(X)%*%X ) # matriz de varianza propuesta
beta <- rep(0,p) # valor inicial beta
######## algoritmo de metropolis
S <- 25000
BETA <- matrix(NA, nrow=S, ncol=p)
ac <- 0
ncat <- floor(S/10)
######## functiones auxiliares
rmvnorm <- function(n,mu,Sigma)
{
p<-length(mu)
res<-matrix(0,nrow=n,ncol=p)
if( n>0 & p>0 ) {
E<-matrix(rnorm(n*p),n,p)
res<-t( t(E%*%chol(Sigma)) +c(mu))
}
res
}
# cadena
set.seed(1)
for(s in 1:S) {
# 1. propuesta
beta.p<- t(rmvnorm(1, beta, var.prop ))
# 2. tasa de aceptacion
lhr <- sum(dpois(y,exp(X%*%beta.p),log=T)) -
sum(dpois(y,exp(X%*%beta),log=T)) +
sum(dnorm(beta.p,pmn.beta,psd.beta,log=T)) -
sum(dnorm(beta, pmn.beta,psd.beta,log=T))
# 3. actualizar
if (runif(1) < exp(lhr)) {
beta <- beta.p
ac <- ac + 1
}
# 4. almacenar
BETA[s,] <- beta
# 5. Progreso
if (s%%ncat == 0) cat(100*round(s/S, 1), "% completed ... \n", sep = "" )
}## 10% completed ...
## 20% completed ...
## 30% completed ...
## 40% completed ...
## 50% completed ...
## 60% completed ...
## 70% completed ...
## 80% completed ...
## 90% completed ...
## 100% completed ...######### fin mcmc
# tasa de aceptacion
100*ac/S## [1] 53.464# diagnosticos
library(coda)
apply(X = BETA, MARGIN = 2, FUN = effectiveSize)## [1] 1663.602 1575.657 1451.526#### grafico diagnostico
blabs <- c(expression(beta[1]),expression(beta[2]),expression(beta[3])) # etiquetas
thin <- c(1,(1:1000)*(S/1000)) # muestreo sistematico
par(mfrow=c(1,3),mar=c(2.75,2.75,.5,.5),mgp=c(1.7,.7,0))
j<-3
plot(thin,BETA[thin,j],type="l",xlab="Iteración",ylab=blabs[j])
abline(h=mean(BETA[,j]) )
acf(BETA[,j],ci.col="gray",xlab="lag")
acf(BETA[thin,j],xlab="lag/10",ci.col="gray")
####par(mfrow=c(1,3), mar=c(2.75,2.75,.5,.5),mgp=c(1.7,.7,0))
plot(density(BETA[,2],adj=2), type="l", main="",
xlab=expression(beta[2]),
ylab=expression(paste(italic("p("),beta[2],"|",italic("y)"),sep="") ) ,
lwd=2,lty=1,col="black")
plot(density(BETA[,3],adj=2), type="l", main="",
xlab=expression(beta[3]),
ylab=expression(paste(italic("p("),beta[3],"|",italic("y)"),sep="") ),
lwd=2,col="black",lty=1)
Xs <- cbind(rep(1,6),1:6,(1:6)^2)
eXB.post <- exp(t(Xs%*%t(BETA )) )
qE <- apply( eXB.post,2,quantile,probs=c(.025,.5,.975))
plot(c(1,6),range(c(0,qE)),type="n",xlab="Edad (años)", ylab="No. Crias")
lines( qE[1,],col="black",lwd=1)
lines( qE[2,],col="black",lwd=2)
lines( qE[3,],col="black",lwd=1)
quantile(BETA[,2], c(.025, .975))## 2.5% 97.5%
## 0.0719287 1.4414667quantile(BETA[,3], c(.025, .975))## 2.5% 97.5%
## -0.26440202 -0.03074222mean( BETA[,2] > 0 )## [1] 0.9846mean( BETA[,3] > 0 )## [1] 0.00468########
# JAGS #
########
library(R2jags)
model <- function() {
# verosimilitud
for (i in 1:n) {
y[i] ~ dpois(theta[i])
log(theta[i]) <- inprod(X[i,], beta) # X[i,1]*beta[1]+X[i,2]*beta[2]+X[i,3]*beta[3]
}
# previa
for (j in 1:p) {
beta[j] ~ dnorm(beta0, phi0)
}
}
# previa
beta0 <- 0
phi0 <- 1/10
# input
model_data <- list(y = y, X = X, n = n, p = p, beta0 = beta0, phi0 = phi0)
# parameters
model_parameters <- c("beta")
# initial values
initial_values <- list(list("beta" = rep(beta0, p)),
list("beta" = rep(beta0, p)),
list("beta" = rep(beta0, p)))
# mcmc settings
niter <- 26000
nburn <- 1000
nthin <- 25
nchain <- length(initial_values)
# mcmc
set.seed(123)
fit <- jags(data = model_data, inits = initial_values,
parameters.to.save = model_parameters, model.file = model,
n.chains = nchain, n.iter = niter, n.thin = nthin, n.burnin = nburn)
print(fit)
# transformar a objecto MCMC
fit_mcmc <- coda::as.mcmc(fit)
# cadena
dim(fit$BUGSoutput$sims.list$beta)
# plots
mcmcplots::traplot(fit_mcmc)
mcmcplots::denplot(fit_mcmc)