How to get prediction for a specific variable in WinBUGS?

Question

I am a new user of WinBUGS and have one question for your help. After running the following code, I got parameters of beta0 through beta4 (stats, density), but I don't know how to get the prediction of the last value of h, which I set to NA to model in the code.

Does anyone can given me a hint? Any advice would be greatly appreciated.

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))

@whuber yes, I think that's correct, but more fundamentally, you need to have things you want to predict (ie, have a posterior distribution for) be distinct from things you've already observed. You can either make the prediction explicitly in winbugs or in postprocessing by using the samples of the betas. — atiretoo, Jul 02 '12 at 22:47
@atiretoo As far as I can tell, the lambdas are exactly what one wants to predict: this is a generalized linear model for a Poisson distribution with log link and the lambdas are the predicted Poisson parameters. They have not been observed. I believe all one needs to do here is to set a monitor on lambda[N]. — whuber, Jul 03 '12 at 13:23
@whuber, I'd rather say monitor `h[N]` instead of `lambda[N]` ... and you get the posterior distribution of the predicted value. — Tomas, Jul 03 '12 at 14:10
@tomek but `h[N]` is not the predicted value: it will be a collection of draws from a set of predicted Poisson distributions. As such it combines variation in the Poisson parameters and variation from those Poisson distributions themselves. What is relevant is the posterior distribution of `lambda[N]`. — whuber, Jul 03 '12 at 14:28
@whuber, I agree with you that `h[N]` involves the Poisson variation, but I don't agree with you that it is not relevant. For predictions, you need exactly this - take the variation into account. This is achieved by looking at the posterior distribution of `h[N]`. `lambda[N]` will be the posterior distribution of the parameter (i.e. mean value), not the predicted value! — Tomas, Jul 03 '12 at 14:38
It is the same difference as between the confidence interval of a distribution vs confidence interval of the mean of that distribution. You need the first thing here. — Tomas, Jul 03 '12 at 14:40
Yes, @Tomas, if you interpret the question to be asking for *prediction intervals* rather than a *prediction* (which is usually interpreted in terms of the fit itself). Because the Poisson variation is trivial to calculate *exactly*, once you have the distribution of `lambda[N]`, you have everything you need for either interpretation of the question. So your suggestion to monitor only `h[N]`, although useful, is restricted in scope. — whuber, Jul 03 '12 at 14:46
@whuber, it won't be restricted. The distribution of `h[N]` tells the whole story, as its mean will converge to `mean(lambda[N])` for big number of MCMC samples. And you also have the credible intervals of the prediction (I guess the OP wants to predict the value, not its mean). — Tomas, Jul 03 '12 at 14:54
@Tomas, I fear you are incorrect. The code (at line 5) explicitly declares `h[N]` to have a Poisson distribution with mean `lambda[N]`. As the MCMC sample count rises, if there were no posterior variation in `lambda[N]`, this will converge to a `lambda[N]` distribution, not a constant. Due to the variation in `lambda[N]`, the distribution of `h[N]` will in fact converge to a mixture of Poissons. — whuber, Jul 03 '12 at 15:09
@whuber, maybe we misunderstood each other. I said that mean(h[N]) will converge to lambda[N] for big number of MCMC samples. This can be proven trivially, as h[N] ~ Poiss(lambda[N]), which means that expected mean(h[N]) must be equal to lambda[N]. — Tomas, Jul 03 '12 at 18:37

score 6 · Answer 1 · edited Apr 13 '17 at 12:44

6

Just add the variable h to the list of the parameters to be monitored. If you are using package like R2WinBUGS, then add variable h to to the list passed to parameters.to.save argument to the bugs function. Then look at your last value in h (the one with NA) - you will get a posterior distribution there.

This is usual way to make predictions in bayesian inference (see also this question). It is nice and simple! No more separation of parameter evaluation and prediction. Everything is done at once. The posterior distrubution of parameters is given by the actual data and propagated to the NA values (as "predictions").

edited Apr 13 '17 at 12:44

Community

1

answered Jul 03 '12 at 07:51

Tomas

5,735
11
52
93

Tomas, thanks for your help. I try to monitor the variable of h in Sample Monitor Tool but it doesn't work. Could you please help me again? The following is the procedure I did in WinBUGS (I do not know how to use R2WinBUGS): 1) select Sample in Sample Monitor Tool 2) type h in the white box marked node 3) click on the button marked set 4) h is not on the list of parameters I want to monitor, while other parameters (beta0, beta1, beta2, beta3, an p) are showed in the list. Do you know how can I add "h" to the list of parameters I want to monitor? Thank again! – Bo Yu Jul 03 '12 at 13:40
@BoYu, I don't know how to do it directly in WinBUGS since I run WinBUGS from R, using the R2WinBUGS package. It is much more practical because you can just save the R script and run it all as a batch, along with producing your own graphs etc. [Look here](http://stats.stackexchange.com/q/22842/5509) for example scripts. – Tomas Jul 03 '12 at 14:07
That said, it will for sure be also possible in WinBUGS itself, but I don't know how (and I guess most people call it from R). – Tomas Jul 03 '12 at 14:09
Tomas, thank you so much for your quick reply. I will try R2WinBUGS package. Best Regards, – Bo Yu Jul 03 '12 at 14:17
1

First of all, thank you, whuber, atiretoo, and Tomas! As whuber already mentioned, yes, it is a generalized linear model, the variable of h is fitted by Poisson distribution with varying rate (lambda) conditioned with different predictors (CF01, CF02, CF03, and IND). The last value of h is what I need to know and is not observed (marked as NA), while all other values of h are observed. I think whuber is right, I need to set lambda as a parameter in Sample Monitor Tool and check the stats of the last value of lambda, and further get what is my prediction of last h. Thanks all. – Bo Yu Jul 03 '12 at 15:24
@ whuber, as you already mentioned that h[N] is not the predicted value, so I need to calculate it based on the probability density of lambda corresponding to the last value of h. Is there any direct way for WinBUGS to provide users the estimates (stats, density) of the last value of h in my case? Thanks. – Bo Yu Jul 03 '12 at 17:49
@BoYu, I'm afraid you are maybe complicating it too much. You don't need to calculate anything in WinBUGS. Just purely use the posterior of h[N]... Maybe this link persuades whuber and you that this approach is correct: http://stats.stackexchange.com/q/22184/5509 – Tomas Jul 03 '12 at 18:43
@ Tomas, thanks for your help and for the link. I will look at the link and let you know my findings. Best Regards. – Bo Yu Jul 03 '12 at 18:53
1

@ Tomas, thank your so much. Yes, you are right! WinBUGS provides the prediction of h[N], including stats and probability density. I got it now. Best Regards, – Bo Yu Jul 03 '12 at 19:08
@BoYu, great - you've seen it yourself! Simplicity is powerful :-) – Tomas Jul 03 '12 at 20:07

How to get prediction for a specific variable in WinBUGS?

1 Answers1

Linked