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

[exercises-ch8-dugongs-solutions0][exercises-ch8-dugongs-solutions1][exercises-ch8-dugongs-solutions2][exercises-ch8-dugongs-solutions3]

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.

[exercises-ch8-dugongs-solutions4]

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?

[exercises-ch8-dugongs-solutions5]
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)
}

[exercises-ch8-dugongs-solutions6][exercises-ch8-dugongs-solutions7]
   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   

[exercises-ch8-dugongs-solutions8]
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.