Como obter previsão para uma variável específica no WinBUGS?

10

Eu sou um novo usuário do WinBUGS e tenho uma pergunta para sua ajuda. Depois de executar o código a seguir, obtive parâmetros de beta0through beta4(estatísticas, densidade), mas não sei como obter a previsão do último valor de h, que defini NAcomo modelo no código.

Alguém pode me dar uma dica? Qualquer conselho seria muito 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
fonte
Você não está apenas pedindo o valor de lambda [N]?
whuber
@ whuber sim, acho que está correto, mas, mais fundamentalmente, você precisa ter as coisas que deseja prever (ou seja, ter uma distribuição posterior para) serem distintas das coisas que você já observou. Você pode fazer a previsão explicitamente no winbugs ou no pós-processamento usando as amostras dos betas.
atiretoo - restabelece monica
@atiretoo Até onde eu sei, as lambdas são exatamente o que se deseja prever: este é um modelo linear generalizado para uma distribuição Poisson com link de log e as lambdas são os parâmetros Poisson previstos. Eles não foram observados. Acredito que tudo o que precisamos fazer aqui é configurar um monitor no lambda [N].
whuber
@ whuber, eu prefiro dizer monitor em h[N]vez de lambda[N]... e você obtém a distribuição posterior do valor previsto.
curioso
@tomek, mas h[N]não é o valor previsto: será uma coleção de empates de um conjunto de distribuições previstas de Poisson. Como tal, combina variação nos parâmetros de Poisson e variação dessas próprias distribuições de Poisson. O que é relevante é a distribuição posterior de lambda[N].
whuber

Respostas:

6

Basta adicionar a variável hà lista dos parâmetros a serem monitorados. Se você estiver usando um pacote como R2WinBUGS, adicione variável hà lista passada ao parameters.to.saveargumento para a bugsfunção Então olhe seu último valor em h(aquele com NA) - você obterá uma distribuição posterior lá.

Esta é a maneira usual de fazer previsões na inferência bayesiana ( veja também esta pergunta ). É simples e agradável! Não há mais separação da avaliação e previsão dos parâmetros. Tudo é feito de uma só vez. A distribuição posterior dos parâmetros é dada pelos dados reais e propagada para os valores de NA (como "previsões").

Curioso
fonte
Tomas, obrigado por sua ajuda. Eu tento monitorar a variável h na Ferramenta de Monitor de Amostra, mas ela não funciona. Você poderia me ajudar de novo? A seguir, é apresentado o procedimento que fiz no WinBUGS (não sei como usar o R2WinBUGS): 1) selecione Amostra na ferramenta Monitor de amostra 2) digite h na caixa branca marcada como nó 3) clique no botão marcado como conjunto 4) h é não está na lista de parâmetros que eu quero monitorar, enquanto outros parâmetros (beta0, beta1, beta2, beta3, p) são mostrados na lista. Você sabe como posso adicionar "h" à lista de parâmetros que quero monitorar? Agradeço mais uma vez!
Bo Yu
@ BoYu, eu não sei como fazê-lo diretamente no WinBUGS desde que eu execute o WinBUGS a partir do R, usando o pacote R2WinBUGS. É muito mais prático, pois você pode salvar o script R e executá-lo em lote, além de produzir seus próprios gráficos etc. Veja aqui exemplos de scripts.
curioso
Dito isto, com certeza também será possível no próprio WinBUGS, mas não sei como (e acho que a maioria das pessoas chama isso de R).
Curioso
11
Em primeiro lugar, obrigado, whuber, atiretoo e Tomas! Como o whuber já mencionou, sim, é um modelo linear generalizado, a variável h é ajustada pela distribuição de Poisson com taxa variável (lambda) condicionada por diferentes preditores (CF01, CF02, CF03 e IND). O último valor de h é o que eu preciso saber e não é observado (marcado como NA), enquanto todos os outros valores de h são observados. Eu acho que whuber está certo, preciso definir lambda como um parâmetro na Ferramenta de Monitor de Amostra e verificar as estatísticas do último valor de lambda, e ainda obter qual é a minha previsão do último h. Obrigado a todos.
Bo Yu
11
@ Tomas, muito obrigado. Sim você está certo! O WinBUGS fornece a previsão de h [N], incluindo estatísticas e densidade de probabilidade. Agora eu entendi. Atenciosamente,
Bo Yu