Bayesian Models for Ecologists

Markov Chain Monte Carlo 2: Multiparameter models, Gibbs and M-H updates

June 07, 2022

---

I. Introduction

If we seek to fit the two-parameter model described in the lecture to data \(y_i\) for \(i=1,\ldots,n\), we begin with a model statement written using statistical notation as

\[\begin{align} y_i &\sim \text{N}(\mu,\sigma^2) \;, \\ \mu &\sim \text{N}(\mu_0,\sigma^2_0) \;, \\ \sigma^2 &\sim \text{IG}(q,r) \;, \end{align}\]

The we write the full posterior distribution as

\[\begin{equation} [\mu,\sigma^2|\mathbf{y}]\propto \prod_{i=1}^n [y_i|\mu,\sigma^2][\mu][\sigma^2] \;. \notag \end{equation}\]

Then we find the full-conditional distributions for unknown parameters \(\mu\) and \(\sigma^2\). These full-conditional distributions are proportional to the multiplicative components right-hand side of the expression above so that

\[\begin{align} [\mu | \cdot] &\propto \prod_{i=1}^n[y_i|\mu,\sigma^2][\mu] \;, \\ &\propto \text{N}(\tilde\mu,\tilde\sigma^2) \;, \end{align}\]

where,

\[\begin{align} \tilde{\sigma}^2 &= \left(\frac{n}{\sigma^2}+\frac{1}{\sigma^2_0}\right)^{-1} \;, \\ \tilde{\mu} &= \tilde{\sigma}^2\left(\frac{\sum_{i=1}^n y_i}{\sigma^2}+\frac{\mu_0}{\sigma^2_0}\right) \;, \end{align}\]

and,

\[\begin{align} [\sigma^2 | \cdot] &\propto \prod_{i=1}^n[y_i|\mu,\sigma^2][\sigma^2] \;, \\ &\propto \text{IG}(\tilde{q},\tilde{r}) \;, \end{align}\]

where,

\[\begin{align} \tilde{q} &= \frac{n}{2}+q \;, \\ \tilde{r} &= \left(\frac{\sum_{i=1}^n (y_i-\mu)^2}{2}+r\right) \;. \end{align}\]

---

The MCMC algorithm

This model results in conjugacy for both \(\mu\) and \(\sigma^2\) because we can work out both full-conditional distributions by hand and the form of the full-conditional distributions matches the priors. Thus, our MCMC algorithm is a Gibbs sampler.

norm.mcmc <- function(y,n.mcmc){
 
  ###
  ###  Set up Variables
  ###

  n=length(y)
  mu.save=rep(0,n.mcmc)
  s2.save=rep(0,n.mcmc)
 
  ###
  ###  Priors and Starting Values
  ###

  mu.0=0
  s2.0=10000

  q=.001
  r=.001

  mu=mean(y)

  ###
  ###  MCMC Loop
  ###

  for(k in 1:n.mcmc){
    
    ###  
    ###  Sample s2
    ###
    
    q.tilde=n/2+q 
    r.tilde=sum((y-mu)^2)/2+r
    s2=1/rgamma(1,q.tilde,r.tilde)
    
    ###  
    ###  Sample mu
    ###
    
    s2.tilde=1/(n/s2+1/s2.0)
    mu.tilde=s2.tilde*(sum(y)/s2+mu.0/s2.0)
    mu=rnorm(1,mu.tilde,sqrt(s2.tilde))
    
    ###  
    ###  Save samples
    ###
    
    mu.save[k]=mu
    s2.save[k]=s2
  
  }
  
  ###  
  ###  Write Output
  ###

  list(mu.save=mu.save,s2.save=s2.save,mu.0=mu.0,s2.0=s2.0,q=q,r=r)

}

The data we will analyze are measured snout-vent-lengths (svl, in cm) of adult Tiger Snakes (Notechis scutatus) in the Christmas Island population in Australia. These are a subset of data that are well-described by Aubret (2012). According to Aubret (2012), the motivation behind this study was: ``Mean adult size has been used as the traditional measure of body size to explain trends of insular gigantism and dwarfism in a wide array of taxa. However, patterns of variation in body size at birth have received surprisingly little attention, leaving open the possibility that adult body-size differences are nonadaptive consequences of selection acting on neonate body size. Here I used an empirical and correlative approach to test this hypothesis in a mosaic of 12 island and mainland snake populations in Australia.’’

Aubret F. (2012). Body size evolution on islands: are adult size variations in tiger snakes a non-adaptive consequence of selection on birth size? , 179: 756-767.

Our goal is to estimate the tiger snake population mean svl and variance using the normal model. We can fit the model by running the algorithm with our data and look at the trace plots for the MCMC realizations of \(\mu\) and \(\sigma^2\).

y=scan("svl.dat")
out=norm.mcmc(y,10000)

layout(matrix(1:2,2,1))
plot(out$mu.save,type="l")
plot(out$s2.save,type="l")

The trace plots look good (i.e., well-mixed and converged) with \(10000\) MCMC iterations. We can now make inference and look at the marginal posterior histograms and the mean and 95% credible intervals for each parameter.

dIG<-function(x,q,r){ # our IG PDF
  (r^q)/gamma(q)*x^(-(q+1))*exp(-r/x)
}

mu.mean=mean(out$mu.save)
mu.l=quantile(out$mu.save,.025)
mu.u=quantile(out$mu.save,.975)

s2.mean=mean(out$s2.save)
s2.l=quantile(out$s2.save,.025)
s2.u=quantile(out$s2.save,.975)

layout(matrix(1:2,1,2))
hist(out$mu.save,breaks=40,prob=TRUE,xlab=bquote(mu),main="")
curve(dnorm(x,out$mu.0,sqrt(out$s2.0)),add=TRUE,col=2)
abline(v=mu.mean,col=3)
rug(c(mu.l,mu.u),col=3,lwd=2)
hist(out$s2.save,breaks=40,prob=TRUE,xlab=bquote(sigma^2),main="")
curve(dIG(x,out$q,out$r),add=TRUE,col=2)
abline(v=s2.mean,col=3)
rug(c(s2.l,s2.u),col=3,lwd=2)

---

Problem

1. Modify the MCMC algorithm above to accommodate the prior: \(\log(\sigma^2)\sim\text{N}(\mu_l,\sigma^2_l)\) and use a random-walk proposal distribution for \(\log(\sigma^2)\).

Answer

norm.MH.mcmc <- function(y,n.mcmc,s.tune){

  ###
  ###  Set up Variables
  ###

  n=length(y)
  mu.save=rep(0,n.mcmc)
  s2.save=rep(0,n.mcmc)
  acc=0

  ###
  ###  Priors and Starting Values
  ###

  mu.0=0
  s2.0=10000

  mu.l=0
  sig.l=1000

  mu=mean(y)
  s2=var(y)

  ###
  ###  MCMC Loop
  ###

  for(k in 1:n.mcmc){

    ###
    ###  Sample log(s2)
    ###

    s2.star=exp(rnorm(1,log(s2),s.tune))
    mh.1=sum(dnorm(y,mu,sqrt(s2.star),log=TRUE))+dnorm(log(s2.star),mu.l,sig.l,log=TRUE)
    mh.2=sum(dnorm(y,mu,sqrt(s2),log=TRUE))+dnorm(log(s2),mu.l,sig.l,log=TRUE)
    mh=exp(mh.1-mh.2)
    if(mh>runif(1)){
      s2=s2.star
      acc=acc+1
    }

    ###
    ###  Sample mu
    ###

    s2.tilde=1/(n/s2+1/s2.0)
    mu.tilde=s2.tilde*(sum(y)/s2+mu.0/s2.0)
    mu=rnorm(1,mu.tilde,sqrt(s2.tilde))

    ###
    ###  Save samples
    ###

    mu.save[k]=mu
    s2.save[k]=s2

  }

  ###
  ###  Write Output
  ###

  list(mu.save=mu.save,s2.save=s2.save,mu.0=mu.0,s2.0=s2.0,mu.l=mu.l,sig.l=sig.l,acc=acc)

}

2. Fit the new model to data using your MCMC algorithm with Metropolis-Hastings update for \(\log(\sigma^2)\). Tune the proposal distribution by adjusting the proposal standard deviation to yield approximately 40% acceptance rate and look at the trace plots for \(\mu\) and \(\sigma^2\).

Answer

out.MH=norm.MH.mcmc(y,10000,.5)
out.MH$acc/10000

## [1] 0.4079

layout(matrix(1:2,2,1))
plot(out.MH$mu.save,type="l")
plot(out.MH$s2.save,type="l")

3. Overlay marginal posterior smoothed density plots based on the output from fitting model 1 versus model 2 to compare the inference.

Answer

plot(density(out$s2.save),type="l",lwd=2,col=1,main="",xlab=bquote(sigma^2))
lines(density(out.MH$s2.save),type="l",lwd=2,col=2)
legend("topright",lwd=2,col=1:2,legend=c("Model 1","Model 2"))