Prior predictive distributions
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. Si prestamos atención a esta distribución o a ciertas propiedades de ella, podemos evaluar si las previas que elegimos tienen sentido, aún sin haber colectado ningún dato.
Por ejemplo, si que queremos modelar la distribución de pesos de monito del monte adultos, podemos asumir que tienen una distribución normal, caracterizada por un valor de media y cierto desvío estándar. Una opción de previas no-informativas sería, por ejemplo \(\mu \sim N(0, 100)\) y \(\sigma \sim U(0,100)\):
m = 0
s = 100
op = par(mfrow = c(1,2), bty = "l", las = 1)
curve(dnorm(x, m, s), xlim = c(-400, 400),
xlab = expression(mu), ylab = "densidad")
curve(dunif(x, 0, 100), xlim = c(-1, 101),
xlab = expression(sigma), ylab = "densidad")
par(op)¿Qué valores de pesos obtendríamos con estas previas? Para responder a eso, podemos simular “datos” a partir de ellas y del modelo de datos.
mus = rnorm(1e4, m, s)
sds = runif(1e4, 0, 100)
psims = rnorm(1e4, mus, sds)
plot(density(psims), main = "", xlab = "peso simulado")
Claramente, estos valores son absurdos. El peso de los monitos adultos es variable pero en mi experiencia de trabajo con esta especie en Llao Llao, puedo decir que nunca vi un monito adulto de menos de 16 gramos o de más de 55 gramos. Entonces, para poner una previa más razonable, podríamos usar una normal con \(\mu = 30\) y \(\sigma = 5\). Para el desvío estándar podemos usar como previa una uniforme entre \(0\) y \(5\).
m = 30
s = 5
op = par(mfrow = c(1,2), bty = "l", las = 1)
curve(dnorm(x, m, s), xlim = c(10, 55))
abline(v = m - 3*s, lty = 3)
abline(v = m + 3*s, lty = 3)
curve(dunif(x, 0, 5), xlim = c(-1, 6))
par(op)Veamos cómo queda nuestra prior predictive ahora…
mus = rnorm(1e4, m, s)
sds = runif(1e4, 0, 5)
lsims = rnorm(1e4, mus, sds)
plot(density(lsims), main = "", xlab = "peso simulado")
Vemos que estos valores son razonables para estos datos y que permiten valores extremos pero quizás no imposibles.
Consideremos ahora el caso de una regresión para el peso de los animales en función de su largo. Asumiendo una relación linear estándar, nuestro modelo de datos es:
\[ \text{peso} \sim N \left(\alpha + \beta \times \text{largo} , \sigma \right) \] En general, nos conviene trabajar con variables predictoras centradas, de manera que la ordenada al origen (\(\alpha\)) represente el valor de la variable de respuesta para valores promedios de la predictora. Reemplazando a “peso” por \(y\) y a “largo” por \(x\) de manera de generalizar la notación, tenemos:
\[ y \sim N \left(\alpha + \beta \left(x - \bar{x} \right), \sigma \right) \] Ahora tenemos que definir previas para \(\alpha, \beta\) y \(\sigma\). Veamos qué pasa con la regresión si usamos \(\alpha \sim N(30, 5)\) y una previa poco informativa para \(\beta\): \(\beta \sim N(0, 100)\)
Carguemos los datos y veamos qué dicen estas previas
datos = read.table("https://github.com/jmmorales/cursoMyD/raw/master/Data/peso.txt",
header = TRUE)
a = rnorm(100, m, s)
b = rnorm(100, 0, 100)
plot( NULL , xlim = range(datos$largo) ,
ylim = c(-30, 100) ,
xlab = "largo",
ylab = "peso" )
abline( h = 16, lty = 2, lwd = 2, col = rgb(0,0,0, 0.5) )
abline( h = 55, lty = 2, lwd = 2, col = rgb(0,0,0, 0.5) )
xbar <- mean(datos$largo)
for ( i in 1:100 ){
curve( a[i] + b[i] * (x - xbar),
from = min(datos$largo),
to=max(datos$largo),
add=TRUE,
col=rgb(0,0,0,0.2) )
}
Con estas previas vemos que las líneas de regresión son muy poco realistas.
Probemos entonces con \(\beta \sim N(0, 3)\)
a = rnorm(100, m, s)
b = rnorm(100, 0, 3)
plot( NULL , xlim = range(datos$largo) ,
ylim = c(-30, 100) ,
xlab = "largo",
ylab = "peso" )
abline( h = 16, lty = 2, lwd = 2, col = rgb(0,0,0, 0.5) )
abline( h = 55, lty = 2, lwd = 2, col = rgb(0,0,0, 0.5) )
xbar <- mean(datos$largo)
for ( i in 1:100 ){
curve( a[i] + b[i] * (x - xbar),
from = min(datos$largo),
to=max(datos$largo),
add=TRUE,
col=rgb(0,0,0,0.2) )
}
Ahora ajustemos el modelo a los datos. Para eso escribimos el modelo en lenguaje de BUGS
cat(file = "ml.bug",
"
model {
# Función de likelihood
for( i in 1:n){
y[i] ~ dnorm (mu[i], tau)
mu[i] <- a + b * x[i]
}
# Previas
a ~ dnorm(21, 0.25)
b ~ dnorm(0, tau_b)
tau <- 1/(sigma*sigma)
sigma ~ dunif(0, 5)
}
")Ahora cargamos los datos y usamos el paquete jagsUI para correr el análisis. Vamos a comparar los resultados con las alternativas de previas para \(\beta\)
library(jagsUI)
# defino las variables
x <- datos$largo - mean(datos$largo)
y <- datos$peso
# vector con los parámetros para que JAGS guarde los valores de las cadenas Markovianas
params <- c("a", "b", "sigma")
# función para generar valores iniciales de los parámetros (no confundir con las previas!)
inits <- function() list(a = runif(1, 16, 27),
b = runif(1, 0, 1),
sigma = runif(1, 0, 5))Ahora llamamos a JAGS para que genere las cadenas Markovianas para las dos opciones de previas
dat <- list(x = x,
y = y,
n = length(datos$largo),
tau_b = (1/100)^2
)
m100.sim <- jags(data = dat,
inits = inits,
parameters.to.save = params,
model.file = "ml.bug",
n.chains = 3,
n.iter = 1000)##
## Processing function input.......
##
## Done.
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 139
## Unobserved stochastic nodes: 3
## Total graph size: 329
##
## Initializing model
##
## Adaptive phase.....
## Adaptive phase complete
##
## No burn-in specified
##
## Sampling from joint posterior, 1000 iterations x 3 chains
##
##
## Calculating statistics.......
##
## Done.dat <- list(x = x,
y = y,
n = length(datos$largo),
tau_b = (1/3)^2
)
m3.sim <- jags(data = dat,
inits = inits,
parameters.to.save = params,
model.file = "ml.bug",
n.chains = 3,
n.iter = 1000)##
## Processing function input.......
##
## Done.
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 139
## Unobserved stochastic nodes: 3
## Total graph size: 329
##
## Initializing model
##
## Adaptive phase.....
## Adaptive phase complete
##
## No burn-in specified
##
## Sampling from joint posterior, 1000 iterations x 3 chains
##
##
## Calculating statistics.......
##
## Done.y vemos los resultados
print(m100.sim)## JAGS output for model 'ml.bug', generated by jagsUI.
## Estimates based on 3 chains of 1000 iterations,
## adaptation = 100 iterations (sufficient),
## burn-in = 0 iterations and thin rate = 1,
## yielding 3000 total samples from the joint posterior.
## MCMC ran for 0.005 minutes at time 2021-12-07 23:06:32.
##
## mean sd 2.5% 50% 97.5% overlap0 f Rhat n.eff
## a 30.165 0.417 29.338 30.167 30.954 FALSE 1 1.003 852
## b 4.244 0.486 3.320 4.239 5.218 FALSE 1 1.001 2739
## sigma 4.967 0.033 4.883 4.976 4.999 FALSE 1 1.001 3000
## deviance 968.845 3.218 964.819 968.110 976.939 FALSE 1 1.001 2052
##
## 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 = 5.2 and DIC = 974.023
## DIC is an estimate of expected predictive error (lower is better).print(m3.sim)## JAGS output for model 'ml.bug', generated by jagsUI.
## Estimates based on 3 chains of 1000 iterations,
## adaptation = 100 iterations (sufficient),
## burn-in = 0 iterations and thin rate = 1,
## yielding 3000 total samples from the joint posterior.
## MCMC ran for 0.008 minutes at time 2021-12-07 23:06:32.
##
## mean sd 2.5% 50% 97.5% overlap0 f Rhat n.eff
## a 30.170 0.407 29.387 30.157 30.979 FALSE 1 1.000 3000
## b 4.115 0.472 3.177 4.117 5.022 FALSE 1 1.001 1445
## sigma 4.964 0.034 4.875 4.974 4.999 FALSE 1 1.008 352
## deviance 968.938 3.177 964.875 968.258 976.602 FALSE 1 1.002 1035
##
## 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 = 5 and DIC = 973.978
## DIC is an estimate of expected predictive error (lower is better).En este caso no parece haber mucha diferencia entre las dos previas. Esto se debe a que tenemos bastantes datos y pocas variables predictoras
plot(density(m100.sim$sims.list$b), main = "",
xlab = expression(beta), bty = "l", ylim = c(0,1))
lines(density(m3.sim$sims.list$b), col=2)
Veamos qué pasa con una muestra más chica:
n = 10
set.seed(1234)
idx = sample(1:length(datos$largo), n)
dat <- list(x = x[idx],
y = y[idx],
n = n,
tau_b = (1/100)^2
)
m100r.sim <- jags(data = dat,
inits = inits,
parameters.to.save = params,
model.file = "ml.bug",
n.chains = 3,
n.iter = 10000)##
## Processing function input.......
##
## Done.
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 10
## Unobserved stochastic nodes: 3
## Total graph size: 47
##
## Initializing model
##
## Adaptive phase.....
## Adaptive phase complete
##
## No burn-in specified
##
## Sampling from joint posterior, 10000 iterations x 3 chains
##
##
## Calculating statistics.......
##
## Done.dat <- list(x = x[idx],
y = y[idx],
n = n,
tau_b = (1/3)^2
)
m3r.sim <- jags(data = dat,
inits = inits,
parameters.to.save = params,
model.file = "ml.bug",
n.chains = 3,
n.iter = 10000)##
## Processing function input.......
##
## Done.
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 10
## Unobserved stochastic nodes: 3
## Total graph size: 47
##
## Initializing model
##
## Adaptive phase.....
## Adaptive phase complete
##
## No burn-in specified
##
## Sampling from joint posterior, 10000 iterations x 3 chains
##
##
## Calculating statistics.......
##
## Done.plot(density(m100r.sim$sims.list$b), main = "",
xlab = expression(beta), bty = "l", ylim = c(0,0.4))
lines(density(m3r.sim$sims.list$b), col=2)
print(m100r.sim)## JAGS output for model 'ml.bug', generated by jagsUI.
## Estimates based on 3 chains of 10000 iterations,
## adaptation = 100 iterations (sufficient),
## burn-in = 0 iterations and thin rate = 1,
## yielding 30000 total samples from the joint posterior.
## MCMC ran for 0.016 minutes at time 2021-12-07 23:06:33.
##
## mean sd 2.5% 50% 97.5% overlap0 f Rhat n.eff
## a 28.408 1.236 26.015 28.403 30.819 FALSE 1 1.000 30000
## b 6.208 1.617 3.036 6.198 9.410 FALSE 1 1.000 15119
## sigma 4.820 0.166 4.388 4.867 4.995 FALSE 1 1.002 1836
## deviance 84.582 4.923 76.787 83.967 95.750 FALSE 1 1.000 23088
##
## 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 = 12.1 and DIC = 96.702
## DIC is an estimate of expected predictive error (lower is better).print(m3r.sim)## JAGS output for model 'ml.bug', generated by jagsUI.
## Estimates based on 3 chains of 10000 iterations,
## adaptation = 100 iterations (sufficient),
## burn-in = 0 iterations and thin rate = 1,
## yielding 30000 total samples from the joint posterior.
## MCMC ran for 0.008 minutes at time 2021-12-07 23:06:34.
##
## mean sd 2.5% 50% 97.5% overlap0 f Rhat n.eff
## a 28.552 1.223 26.143 28.557 30.949 FALSE 1 1.000 30000
## b 4.814 1.432 2.017 4.819 7.591 FALSE 1 1.000 30000
## sigma 4.821 0.165 4.390 4.868 4.995 FALSE 1 1.001 4649
## deviance 84.549 4.949 76.827 83.923 95.918 FALSE 1 1.000 30000
##
## 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 = 12.2 and DIC = 96.795
## DIC is an estimate of expected predictive error (lower is better).Como dice Richard McElreath, las previas no-informativas permiten que el modelo “se entusiasme” con los datos que observa de manera que es probable que estemos ajustando “ruido” y tengamos poca capacidad predictiva. Por eso es que previas más informativas y que producen predicciones previas razonables a veces se las llama “regularizadoras”.
Posterior predictive checks
La distribución posterior predictiva es una distribución de datos que podríamos obtener considerando que el modelo de datos es el correcto y teniendo en cuenta la incertidumbre en los parámetros una vez que ya vimos un set de datos. Es decir \[ p(y^{\text{rep}}|y) = \int p(y^{\text{rep}}|\theta) p(\theta| y) d \theta \]
Para aproximar esa distribución, usamos muestras de la posterior conjunta y el modelo de datos
nsims <- 1000
n = length(y)
# matriz donde guardamos los datos simulados
y_rep <- matrix(NA, nrow = nsims, ncol = n)
nsample <- length(m3.sim$sims.list$a) # muestras de la posterior disponibles
# vectores para guardar las sumas de residuos
sum_res = numeric(nsims)
sum_res_rep = numeric(nsims)
for(i in 1:nsims){
idx = sample.int(nsample, 1) # número de muestra que vamos a usar
m_rep = m3.sim$sims.list$a[idx] + m3.sim$sims.list$b[idx] * x
sd_rep = m3.sim$sims.list$sigma[idx]
# simulamos los datos
y_rep[i, ] = rnorm(n, m_rep, sd_rep)
# suma de residuos
sum_res[i] = sum(y - m_rep)
sum_res_rep[i] = sum(y_rep[i,] - m_rep)
}Podemos hacer comparaciones visuales entre los datos observados y los simulados
op = par(mfrow = c(5,5), mar = c(2,1,1,2))
hist(datos$peso, main = "", xlab = "", ylab = "" , xlim = c(5,65), ylim = c(0,25), col = 2, 20)
for(i in 1: 24){
hist(y_rep[i,],
main = "",
xlab = "",
ylab = "",
xlim = c(5,65),
ylim = c(0,25),
20)
} 
par(op)Podemos definir valores de “p” sobre alguna cantidad o “test”. Por ejemplo, veamos el valor de peso máximo:
test = max(datos$peso)
idx = numeric(nsims)
for(i in 1:nsims){
if(max(y_rep[i,]) > test){
idx[i] = 1
}
}
p_value = sum(idx)/nsimsEstas comparaciones también las podemos hacer con gráficos. Por ejemplo, veamos el valor de peso promedio
hist(rowMeans(y_rep), 100, main = "", xlab = "mean y_rep")
abline(v = mean(y), lwd = 3)
y el desvío estándar:
# calculamos el desvío de los datos simulados
sd_sim = numeric(nsims)
for(i in 1:nsims) sd_sim[i] <- sd(y_rep[i,])
hist(sd_sim, 100, main = "", xlab = "sd y_rep", xlim = c(4.5, 8))
abline(v = sd(y), lwd = 3)
Vemos que los datos tienen más variabilidad que la predicha por el modelo. Es decir que hay otras fuentes de variabilidad que no están contempladas por la relación lineal entre largo y peso.
También podemos ver qué pasa con las sumas de residuos que calculamos. Para los datos observados y predichos
plot(sum_res, sum_res_rep, asp = 1)
abline(0, 1)
test = numeric(nsims)
for(i in 1:nsims) if(sum_res[i] > sum_res_rep[i]) test[i] = 1
mean(test)## [1] 0.757Otra forma de ver si los datos observados son compatibles con los simulados es graficar los intervalos de credibilidad del pp-check en función de \(x\):
library(coda)
tmp <- HPDinterval(as.mcmc(y_rep))
idx <- sort(datos$largo, index.return = TRUE)
op <- par(las = 1, bty = "n")
plot(datos$largo, datos$peso, pch = 16, col = gray(0.5, 0.5),
cex = 1.5, ylim = c(10,55))
lines(datos$largo[idx$ix], tmp[idx$ix, 1], lwd = 3, col = "darkgray")
lines(datos$largo[idx$ix], tmp[idx$ix, 2], lwd = 3, col = "darkgray")
par(op)Ejercicios
1- Extender el modelo de regresión para incluir además del largo, al sexo y el ancho de cola como predictoras. Los monitos guardan reservas en la base de la cola, por eso es de esperarse que el ancho de cola esté relacionado al peso. Como largo y ancho de cola tienen rangos de valores diferentes conviene estandarizar estas dos variables antes de incluirlas en la regresión.
2- Una vez que estandarizaron largo y ancho de cola, qué previas para los coeficientes de regresión les parecen razonables?
3- Realizar al menos 3 pp-checks para el nuevo modelo de regresión