¿Cómo obtener predicciones para una variable específica en WinBUGS?

10

Soy un nuevo usuario de WinBUGS y tengo una pregunta por su ayuda. Después de ejecutar el siguiente código, obtuve parámetros de beta0through beta4(estadísticas, densidad), pero no sé cómo obtener la predicción del último valor de h, que configuré NApara modelar en el código.

¿Alguien me puede dar una pista? Cualquier consejo sería muy apreciado.

model {
for(i in 1: N) {
CF01[i] ~ dnorm(0, 20)
CF02[i]  ~ dnorm(0, 1)
h[i] ~ dpois (lambda [i])
log(lambda [i]) <- beta0 + beta1*CF03[i] + beta2*CF02[i] + beta3*CF01[i] + beta4*IND[i]
}
beta0 ~ dnorm(0.0, 1.0E-6)
beta1 ~ dnorm(0.0, 1.0E-6)
beta2 ~ dnorm(0.0, 1.0E-6)
beta3 ~ dnorm(0.0, 1.0E-6)
beta4  <- log(p)
p ~ dunif(lower, upper)
}

INITS
list(beta0 = 0, beta1 = 0, beta2 = 0, beta3 = 0, p = 0.9)

DATA(LIST)
list(N = 154, lower = 0.80, upper = 0.95,

h = c(1, 4, 1, 2, 1, 2, 1, 1, 1, 3, 3, 0, 0, 0, 2, 0, 1, 0, 4, 2,
3, 0, 2, 1, 1, 2, 2, 2, 3, 4, 2, 3, 1, 0, 1, 3, 3, 3, 1, 0, 1,
0, 5, 2, 1, 2, 1, 3, 3, 1, 1, 0, 2, 2, 0, 3, 0, 0, 3, 2, 2, 2,
1, 0, 3, 3, 1, 1, 1, 2, 1, 0, 1, 2, 1, 2, 0, 2, 1, 0, 0, 2, 5,
0, 2, 1, 0, 2, 1, 2, 2, 2, 0, 3, 2, 1, 3, 3, 3, 3, 0, 1, 3, 3,
3, 1, 0, 0, 1, 2, 1, 0, 1, 4, 1, 1, 1, 1, 2, 1, 3, 0, 0, 1, 1,
1, 1, 0, 2, 1, 0, 0, 1, 1, 5, 1, 1, 1, 3, 0, 1, 1, 1, 0, 2, 1,
0, 3, 3, 0, 0, 1, 2, 6, NA),

CF03 = c(-1.575, 0.170, -1.040, -0.010, -0.750,
0.665, -0.250, 0.145, -0.345, -1.915, -1.515,
0.215, -1.040, -0.035, 0.805, -0.860, -1.775,
1.725, -1.345, 1.055, -1.935, -0.160, -0.075,
-1.305, 1.175, 0.130, -1.025, -0.630, 0.065,
-0.665, 0.415, -0.660, -1.145, 0.165, 0.955,
-0.920, 0.250, -0.365, 0.750, 0.045, -2.760,
-0.520, -0.095, 0.700, 0.155, -0.580, -0.970,
-0.685, -0.640, -0.900, -0.250, -1.355, -1.330,
0.440, -1.505, -1.715, -0.330, 1.375, -1.135,
-1.285, 0.605, 0.360, 0.705, 1.380, -2.385, -1.875,
-0.390, 0.770, 1.605, -0.430, -1.120, 1.575, 0.440,
-1.320, -0.540, -1.490, -1.815, -2.395, 0.305,
0.735, -0.790, -1.070, -1.085, -0.540, -0.935,
-0.790, 1.400, 0.310, -1.150, -0.725, -0.150,
-0.640, 2.040, -1.180, -0.235, -0.070, -0.500,
-0.750, -1.450, -0.235, -1.635, -0.460, -1.855,
-0.925, 0.075, 2.900, -0.820, -0.170, -0.355,
-0.170, 0.595, 0.655, 0.070, 0.330, 0.395, 1.165,
0.750, -0.275, -0.700, 0.880, -0.970, 1.155, 0.600,
-0.075, -1.120, 1.480, -1.255, 0.255, 0.725,
-1.230, -0.760, -0.380, -0.015, -1.005, -1.605,
0.435, -0.695, -1.995, 0.315, -0.385, -0.175,
-0.470, -1.215, 0.780, -1.860, -0.035, -2.700,
-1.055, 1.210, 0.600, -0.710, 0.425, 0.155, -0.525,
-0.565),

CF02 = c(NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA, 0.38, 0.06, -0.94,
-0.02, -0.28, -0.78, -0.95, 2.33, 1.43, 1.24, 1.26,
-0.75, -1.5, -2.09, 1.01, -0.05, 2.48, 2.48, 0.46,
0.46, -0.2, -1.11, 0.52, -0.37, 0.58, 0.86, 0.59,
-0.12, -1.33, 1.4, -1.84, -1.4, -0.76, -0.23,
-1.78, -1.43, 1.2, 0.32, 1.87, 0.43, -1.71, -0.54,
-1.25, -1.01, -1.98, 0.52, -1.07, -0.44, -0.24,
-1.31, -2.14, -0.43, 2.47, -0.09, -1.32, -0.3,
-0.99, 1.1, 0.41, 1.01, -0.19, 0.45, -0.07, -1.41,
0.87, 0.68, 1.61, 0.36, -1.06, -0.44, -0.16, 0.72,
-0.69, -0.94, 0.11, 1.25, 0.33, -0.05, 0.87, -0.37,
-0.2, -2.22, 0.26, -0.53, -1.59, 0.04, 0.16, -2.66,
-0.21, -0.92, 0.25, -1.36, -1.62, 0.61, -0.2, 0,
1.14, 0.27, -0.64, 2.29, -0.56, -0.59, 0.44, -0.05,
0.56, 0.71, 0.32, -0.38, 0.01, -1.62, 1.74, 0.27, 0.97,
1.22, -0.21, -0.05, 1.15, 1.49, -0.15, 0.05, -0.87,
-0.3, -0.08, 0.5, 0.84, -1.67, 0.69, 0.47, 0.44,
-1.35, -0.24, -1.5, -1.32, -0.08, 0.76, -0.57,
-0.84, -1.11, 1.94, -0.68),

CF01 = c(NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, -0.117, -0.211, -0.333, -0.229, -0.272,
-0.243, -0.148, 0.191, -0.263, -0.239, -0.168,
-0.381, -0.512, -0.338, -0.296, 0.067, 0.104,
-0.254, -0.167, -0.526, -0.096, -0.43, 0.013,
-0.438, -0.297, -0.131, -0.098, -0.046, -0.063,
-0.194, -0.155, -0.645, -0.603, -0.374, -0.214,
-0.165, -0.509, -0.171, -0.442, -0.468, -0.289,
-0.427, -0.519, -0.454, 0.046, -0.275, -0.401,
-0.542, -0.488, -0.52, -0.018, -0.551, -0.444,
-0.254, -0.286, 0.048, -0.03, -0.015, -0.219,
-0.029, 0.059, 0.007, 0.157, 0.141, -0.035, 0.136,
0.526, 0.113, 0.22, -0.022, -0.173, 0.021, -0.027,
0.261, 0.082, -0.266, -0.284, -0.097, 0.097, -0.06,
0.397, 0.315, 0.302, -0.026, 0.268, -0.111, 0.084,
0.14, -0.073, 0.287, 0.061, 0.035, -0.022, -0.091,
-0.22, -0.021, -0.17, -0.184, 0.121, -0.192,
-0.24, -0.283, -0.003, -0.45, -0.138, -0.143,
0.017, -0.245, 0.003, 0.108, 0.015, -0.219, 0.09,
-0.22, -0.004, -0.178, 0.396, 0.204, 0.342, 0.079,
-0.034, -0.122, -0.24, -0.125, 0.382, 0.072, 0.294,
0.577, 0.4, 0.213, 0.359, 0.074, 0.388, 0.253, 0.167),

IND = c(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, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0))
bayesian bugs prediction winbugs Bo Yu
fuente
¿No solo estás preguntando por el valor de lambda [N]?
whuber
@whuber sí, creo que es correcto, pero más fundamentalmente, debe tener cosas que desea predecir (es decir, tener una distribución posterior para) diferentes de las cosas que ya ha observado. Puede hacer la predicción explícitamente en winbugs o en el postprocesamiento utilizando las muestras de las versiones beta.
atiretoo - reinstalar a monica el
@atiretoo Hasta donde puedo decir, las lambdas son exactamente lo que uno quiere predecir: este es un modelo lineal generalizado para una distribución de Poisson con enlace de registro y las lambdas son los parámetros predichos de Poisson. No han sido observados. Creo que todo lo que hay que hacer aquí es configurar un monitor en lambda [N].
whuber
@whuber, prefiero decir monitor en h[N]lugar de lambda[N]... y obtienes la distribución posterior del valor predicho.
Curioso
@tomek pero h[N]no es el valor predicho: será una colección de sorteos de un conjunto de distribuciones de Poisson pronosticadas. Como tal, combina la variación en los parámetros de Poisson y la variación de esas distribuciones de Poisson. Lo relevante es la distribución posterior de lambda[N].
whuber

Respuestas:

6

Simplemente agregue la variable ha la lista de parámetros a monitorear. Si está utilizando un paquete como R2WinBUGS, agregue la variable ha la lista pasada como parameters.to.saveargumento a la bugsfunción. Luego observe su último valor en h(el que tiene NA): obtendrá una distribución posterior allí.

Esta es la forma habitual de hacer predicciones en inferencia bayesiana ( ver también esta pregunta ). ¡Es agradable y simple! No más separación de evaluación y predicción de parámetros. Todo se hace a la vez. La distribución posterior de los parámetros viene dada por los datos reales y se propaga a los valores de NA (como "predicciones").

Curioso
fuente
Tomás, gracias por tu ayuda. Intento monitorear la variable de h en Sample Monitor Tool pero no funciona. ¿Podrías ayudarme de nuevo? El siguiente es el procedimiento que hice en WinBUGS (no sé cómo usar R2WinBUGS): 1) seleccione Sample en Sample Monitor Tool 2) escriba h en el cuadro blanco marcado nodo 3) haga clic en el botón marcado set 4) h es no en la lista de parámetros que quiero monitorear, mientras que otros parámetros (beta0, beta1, beta2, beta3, p) se muestran en la lista. ¿Sabes cómo puedo agregar "h" a la lista de parámetros que quiero monitorear? ¡Gracias de nuevo!
Bo Yu
@BoYu, no sé cómo hacerlo directamente en WinBUGS ya que ejecuto WinBUGS desde R, usando el paquete R2WinBUGS. Es mucho más práctico porque puedes guardar el script R y ejecutarlo todo como un lote, junto con la producción de tus propios gráficos, etc. Mira aquí, por ejemplo, scripts.
Curioso
Dicho esto, seguramente también será posible en WinBUGS, pero no sé cómo (y supongo que la mayoría de la gente lo llama desde R).
Curioso
1
En primer lugar, ¡gracias, whuber, atiretoo y Tomas! Como Whuber ya mencionó, sí, es un modelo lineal generalizado, la variable de h se ajusta mediante la distribución de Poisson con tasa variable (lambda) condicionada con diferentes predictores (CF01, CF02, CF03 e IND). El último valor de h es lo que necesito saber y no se observa (marcado como NA), mientras que se observan todos los demás valores de h. Creo que Whuber tiene razón, necesito configurar lambda como parámetro en Sample Monitor Tool y verificar las estadísticas del último valor de lambda, y obtener aún más cuál es mi predicción de la última h. Gracias a todos.
Bo Yu
1
@ Tomás, muchas gracias. ¡Sí, tiene usted razón! WinBUGS proporciona la predicción de h [N], incluidas las estadísticas y la densidad de probabilidad. Lo entiendo ahora. Saludos cordiales,
Bo Yu