21 Paquete INLA
El paquete INLA (integrated nested Laplace approximation) se utiliza para hacer inferencia bayesiana aproximada de diversos modelos. Al visitar este enlace se encontrará un gran contenido que le podría servir para conocer más sobre INLA.
En el enlace de abajo se puede acceder a un libro muy interesante sobre modelos mixtos con INLA.
https://becarioprecario.bitbucket.io/inla-gitbook/index.html
Ejemplo: modelo Poisson con intercepto aleatorio
En este ejemplo vamos a simular \(n_i=50\) observaciones para \(G=10\) grupos (en total 500 obs) que tengan la estructura mostrada abajo. El objetivo del ejemplo es ilustrar el uso de la función inla para estimar los parámetros del modelo.
\[\begin{align*} y_{ij} | b_0 &\sim Poisson(\mu_{ij}) \\ \log(\mu_{ij}) &= 4 - 6 x_{ij} + b_{0i} \\ b_{0} &\sim N(0, \sigma^2_{b0}=4) \\ x_{ij} &\sim U(0, 1) \end{align*}\]
El vector de parámetros de este modelo es \(\boldsymbol{\Theta}=(\beta_0=4, \beta_1=-6, \sigma^2_{b0}=4)^\top\).
El código para simular las observaciones se muestra a continuación. Se fijó la semilla con set.seed(1234) para que el lector pueda replicar el ejemplo.
gen_data_1 <- function() {
ni <- 50
G <- 10
nobs <- ni * G
grupo <- factor(rep(x=1:G, each=ni))
obs <- rep(x=1:ni, times=G)
x <- runif(n=nobs, min=0, max=1)
b0 <- rnorm(n=G, mean=0, sd=sqrt(4)) # Intercepto aleatorio
b0 <- rep(x=b0, each=ni) # El mismo intercepto aleatorio pero repetido
media <- exp(4 - 6 * x + b0)
y <- rpois(n=nobs, lambda=media)
datos <- data.frame(obs, grupo, b0, x, y)
datos
}
set.seed(1234)
datos <- gen_data_1()
head(datos)## obs grupo b0 x y
## 1 1 1 0.8738612 0.1137034 74
## 2 2 1 0.8738612 0.6222994 2
## 3 3 1 0.8738612 0.6092747 6
## 4 4 1 0.8738612 0.6233794 3
## 5 5 1 0.8738612 0.8609154 1
## 6 6 1 0.8738612 0.6403106 0Para estimar los parámetros del modelo con INLA se usa el siguiente código.
Aplicando la función summary sobre el objeto fit1 se obtiene:
Time used:
Pre = 0.321, Running = 1.7, Post = 0.0978, Total = 2.12
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 4.402 0.555 3.292 4.403 5.509 NA 0
x -6.007 0.048 -6.102 -6.007 -5.913 NA 0
Random effects:
Name Model
grupo IID model
Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for grupo 0.40 0.175 0.138 0.374 0.787 NA
Marginal log-Likelihood: -1110.13
is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')De la salida anterior se tiene que \(\hat{\boldsymbol{\Theta}}=(\hat{\beta_0}=4.402, \hat{\beta_1}=-6.007, \hat{\sigma}^2_{b0}=(1/0.40)^2=2.5^2)^\top\) mientras que \(\boldsymbol{\Theta}=(\beta_0=4, \beta_1=-6, \sigma^2_{b0}=4)^\top\).
Ejemplo: modelo Poisson con intercepto y pendiente aleatoria
En este ejemplo vamos a simular \(n_i=50\) observaciones para \(G=10\) grupos (en total 500 obs) que tengan la estructura mostrada abajo. El objetivo del ejemplo es ilustrar el uso de la función inla para estimar los parámetros del modelo.
\[\begin{align*} y_{ij} | b_0, b_1 &\sim Poisson(\mu_{ij}) \\ \log(\mu_{ij}) &= 2 - 3 x_{ij} + b_{0i} + b_{1i} x_{ij} \\ \left ( \begin{matrix} b_{0} \\ b_{1} \end{matrix} \right ) &\sim N\left ( \left [ \begin{matrix} 0 \\ 0 \end{matrix} \right ], \left [ \begin{matrix} \sigma^2_{b0}=0.7 & \sigma_{b01}=0.5 \\ \sigma_{b01} & \sigma^2_{b1}=0.6 \end{matrix} \right ] \right ) \\ x_{ij} &\sim U(0, 1) \end{align*}\]
El vector de parámetros de este modelo es \(\boldsymbol{\Theta}=(\beta_0=2, \beta_1=-3, \sigma_{b0}^2=0.7, \sigma_{b1}^2=0.6, \sigma_{b01}=0.5)^\top\).
El código para simular las observaciones se muestra a continuación.
gen_data_2 <- function() {
ni <- 50
G <- 30
nobs <- ni * G
grupo <- factor(rep(x=1:G, each=ni))
obs <- rep(x=1:ni, times=G)
x <- runif(n=nobs, min=0, max=1)
Sigma <- matrix(c(0.7, 0.5, # Matriz de var-cov
0.5, 0.6), ncol=2, nrow=2)
b <- MASS::mvrnorm(n=G, mu=rep(0, 2), Sigma=Sigma) # Simulando b0 y b1
b <- apply(b, MARGIN=2, function(c) rep(c, each=ni)) # Replicando
b0 <- as.vector(b[, 1]) # Separando los b0
b1 <- as.vector(b[, 2]) # Separando los b1
media <- exp(2 - 3 * x + b0 + b1 * x)
y <- rpois(n=nobs, lambda=media)
datos <- data.frame(obs, grupo=as.factor(grupo), b0, b1, x, y)
datos
}
set.seed(1234)
datos <- gen_data_2()
head(datos)## obs grupo b0 b1 x y
## 1 1 1 0.4942715 -0.3878164 0.1137034 10
## 2 2 1 0.4942715 -0.3878164 0.6222994 2
## 3 3 1 0.4942715 -0.3878164 0.6092747 2
## 4 4 1 0.4942715 -0.3878164 0.6233794 3
## 5 5 1 0.4942715 -0.3878164 0.8609154 0
## 6 6 1 0.4942715 -0.3878164 0.6403106 1Debemos preparar la base de datos para que se pueda ajustar un modelo con intercepto y pendiente aleatoria. Para mayores detalles de la preparación se puede consultar Wang, Yue, and Faraway (2018) página 110.
n_groups <- length(levels(datos$grupo))
datos$numid <- as.numeric(datos$grupo)
datos$slopeid <- datos$numid + n_groupsPara estimar los parámetros del modelo con INLA se usa el siguiente código.
library(INLA)
formula <- y ~ x + f(numid, model="iid2d", n = 2*n_groups) +
f(slopeid, x, copy="numid")
fit2 <- inla(formula,
data = datos, family = "poisson",
control.predictor = list(compute = TRUE))Aplicando la función summary sobre el objeto fit2 se obtiene:
Time used:
Pre = 0.251, Running = 0.404, Post = 0.243, Total = 0.898
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 1.930 0.133 1.665 1.931 2.190 1.931 0
x -3.069 0.186 -3.453 -3.063 -2.717 -3.063 0
Random effects:
Name Model
numid IID2D model
slopeid Copy
Model hyperparameters:
mean sd 0.025quant 0.5quant
Precision for numid (component 1) 2.170 0.601 1.200 2.099
Precision for numid (component 2) 1.427 0.500 0.679 1.349
Rho1:2 for numid 0.688 0.122 0.402 0.705
0.975quant mode
Precision for numid (component 1) 3.549 1.970
Precision for numid (component 2) 2.623 1.207
Rho1:2 for numid 0.873 0.742
Marginal log-Likelihood: -2437.24
is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
De la salida anterior se tiene que \(\hat{\boldsymbol{\Theta}}=(\hat{\beta_0}=1.930, \hat{\beta_1}=-3.069, \hat{\sigma}^2_{b0}=1/2.170=0.4608, \hat{\sigma}^2_{b1}=1/1.427=0.7001, \hat{\sigma}_{b01}=0.688 * (1/2.170) * (1/1.427)=0.2222)^\top\) mientras que \(\boldsymbol{\Theta}=(\beta_0=2, \beta_1=-3, \sigma_{b0}^2=0.7, \sigma_{b1}^2=0.6, \sigma_{b01}=0.5)^\top\).