Evaluación y Comparación de Modelos

Objetivos:

Usar posterior predictive checks para evaluar si los modelos que ajustamos son apropiados para nuestros datos
Comprar modelos alternativos

Ejemplo: independencia en una secuencia de tipo presencia/ausencia

Colectamos datos en una transecta donde registramos cada \(10\) metros la presencia (\(1\)) o ausencia (\(0\)) de una especie de interés en parcelas de \(1\) metro cuadrado y obtenemos el siguiente resultado:

y <- c(1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0)

Como una primera aproximación, podríamos pensar que un modelo apropiado para estos datos sería una Binomial con probabilidad de éxito (presencia) fijo. Como vimos en el práctico de análisis Bayesianos, existe una solución analítica para la posterior de la probabilidad de éxito de una Binomial si usamos como previa una distribución Beta. Usando una previa no-informativa tenemos:

op <- par(cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, las = 1, bty = "n")
curve(dbeta(x, 1 + sum(y), 1 + length(y) - sum(y)), lwd = 2, xlab = "Probabilidad de presencia", 
    ylab = "Densidad")

par(op)

¿Qué tan apropiado es este modelo para nuestro datos? Si miramos los datos, parece haber cierta autocorrelación en la presencia de la especie (los \(0\) y \(1\) están amontonados), pero nuestro modelo asume independencia entre parcelas. ¿Es posible que la estructura espacial que vemos sea producto del azar? En un contexto Bayesiano podemos responder a esta pregunta teniendo en cuenta no solo propiedades de los datos y del modelo de datos (distribución Binomial) sino también la incertidumbre alrededor de las estimaciones de los parámetros (en este caso tenemos un solo parámetro, la probabilidad de presencia).

Cuando hacemos un posterior predictive check, queremos ver si nuestros datos son consistentes con el modelo que ajustamos a ellos, teniendo en cuenta la incertidumbre existente acerca de qué valores tienen los parámetros. La idea es bastante intuitiva: si el modelo es apropiado para nuestros datos, entonces nuestros datos deberían ser coherentes con datos generados (simulados) a partir de ese modelo. Hay varias formas de definir si los datos observados son consistentes o son coherentes con el modelo que ajustamos. Una forma es definir como “test quantity” una propiedad de los datos que queremos que el modelo capture. Por ejemplo, para los datos que estamos considerando podemos calcular cuántas veces hay un cambio de presencia a ausencia en la transecta. En los datos originales, encontramos apenas tres cambios y queremos ver qué tan probable es encontrar ese valor si los datos provienen de una distribución Binomial con probabilidad de éxito estimada a partir de los datos. Para eso podemos ver si la cantidad de cambios observada cae dentro de la distribución predictiva posterior de este estadístico.

¿Pero, qué es la distribución predictiva posterior?

Empecemos desde el principio… La regla de Bayes dice que

\[ p(\theta|y) = \frac{p(y| \theta) p(\theta)}{p(y)} \]

Donde \(p(y)\) es la probabilidad total (o marginal) de los datos, que es algo que podemos calcular así:

\[ p(y) = \int p(y|theta) p(\theta) d \theta \]

Esa es la distribución predictiva previa es decir, una distribución que predice la probabilidad de datos que no hemos observado en base al modelo de datos (likelihood) y las previas. Una vez que observamos un set de datos (\(y_o\)), podemos calcular la distribución predictiva posterior para datos que todavía no observamos pero en base a lo que aprendimos con los datos \(y_o\):

\[ p(y^{\text{rep}}|y_o) = \int p(y^{\text{rep}}|\theta) p(\theta| y_o) d \theta \]

donde \(y^{\text{rep}}\) es una nueva observación o un set de datos nuevos. Ahora podemos ver qué propiedades tendrían estos \(y^{\text{rep}}\) y ver si los datos observados son consistentes con estas predicciones. Si usamos un “test quantity” \(T(y, \theta)\), podemos calcular una especie de valor de p (Bayesian p-value, aunque suene perverso) para lo cual tenemos que resolver esta cuenta:

\[ p_B \left( T(y^{\text{rep}}, \theta ) \geq T(y, \theta) |y \right) \]

Para lo cual hay que calcular:

\[ p_B = \int \int I_{T(y^{\text{rep}, \theta}) \geq T(y, \theta)} p(y^{\text{rep}}| \theta) p(\theta|y) d y^{\text{rep}} d \theta \]

donde \(I\) es una función indicadora que vale \(1\) si la condición se cumple y \(0\) si no se cumple. Si bien esta equación parece formidable, veremos que es fácil resolverla usando Monte Carlo. Lo que tenemos que hacer es obtener valores de la posterior y a partir de ellos simular valores de \(y^{\text{rep}}\), calcular \(T\) para cada una de estas réplicas y luego ver qué proporción de estos son más extremos o iguales que los observados.

Para hacer estas cuentas con el ejemplo de presencia/ausencia de una especie en una transecta, primero vamos a escribir una función que cuente cuántos cambios hay en una secuencia:

countsw <- function(y) {
    count <- 0
    for (i in 2:length(y)) {
        if (y[i] != y[i - 1]) 
            count <- count + 1
    }
    return(count)
}

Con esta función podemos ver cuál es el valor observado de número de cambios y luego hacer unas cuantas simulaciones de datos a partir de la posterior de la probabilidad de presencia.

T_o <- countsw(y)  # número de cambios observado

nsims <- 10000
theta <- rbeta(nsims, 1 + sum(y), 1 + length(y) - sum(y))
y_rep <- matrix(NA, nrow = nsims, ncol = length(y))

for (i in 1:nsims) {
    y_rep[i, ] <- rbinom(length(y), 1, theta[i])
}

Ts <- numeric(nsims)
for (i in 1:nsims) {
    Ts[i] <- countsw(y_rep[i, ])  # núm cambios
}

Ahora podemos ver cómo es la distribución del número de cambios para los datos simulados y ver dónde cae la cantidad de cambios observada en el set de datos original

op <- par(cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, las = 1, bty = "n")
plot(table(Ts)/length(Ts), col = "gray", lwd = 4, ylab = "Probabilidad")
arrows(T_o, 0.1, T_o, 0.025, angle = 20, lwd = 3)

par(op)

Y podemos calcular un valor de p bayesiano \O_o//

p_B <- length(which(Ts <= T_o))/nsims
p_B

## [1] 0.0283

¿Qué quiere decir este valor respecto a la pregunta original de si el modelo que usamos es apropiado para estos datos?

Más allá de este ejemplo de cómo hacer un posterior predictive check, es importante que vean con qué facilidad podemos simular datos del modelo teniendo en cuenta lo que aprendimos de los datos (las posteriores). Con los datos que simulamos a partir de las posteriores podemos hacer cualquier cosa que nos resulte interesante…

Vamos a usar los modelos que ajustamos en el Práctico de regresiones.

Polinización de Pomelos

Primero cargamos los datos y corremos el modelo otra vez:

# Opción 1:
apis <- read.table("https://sites.google.com/site/modelosydatos/polen_Apis.csv", 
    header = T, sep = ",")

# Opción 2: apis <- read.table('polen_Apis.csv', header = T, sep = ',')

Asumamos que la relación entre estas variables es lineal:

cat(file = "ml.bug", "
      model  {
           # Función de likelihood
           for( i in 1:n){
               granos[i] ~ dpois(lambda[i])
               lambda[i] <- b0 + b1 * visitas[i]
           }
           
           # Previas
           b0 ~ dnorm(0,0.01)
           b1 ~ dnorm(0,0.01)
    }")

Ahora vamos a definir qué cosas le vamos a pasar a JAGS, qué cosas queremos que guarde, y creamos una función para generar valores iniciales de las cadenas Markovianas. También definimos cuántas iteraciones queremos correr por cadena, cuántos valores vamos a ir descartado en una secuencia (“thin”), qué largo tiene el “burn in” y cuántas cadenas queremos correr:

# Defino las variables:
granos <- apis$granos
visitas <- apis$visitas
n <- length(granos)  # número de observaciones

# Armo una lista con los nombres de los datos
data <- list("granos", "visitas", "n")

# Armo un vector con los nombres de los parámetros
params <- c("b0", "b1")

# función para generar valores iniciales
inits <- function() list(b0 = rnorm(1, 20, 1), b1 = rnorm(1, 10, 0.5))

ni <- 5000
nt <- 1
nb <- 2500
nc <- 3

Ahora llamamos a JAGS usando la función jags del paquete jagsUI:

library(jagsUI)
ml.sim <- jags(data = data, inits = inits, parameters.to.save = params, model.file = "ml.bug", 
    n.chains = nc, n.iter = ni, n.burnin = nb, n.thin = nt)

## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 190
##    Unobserved stochastic nodes: 2
##    Total graph size: 407
## 
## Initializing model
## 
## Adaptive phase..... 
## Adaptive phase complete 
##  
## 
##  Burn-in phase, 2500 iterations x 3 chains 
##  
## 
## Sampling from joint posterior, 2500 iterations x 3 chains 
##  
## 
## Calculating statistics....... 
## 
## Done.

Como siempre, chequeamos si convergieron las cadenas y si el n.eff es suficiente.

print(ml.sim)

## JAGS output for model 'ml.bug', generated by jagsUI.
## Estimates based on 3 chains of 5000 iterations,
## adaptation = 100 iterations (sufficient),
## burn-in = 2500 iterations and thin rate = 1,
## yielding 7500 total samples from the joint posterior. 
## MCMC ran for 0.029 minutes at time 2019-03-06 20:03:23.
## 
##              mean    sd     2.5%      50%    97.5% overlap0 f  Rhat n.eff
## b0         14.860 0.362   14.159   14.855   15.588    FALSE 1 1.000  3933
## b1          6.858 0.264    6.340    6.858    7.374    FALSE 1 1.002  1062
## deviance 5418.518 1.984 5416.558 5417.915 5423.903    FALSE 1 1.000  7500
## 
## Successful convergence based on Rhat values (all < 1.1). 
## Rhat is the potential scale reduction factor (at convergence, Rhat=1). 
## For each parameter, n.eff is a crude measure of effective sample size. 
## 
## overlap0 checks if 0 falls in the parameter's 95% credible interval.
## f is the proportion of the posterior with the same sign as the mean;
## i.e., our confidence that the parameter is positive or negative.
## 
## DIC info: (pD = var(deviance)/2) 
## pD = 2 and DIC = 5420.486 
## DIC is an estimate of expected predictive error (lower is better).

Nuevamente chequeamos si convergieron las cadenas y si el n.eff es aceptable. Como parece que está todo en orden vamos a hacer un posterior predictive check. La idea atrás de estos chequeos es ver si datos simulados con este modelo se parecen a los datos que observamos. Vamos a simular datos a partir de nuestro modelo teniendo en cuenta la incertidumbre en los valores de los parámetros. Para eso vamos a usar nuestra función de likelihood (modelo de datos) y las posteriores de los parámetros que gobiernan la relación entre el número de visitas y los granos de polen depostiados.

El “modelo de datos” es:

\[ \begin{aligned} \text{granos}_{i} \sim \text{Pois} (\lambda_i) \\ \lambda_i = \beta_0 + \beta_1 \text{visitas}_i \\ \end{aligned} \]

Para simular datos a partir de este modelo vamos a usar la función rpois de R y para calcular \(\lambda\) vamos a tomar muestras de la posterior de b0 y b1.

nsims <- 1000  # vamos a simular 1000 sets de 'datos''
polen_rep <- matrix(NA, nrow = nsims, ncol = n)  # armamos una matriz donde guardamos los datos simulados
nsample <- length(ml.sim$sims.list$b0)  # número de muestras de la posterior deiponibles 

for (i in 1:nsims) {
    idx <- sample.int(nsample, 1)  # número de muestra que vamos a usar
    polen_rep[i, ] <- rpois(n, lambda = ml.sim$sims.list$b0[idx] + ml.sim$sims.list$b1[idx] * 
        visitas)  # simulamos los datos con la función rpois
}

Ahora podemos calcular por ejemplo el promedio de polen depositado de los datos simulados y ver si se parecen a la de los datos.

op <- par(cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, las = 1, bty = "n")
plot(density(rowMeans(polen_rep)), lwd = 2, main = "", xlab = "Promedio de granos de polen")
abline(v = mean(granos), lty = 2, lwd = 3)

par(op)

Les parece que el modelo logra representar bien el promedio de granos de polen depositado?

Veamos ahora qué pasa con la variabilidad alrededor de la media. Por ejemplo, usando el desvío

sd_s <- numeric(nsims)  # generamos un vector de desvios
sd_o <- sd(granos)  # calcuamos el desvio de los datos

for (i in 1:nsims) sd_s[i] <- sd(polen_rep[i, ])  # calculamos el desvío de los datos simulados 

op <- par(cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, las = 1, bty = "n")
plot(density(sd_s), lwd = 2, main = "", xlab = "Granos de polen (sd)", xlim = c(10, 
    25))
abline(v = mean(granos), lty = 2, lwd = 3)
par(op)
abline(v = sd_o, lwd = 2)

Parece que este modelo no es muy apropiado!

Otra cosa que podemos probar es una comparación gráfica a ver si hay partes de la relación entre visitas y polen depositado que ajusta mejor que otra. Por ejemplo, podemos hacer un gráfico con los datos observados y ver cómo quedan respecto a la posterior predictiva. Una buena forma de ver si los datos son posibles dentro de la posterior predictiva es estimar los intervalos de credibilidad de los datos simulados y luego ver si las observaciones caen dentro de estos intervalos. Para calcular los intervalos de credibilidad vamos a usar la función HPDinterval del paquete coda:

library(coda)
tmp <- HPDinterval(as.mcmc(polen_rep))
idx <- sort(visitas, index.return = TRUE)

op <- par(cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, las = 1, bty = "n")
plot(visitas, granos, pch = 21, bg = "gray", cex = 1.5)
lines(visitas[idx$ix], tmp[idx$ix, 1], lwd = 3, col = "darkgray")
lines(visitas[idx$ix], tmp[idx$ix, 2], lwd = 3, col = "darkgray")

par(op)

¿Qué tal el ajuste?

Ejercicio:

Realicen estas mismas comparaciones entre datos simulados y datos observados para el modelo alternativo que propusieron cuando hicimos el práctico de regresión (o con otro que se les ocurra ahora).
¿Qué modelo parece más apropiado?
¿Cómo comparan los valores de DIC?

Ejemplo de frambuesas

Ahora vamos a hacer un posterior predictive check para el ejemplo de las frambuesas.

datos <- read.table("https://sites.google.com/site/modelosydatos/Datosframbuesas.txt", 
    header = TRUE)

Vamos a probar con el modelo de “partial pooling” en la ordenada al origen de la regresión logística para la probabilidad de daño de pistilos:

cat(file = "mbin.bug", "
model {
        for(i in 1: n) {
           rotos[i] ~ dbin(p[i], N[i])
           logit(p[i]) <- b0[Nsitio[i]] + b1 * freq[i]
        }
    
       for( j in 1: S){
          b0[j] ~ dnorm(mub0, taub0)
    }
    
    b1 ~ dnorm(0, 0.01)
    mub0  ~ dnorm(0, 0.01)
    taub0 <- 1/(sdb0 * sdb0)
    sdb0 ~ dunif(0, 10)
    }
    ")

Luego definimos los datos, la función para los valores iniciales y los parámetros que queremos guardar.

rotos <- datos$Rotos
freq <- datos$Frecuencia
n <- length(rotos)
N <- datos$Total
Nsitio <- as.numeric(datos$Sitio)
S <- max(Nsitio)

data <- list("rotos", "freq", "N", "n", "S", "Nsitio")
params <- c("b0", "b1", "mub0", "sdb0")
inits <- function() {
    list(b0 = rnorm(S, 0, 1), b1 = rnorm(1, 0, 1), mub0 = rnorm(1, 0, 1), sdb0 = runif(1, 
        0, 10))
}

Ahora llamamos a JAGS

library(jagsUI)
mb.sim <- jags(data, inits, params, model.file = "mbin.bug", n.chains = 3, n.iter = 10000, 
    n.burnin = 5000, n.thin = 3)

## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 477
##    Unobserved stochastic nodes: 19
##    Total graph size: 1986
## 
## Initializing model
## 
## Adaptive phase..... 
## Adaptive phase complete 
##  
## 
##  Burn-in phase, 5000 iterations x 3 chains 
##  
## 
## Sampling from joint posterior, 5000 iterations x 3 chains 
##  
## 
## Calculating statistics....... 
## 
## Done.

Vemos que hayan convergido y que los tamaños efectivos de las muestras de las posteriores (n.eff) sean razonables.

print(mb.sim)

## JAGS output for model 'mbin.bug', generated by jagsUI.
## Estimates based on 3 chains of 10000 iterations,
## adaptation = 100 iterations (sufficient),
## burn-in = 5000 iterations and thin rate = 3,
## yielding 4998 total samples from the joint posterior. 
## MCMC ran for 0.22 minutes at time 2019-03-06 20:03:27.
## 
##              mean    sd     2.5%      50%    97.5% overlap0     f  Rhat
## b0[1]      -3.860 0.629   -5.158   -3.856   -2.640    FALSE 1.000 1.001
## b0[2]      -2.936 0.344   -3.619   -2.939   -2.266    FALSE 1.000 1.001
## b0[3]      -4.197 0.531   -5.380   -4.139   -3.276    FALSE 1.000 1.023
## b0[4]      -3.961 0.507   -4.992   -3.941   -2.997    FALSE 1.000 1.001
## b0[5]      -4.276 0.761   -5.787   -4.265   -2.806    FALSE 1.000 1.000
## b0[6]      -4.578 0.451   -5.516   -4.562   -3.750    FALSE 1.000 1.011
## b0[7]      -4.793 0.462   -5.751   -4.776   -3.921    FALSE 1.000 1.000
## b0[8]      -3.075 0.284   -3.635   -3.077   -2.505    FALSE 1.000 1.012
## b0[9]      -2.900 0.347   -3.645   -2.871   -2.265    FALSE 1.000 1.033
## b0[10]     -4.591 0.445   -5.525   -4.570   -3.773    FALSE 1.000 1.005
## b0[11]     -1.860 0.870   -3.529   -1.880   -0.142    FALSE 0.982 1.009
## b0[12]     -5.556 0.804   -7.180   -5.550   -3.983    FALSE 1.000 1.000
## b0[13]     -2.732 0.325   -3.404   -2.710   -2.154    FALSE 1.000 1.012
## b0[14]     -4.324 0.572   -5.450   -4.317   -3.217    FALSE 1.000 1.040
## b0[15]     -2.869 0.811   -4.492   -2.868   -1.265    FALSE 1.000 1.001
## b0[16]     -3.423 0.888   -5.187   -3.415   -1.701    FALSE 1.000 1.001
## b1          1.606 0.281    1.055    1.606    2.167    FALSE 1.000 1.001
## mub0       -3.741 0.519   -4.790   -3.729   -2.742    FALSE 1.000 1.000
## sdb0        1.161 0.282    0.747    1.118    1.834    FALSE 1.000 1.007
## deviance 1276.303 6.099 1266.382 1275.560 1289.787    FALSE 1.000 1.015
##          n.eff
## b0[1]     4998
## b0[2]     4441
## b0[3]      126
## b0[4]     3327
## b0[5]     4998
## b0[6]      241
## b0[7]     2859
## b0[8]      203
## b0[9]       69
## b0[10]     604
## b0[11]     258
## b0[12]    4998
## b0[13]     528
## b0[14]      81
## b0[15]    3023
## b0[16]    4998
## b1        3489
## mub0      4982
## sdb0       308
## deviance   205
## 
## Successful convergence based on Rhat values (all < 1.1). 
## Rhat is the potential scale reduction factor (at convergence, Rhat=1). 
## For each parameter, n.eff is a crude measure of effective sample size. 
## 
## overlap0 checks if 0 falls in the parameter's 95% credible interval.
## f is the proportion of the posterior with the same sign as the mean;
## i.e., our confidence that the parameter is positive or negative.
## 
## DIC info: (pD = var(deviance)/2) 
## pD = 18.4 and DIC = 1294.73 
## DIC is an estimate of expected predictive error (lower is better).

Vamos a empezar simulando datos como hicimos con el modelo anterior. Acá hay que recordar que tenemos distintos sitios y una ordenada al origen para cada uno de ellos. Entonces podemos simular un set de datos para cada sitio y comparar con los datos de ese sitio en particular. Lo que vamos a hacer es simular \(1000\) sets de datos para cada sitio y calcular la media y desvío de cada uno, y para esas \(1000\) medias y \(1000\) desvíos vamos a calcular un intervalo de credibilidad.

nsims <- 1000  # vamos a simular 1000 sets de datos para cada sitio
nsample <- length(mb.sim$sims.list$b1)
datamean <- numeric(S)  # vector para las medias de sitio (datos reales)
datasd <- numeric(S)  # vector para desvio de sitio (datos reales)
low <- numeric(S)  # vector para el limite inferior del intervalo 
# de credibilidad para las medias (datos simulados)
up <- numeric(S)  # vector que guarda el limite superior del intervalo de
# credibilidad para las medias (datos simulados)
slow <- numeric(S)  # idem desvios
sup <- numeric(S)

for (j in 1:S) {
    xfreq <- freq[which(Nsitio == j)]
    Rpistilos <- rotos[which(Nsitio == j)]
    Tpistilos <- N[which(Nsitio == j)]
    rep <- matrix(NA, nrow = nsims, ncol = length(xfreq))
    
    for (i in 1:nsims) {
        idx <- sample.int(nsample, 1)
        rep[i, ] <- rbinom(length(xfreq), size = Tpistilos, prob = plogis(mb.sim$sims.list$b0[idx, 
            j] + mb.sim$sims.list$b1[idx] * xfreq))
    }
    
    sim_mean <- numeric(nsims)
    for (s in 1:nsims) sim_mean[s] <- mean(rep[s, ])
    sim_sd <- numeric(nsims)
    for (s in 1:nsims) sim_sd[s] <- sd(rep[s, ])
    datamean[j] <- mean(Rpistilos)
    datasd[j] <- sd(Rpistilos)
    low[j] <- HPDinterval(as.mcmc(sim_mean), prob = 0.95)[1]
    up[j] <- HPDinterval(as.mcmc(sim_mean), prob = 0.95)[2]
    slow[j] <- HPDinterval(as.mcmc(sim_sd), prob = 0.95)[1]
    sup[j] <- HPDinterval(as.mcmc(sim_sd), prob = 0.95)[2]
}

length(which(datamean > low & datamean < up))/S

## [1] 1

length(which(datasd > slow & datasd < sup))/S

## [1] 0.1875

op <- par(cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, las = 1, bty = "n")
plot(1:16, datamean, ylim = c(0, 6), xlim = c(0, 20), pch = 16, xlab = "Sitios", 
    ylab = "Promedio de pistilos rotos")
arrows(1:16, low, 1:16, up, code = 0, col = "gray30")

par(op)

op <- par(cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, las = 1, bty = "n")
plot(1:16, datasd, ylim = c(0, 4), xlim = c(0, 20), pch = 16, xlab = "Sitios", 
    ylab = "sd pistilos rotos")
arrows(1:16, slow, 1:16, sup, code = 0, col = "gray30")

¿Cómo evalúan estos ajustes?

Ejercicios

Armar un posterior predictive check para el modelo que considera variabilidad entre plantas.
Lo mismo para los modelos de Aliso

Juan Manuel Morales. 6 de Septiembre del 2015. Última actualización: 2019-03-06