Chapter 8 Exercises
Dugongs
Solutions
1. Write BUGS code to fit the dugongs data of Example 6.3.1, specifying appropriate flat priors for alpha, beta, gamma, and log(sigma). Obtain posterior summaries and densities for each of these parameters, a plot of the model fit, and an estimate of the number of effective parameters via the DIC Tool. How does this compare with the true number?
model {
for(i in 1:N) {
y[i] ~ dnorm(mu[i], inv.sigma2)
mu[i] <- alpha - beta*pow(gamma, x[i])
}
alpha ~ dunif(0, 100)
beta ~ dunif(0, 100)
gamma ~ dunif(0, 1)
inv.sigma2 <- 1/pow(sigma, 2)
log(sigma) <- log.sigma
log.sigma ~ dunif(-10, 10)
}
list(alpha = 3, beta = 2, gamma = 0.9, log.sigma = -5)
list(x = c(1.0, 1.5, 1.5, 1.5, 2.5, 4.0, 5.0, 5.0, 7.0,
8.0, 8.5, 9.0, 9.5, 9.5, 10.0, 12.0, 12.0, 13.0,
13.0, 14.5, 15.5, 15.5, 16.5, 17.0, 22.5, 29.0, 31.5),
y = c(1.80, 1.85, 1.87, 1.77, 2.02, 2.27, 2.15, 2.26, 2.47,
2.19, 2.26, 2.40, 2.39, 2.41, 2.50, 2.32, 2.32, 2.43,
2.47, 2.56, 2.65, 2.47, 2.64, 2.56, 2.70, 2.72, 2.57),
N = 27)
node mean sd MC error 2.5% median 97.5% start sample
alpha 2.651 0.07145 0.001808 2.529 2.645 2.807 1001 40000
beta 0.973 0.07641 9.745E-4 0.8249 0.9718 1.13 1001 40000
gamma 0.8619 0.03206 7.146E-4 0.7884 0.8652 0.9144 1001 40000
log.sigma -2.328 0.1485 0.001083 -2.596 -2.336 -2.016 1001 40000
Dbar = post.mean of -2logL; Dhat = -2LogL at post.mean of stochastic nodes
Dbar Dhat pD DIC
y -49.072 -52.790 3.719 -45.353
total -49.072 -52.790 3.719 -45.353
The true number of parameters is 4, and so pD is quite close.
2. Check the shape of the posterior distribution for gamma: do you think the uniform prior given to gamma is very influential? Try running again but giving gamma a prior more concentrated near one, say a Beta(8, 2). Does this make much difference? How might you summarise the effect of using a more 'informative' prior distribution?
Posterior mass is not too close to the limits (0, 1) imposed by the uniform prior, so the prior doesn't appear to be too influential...
model {
for(i in 1:N) {
y[i] ~ dnorm(mu[i], inv.sigma2)
mu[i] <- alpha - beta*pow(gamma, x[i])
}
alpha ~ dunif(0, 100)
beta ~ dunif(0, 100)
gamma ~
dbeta(8, 2)
gamma.prior ~ dbeta(8, 2) # allows sampling from prior
inv.sigma2 <- 1/pow(sigma, 2)
log(sigma) <- log.sigma
log.sigma ~ dunif(-10, 10)
}
node mean sd MC error 2.5% median 97.5% start sample
gamma 0.8637 0.03045 6.954E-4 0.7957 0.8667 0.9143 1001 40000
The Beta(8, 2) prior has very little effect on the posterior.
3. Try fitting a growth curve that is quadratic in age.
[Hint: it is helpful, for the sake of numerical stability, to standardise the x variable.]
Does pD get the number of parameters approximately right? Is the model preferable to the original model according to DIC? Plot the fitted line to examine the fit visually. Is this a sensible model?
model {
for(i in 1:N) {
y[i] ~ dnorm(mu[i], inv.sigma2)
z[i] <- (x[i] - mean(x[]))/sd(x[])
mu[i] <- alpha + beta*z[i] + gamma*z[i]*z[i]
}
alpha ~ dunif(
-100
, 100)
beta ~ dunif(
-100
, 100)
gamma ~ dunif(
-100
,
100
)
inv.sigma2 <- 1/pow(sigma, 2)
log(sigma) <- log.sigma
log.sigma ~ dunif(-10, 10)
}
list(alpha = 3, beta = 2, gamma = 0.9, log.sigma = -5)
list(x = c(1.0, 1.5, 1.5, 1.5, 2.5, 4.0, 5.0, 5.0, 7.0,
8.0, 8.5, 9.0, 9.5, 9.5, 10.0, 12.0, 12.0, 13.0,
13.0, 14.5, 15.5, 15.5, 16.5, 17.0, 22.5, 29.0, 31.5),
y = c(1.80, 1.85, 1.87, 1.77, 2.02, 2.27, 2.15, 2.26, 2.47,
2.19, 2.26, 2.40, 2.39, 2.41, 2.50, 2.32, 2.32, 2.43,
2.47, 2.56, 2.65, 2.47, 2.64, 2.56, 2.70, 2.72, 2.57),
N = 27)
node mean sd MC error 2.5% median 97.5% start sample
alpha 2.424 0.02555 2.055E-4 2.374 2.424 2.475 1001 40000
beta 0.3131 0.02565 1.928E-4 0.2625 0.3131 0.3637 1001 40000
gamma -0.09373 0.01603 1.439E-4 -0.1255 -0.09363 -0.06214 1001 40000
log.sigma -2.265 0.1494 8.838E-4 -2.538 -2.272 -1.952 1001 40000
Dbar = post.mean of -2logL; Dhat = -2LogL at post.mean of stochastic nodes
Dbar Dhat pD DIC
y -45.648 -49.672 4.024 -41.625
total -45.648 -49.672 4.024 -41.625
The value of pD is very close to the true number of parameters. The value of DIC is larger than previously (-45.353), suggesting that the quadratic model does not fit as well. The model fit looks good at first glance, but the model is not really sensible since it does not tend to an asymptote -- indeed, it predicts that dugong lengths start decreasing again after age ~25.