0
$\begingroup$

Can anyone give me a quick walkthrough of what this code is doing? I think the line

log( mu[i] ) <- alpha + beta * (year[i]-1970) + log( pop[i] )

is updating mu for each iteration and then fatal[i] and pred[i] sample from this. I am not clear on exactly what is the difference between the two. How does the data for fatal fit into the model? I know how the step function works but I am not clear on what is happening in the line

probpred[i] <- step(pred[i]-fatal[i]) 

Any help would be appreciated. I have had a look at the help manual but it is a little brief for my understanding. I am brand new to Bayesian analysis and WinBUGS.

model
{
    for (i in 1:5) {
        fatal[i] ~ dpois(mu[i])
        pred[i] ~ dpois(mu[i])
        probpred[i] <- step(pred[i]-fatal[i])
        log( mu[i] ) <- alpha + beta * (year[i]-1970) + log( pop[i] )
        }
        alpha ~ dnorm(0,0.001)
        beta  ~ dnorm(0,0.001)
}

list(alpha=5.74, beta=-0.03)

list(pop=c(124.07, 126.63, 130.67, 133.03, 135.04), 
fatal=c(3798, 3590, 3422, 3679, NA ),
year = c(1970, 1971, 1972, 1973, 1974 )  )   
$\endgroup$
0
$\begingroup$

Briefly;

the number of fatalities in a given year (the variable fatal) follows a poisson ditribution, which has a linear trend (alpha + beta*(year-1970)) with an extra hazard of the log of the size of the population.

Incidentally, in winBUGS, anything assigned with <- is deterministic; everything assigned with ~ is probabilistic.

The prior distribution for the linear parameters alpha and beta is a "vague" Normal distribution with variance 1000. Note that most distributions in winBUGS are parametrised differently to what you might be used to.

The stuff after the closing brace are the initial values to start the MCMC runs at, and the data on which the model is to be fit.

The fatal data influences mu via the likelihood.

If looks like probpred (and hence pred) was only created for tracking purposes, to see how often the actual # deaths was more or less than a random sample from the predictive distribution.

$\endgroup$
  • $\begingroup$ Ok, so fatal[i] and pred[i] are not sampling from the same distribution: fatal[i] samples from the posterior distribution, influenced by the data, and pred[i] samples from the posterior predictive distribution. Is that correct? What about the NA in fatal[5]? Will fatal[5] and pred[5] sample from the same distribution? Also, does probpred[i] <- step(pred[i]-fatal[i]) use the observed data of does it also sample from the posterior distribution? $\endgroup$ – Jon Aug 30 '16 at 9:49
  • $\begingroup$ That was probably poor wording on my part. fatal is data, it isn't sampled (in the computational sense) from anything. The line fatal[i] ~ dpois(mu[i]) is a statement about mu, not about fatal. pred on the other hand is sampled from that distribution. It often helps to draw out the DAG corresponding to a winBUGS program. I think there is a utility that does that for you. $\endgroup$ – JDL Aug 30 '16 at 10:00
  • $\begingroup$ So if I monitor fatal the missing value for 1976 will give the expected number of fatalities in that year and this will come from the posterior distribution which will use the observed data from previous years and the prior assumption for the parameters, and this is different to the posterior predicted number of fatalities in 1976. $\endgroup$ – Jon Aug 30 '16 at 10:25
  • $\begingroup$ yes, if you look at the trace for fatal[1976] after the MCMC has converged, it should provide samples from the posterior distribution of fatal[1976] given all of the other data. $\endgroup$ – JDL Aug 30 '16 at 10:27
  • $\begingroup$ OK, and then if you wanted to estimate the probability that the number of fatalities in 1976 was less than say 3000 you would construct a statement like: probpred[i] <- step(pred[i]-3000) and look at the value for probpred[5]. $\endgroup$ – Jon Aug 30 '16 at 10:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.