Data augmentation

The linkage data example

1 The linkage data - Section 2 of Tanner and Wong (1987)

Martin A. Tanner and Wing Hung Wong (1987) The Calculation of Posterior Distributions by Data Augmentation, Journal of the American Statistical Association, Vol. 82, No. 398, pages 528-540. URL: https://www.jstor.org/stable/2289457

From a genetic linkage model, it is believed that 197 animals are distributed multinomially into four categories, ie. \[ y = (y_1,y_2,y_3,y_4)=(125,18,20,34), \] with cell probabilities specified by \[ \left(\frac{1}{2}+\frac{\theta}{4},\frac{1-\theta}{4},\frac{1-\theta}{4},\frac{\theta}{4}\right). \] The likelihood is of the form \[ p(y|\theta) \propto (2+\theta)^{y_1}(1-\theta)^{y_2+y_3}\theta^{y_4}, \] To illustrate the data augmentation algorithm, \(y\) is augmented by splitting the first cell into two cells, once of which having cell probability \(1/2\), the other having cell probability \(\theta/4\). Then the augmented data set is given by \[ x=(x_1,x_2,x_3,x_4,x_5), \] where \(x_1+x_2=125\), \(x_3=y_2\), \(x_4=y_3\), and \(x_5=y_4\). The augmented likelihood is of the form \[ p(x|\theta) \propto \theta^{x_2+x_5}(1-\theta)^{x_3+x_4}. \] When the prior of \(\theta\) is \(U(0,1)\), it follows that \(\theta|x \sim Beta(x_2+x_5+1,x_3+x_4+1)\). Also, the distribution of \(x_2|\theta \sim Binomial(y_1,\theta/(\theta+2))\).

1.1 Data augmentation algorithm

For given initial value, \(\theta^{(0)}\), sample iteratively through the following distributions, for \(j=1,2,\ldots,M\), for large \(M\):

1. Sample \(x_2^{(j)}\) from \(Bin(y_1,\theta^{(j-1)}/(2+\theta^{(j-1)}))\)

2. Sample \(\theta^{(j)}\) from \(Beta(x_2^{(j)}+y_4+1,y_2+y_3+1)\)

The set of draws \(\{(x_2,\theta)^{(j)}\}_{j=1}^M\) approximates the augmented posterior \(p(x_2,\theta|y)\). In addition, \(\{\theta^{(j)}\}_{j=1}^M\) approximates \(p(\theta|y)\)

1.2 Implementation

Let \(\theta^{(0)}=1/2\), or any other values in \((0,1)\),

y = c(125,18,20,34)

theta = 0.5
M = 10000
draws = matrix(0,M,2)
for (iter in 1:M){
  x2 = rbinom(1,y[1],theta/(2+theta))
  theta = rbeta(1,x2+y[4]+1,y[2]+y[3]+1)
  draws[iter,] = c(theta,x2)
}
theta.d = draws[,1]
x2.d = draws[,2]

par(mfrow=c(1,3))
ts.plot(theta.d,xlab="iterations",ylab="",main=expression(theta))
acf(theta.d,main="")
hist(theta.d,prob=TRUE,xlab=expression(theta),main="")

quantile(theta.d,c(0.025,0.5,0.975))

##      2.5%       50%     97.5% 
## 0.5204466 0.6239384 0.7182968

2 Posterior summaries

I1 = mean(theta.d)
I2 = mean((x2.d+y[4]+1)/(x2.d+y[2]+y[3]+y[4]+2))
c(I1,I2)

## [1] 0.6229764 0.6226399

N = 1000
thetas = seq(0,1,length=N)
post = rep(0,N)
for (i in 1:N)
  post[i] = mean(dbeta(thetas[i],x2.d+y[4]+1,y[2]+y[3]+1))

hist(theta.d,prob=TRUE,xlab=expression(theta),main="")
lines(density(theta.d),col=1,lwd=2)
lines(thetas,post,col=2,lwd=2)
legend("topleft",legend=c("Kernel density approx.","Raoblackwell approx."),col=1:2,lwd=2,bty="n")