3 Week 8: Gibbs Sampler – Data augmentation

Data augmentation is a useful tool for obtaining samples from the posterior distribution of parameters where $\pi (𝜽|𝐱)$ is not in a convenient form for analysis. The Gibbs sampler with data augmentation is primarily used in the following two situations:

1. 1.

   In real world problems there is often missing data. Although in theory the distribution of the observed data can be obtained by integrating out the missing data, this is often difficult or impossible to do.

2. 2.

   The likelihood of the data is not tractable in its current form, but conditional upon the collection of unobserved (extra) data, the likelihood becomes tractable. This corresponds to the situation we have observed with the EM algorithm and there are similarities (as well as obvious differences) between the EM algorithm and Gibbs sampler in data augmentation problems.

More generally, data augmentation is used extensively within MCMC (Markov chain Monte Carlo) algorithms to assist with obtaining samples from $\pi (𝜽|𝐱)$.

The generic set up is as follows. Let $𝐱$ denote the observed data and let $𝜽$ denote the model parameters. Then

$$\pi (𝜽|𝐱)\propto \pi (𝐱|𝜽)\times \pi (𝜽)$$

and this can be difficult to work with if the likelihood $\pi (𝐱|𝜽)$ is not in a convenient form. Now suppose that there is extra information $𝐲$ such that joint distribution of $𝐱$ and $𝐲$, $\pi (𝐱,𝐲|𝜽)$ is more convenient to work with. We can then look to construct a Gibbs sampler (or alternative MCMC algorithm) to obtain samples from $\pi (𝜽,𝐲|𝐱)$. It is then trivial to obtain samples from the marginal distribution $\pi (𝜽|𝐱)$ by simply ignoring the $𝐲$ values. In other words, focus on the marginal distribution $\pi (𝜽|𝐱)$. Note that

	$\pi (𝜽,𝐲|𝐱)$	$=$	${\displaystyle \frac{\pi (𝐱,𝐲|𝜽)\pi (𝜽)}{\pi (𝐱)}}$
		$\propto$	$\pi (𝐱,𝐲|𝜽)\pi (𝜽).$

Then the Gibbs sampler alternates between updating the parameters and augmented data.

1. 1.

   Update $𝜽$ given $𝐱$ and $𝐲$. ie. Use $\pi (𝜽|𝐱,𝐲)$.

2. 2.

   Update $𝐲$ given $𝐱$ and $𝜽$. ie. Use $\pi (𝐲|𝐱,𝜽)$.

Both the updates of $𝜽$ and $𝐲$ will often be broken down into a number of steps. Note the similarities to the EM algorithm with step 1 replacing the M-step (updating the parameters given the observed and augmented data) and step 2 replacing the E-step (updating the augmented data given the observed data and parameters).

We illustrate data augmentation within the Gibbs sampler using a range of examples.

We begin by illustrating data augmentation Gibbs sampling with a very simple example, which is a simplified version of the Normal mixtures introduced in the Week 6 notes for the EM algorithm.

Suppose that $X_1,X_2,\mathrm{\dots },X_n$ are iid (independent and identically distributed) from the mixture density

$$f(x)=\frac{1}{2\sqrt{2\pi }}\left(e^{-x^2/2}+e^{-{(x-\theta )}^2/2}\right).$$

That is,

Thus the mixture can be constructed as follows. For each observation toss a fair coin and let $Y$ denote the outcome of the coin toss with $Y=0$ if the coin shows a head and $Y=1$ otherwise. If the $i^{th}$ coin toss shows a head $(Y_i=0)$ draw $X_i$ from $N(0,1)$, otherwise draw $X_i$ from $N(\theta ,1)$.

The likelihood can be written explicitly:-

$$L(\theta ;𝐱)=\pi (𝐱|𝜽)=\prod _{i=1}^n\frac{1}{2\sqrt{2\pi }}\left(e^{-x_i^2/2}+e^{-{(x_i-\theta )}^2/2}\right),$$

where $𝐱=(x_1,x_2,\mathrm{\dots },x_n)$. It is difficult to maximise (for MLE) or to compute the posterior distribution for $\theta$ using the given likelihood. However, if we knew the results of the coin tosses (ie. If we knew $Y_1,Y_2,\mathrm{\dots }{Y}_n)$, then we would know which observations come from the $N(\theta ,1)$ and the posterior distribution (or MLE) would be straightforward to derive.

Let $𝐲=(y_1,y_2,\mathrm{\dots },y_n)$ denote the outcome of the $n$ coin tosses and let $𝒜=\{i;y_i=1\}$, the set of coin tosses that are tails. Thus for $i\in 𝒜$, $X_i$ is drawn from $N(\theta ,1)$.

Original likelihood:-

$L(\theta ;𝐱)={\displaystyle \prod _{i=1}^n}{\displaystyle \frac{1}{2\sqrt{2\pi }}}\left(e^{-x_i^2/2}+e^{-{(x_i-\theta )}^2/2}\right).$

Augmented likelihood:-

	$L(\theta ;𝐱,𝐲)=\pi (𝐱,𝐲|\theta )$	$=$	${\displaystyle \prod _{i\in 𝒜}}\left({\displaystyle \frac{1}{2}}\times {\displaystyle \frac{1}{\sqrt{2\pi }}}{e}^{-{(x_i-\theta )}^2/2}\right)\times {\displaystyle \prod _{i\in {𝒜}^C}}\left({\displaystyle \frac{1}{2}}\times {\displaystyle \frac{1}{\sqrt{2\pi }}}{e}^{-x_i^2/2}\right)$
		$\propto$	$\exp \left(-{\displaystyle \sum _{i\in 𝒜}}{(x_i-\theta )}^2/2\right).$

Therefore if we take a $N(0,1)$ prior for $\theta$, we have that

	$\pi (\theta |𝐱,𝐲)$	$\propto$	$L(\theta ;𝐱,𝐲)\times \pi (\theta )$
		$\propto$	$\exp \left(-{\displaystyle \sum _{i\in 𝒜}}{(x_i-\theta )}^2/2\right)\times \exp (-{\theta }^2/2)$
		$=$	$\exp \left(-{\displaystyle \frac{1}{2}}\left\{m{\theta }^2-2\theta {\displaystyle \sum _{i=1}^m}{x}_i+{\displaystyle \sum _{i\in 𝒜}}{x}_i^2\right\}\right)\times \exp \left(-{\displaystyle \frac{1}{2}}{\theta }^2\right)$
		$\propto$	$\exp \left(-{\displaystyle \frac{1}{2}}\left\{(m+1){\theta }^2-2\theta {\displaystyle \sum _{i\in 𝒜}}{x}_i\right\}\right)$

Thus the posterior distribution of $\theta$ is

$$N(\frac{1}{m+1}\sum _{i\in 𝒜}{x}_i,\frac{1}{m+1}),$$

where $m=|𝒜|$.

On the other hand let $Y_i$ denote the outcome of the $i^{th}$ coin toss. Then it is straightforward to show that

	$P(Y_i=1|\theta ,\text{everything else})$	$=$	${\displaystyle \frac{P(Y_i=1,x_i)}{f(x_i)}}$
		$=$	${\displaystyle \frac{P(Y_i=1,x_i)}{P(Y_i=1,x_i)+P(Y_i=0,x_i)}}$
		$=$	${\displaystyle \frac{\exp (-{(x_i-\theta )}^2/2)/(2\sqrt{2\pi })}{\exp (-{(x_i-\theta )}^2/2)/(2\sqrt{2\pi })+\exp (-x_i^2/2)/(2\sqrt{2\pi })}}$
		$=$	${\displaystyle \frac{\exp (-{(x_i-\theta )}^2/2)}{\exp (-{(x_i-\theta )}^2/2)+\exp (-x_i^2/2)}}.$

The Gibbs sampling algorithm is:-

1. 1.

   Initial value for ${\theta }_0$. (A reasonable starting value is ${\theta }_0=\frac{1}{n}{\sum }_{i=1}^n{x}_i$, the sample mean.)

2. 2.

   For $i=1,2,\mathrm{\dots },n$, set $Y_i=1$ with probability $\frac{\exp (-{(x_i-\theta )}^2/2)}{\exp (-{(x_i-\theta )}^2/2)+\exp (-x_i^2/2)}$, and set $Y_i=0$ otherwise. This updates the set $𝒜$.

3. 3.

   Sample $\theta \sim N(\frac{1}{m+1}{\sum }_{i\in 𝒜}{x}_i,\frac{1}{m+1})$.

We revisit the genetic linkage example from Week 6, the EM algorithm. Whilst we are now focused on the posterior distribution rather than the MLE of the parameter $\theta$, we use exactly the same data augmentation as before.

The data consist of the genetic linkage of 197 animals and the data are divided into four genetic categories, labeled 1 through to 4, see [1]. The probabilities that an animal belongs to each of the four categories are $\frac{1}{2}+\frac{\theta }{4},\frac{1-\theta }{4},\frac{1-\theta }{4},\frac{\theta }{4}$, respectively. Let $𝐱=(x_1,x_2,x_3,x_4)=(125,18,20,34)$ denote the total number of animals in each category.

We are interested in estimating $\theta$. Note that for the category probabilities to be valid we require that . Therefore we shall assign an (uninformative) uniform prior to $\theta$. This is a multinomial experiment (4 different outcomes), so

	$L(\theta ;𝐱)$	$=$	${\displaystyle \frac{(x_1+x_2+x_3+x_4)!}{x_1!x_2!x_3!x_4!}}{p}_1^{x_1}{p}_2^{x_2}{p}_3^{x_3}{p}_4^{x_4}$
		$=$	${\displaystyle \frac{197!}{125!18!20!34!}}{\left({\displaystyle \frac{1}{2}}+{\displaystyle \frac{\theta }{4}}\right)}^{125}{\left({\displaystyle \frac{1}{4}}(1-\theta )\right)}^{18}{\left({\displaystyle \frac{1}{4}}(1-\theta )\right)}^{20}{\left({\displaystyle \frac{\theta }{4}}\right)}^{34}.$

The posterior distribution of $\theta$ is not of a particularly nice form,

$$p(\theta |𝐱)\propto {(2+\theta )}^{x_1}{(1-\theta )}^{x_2+x_3}{\theta }^{x_4}.$$

Suppose that the observed cell $x_1$ could be divided into two subcategories $A$ and $B$. Suppose that $y$ $(x_1-y)$ is the number of animals in subcategory $A$ $(B)$ with cell probability $p_A=\theta /4$ $(p_B=1/2)$. This would give an augmented data set $(𝐱,y)$. Then

	$L(\theta ;𝐱,y)$	$=$	${\displaystyle \frac{(x_1+x_2+x_3+x_4)!}{(x_1-y)!y!x_2!x_3!x_4!}}{p}_A^y{p}_B^{x_1-y}{p}_2^{x_2}{p}_3^{x_3}{p}_4^{x_4}$
		$\propto$	${\left({\displaystyle \frac{\theta }{4}}\right)}^y{\left({\displaystyle \frac{1}{4}}(1-\theta )\right)}^{18}{\left({\displaystyle \frac{1}{4}}(1-\theta )\right)}^{20}{\left({\displaystyle \frac{\theta }{4}}\right)}^{34}.$

Thus

$$p(\theta |𝐱,y)\propto {(1-\theta )}^{38}{\theta }^{34+y}.$$

Therefore the posterior density is proportional to a Beta density giving

$$\theta |𝐱,y\sim \mathrm{Beta}(y+35,39).$$

What is $y\mathrm{|}\theta \mathrm{,}𝐱$?

There are 125 $(=x_1)$ animals in categories $A$ and $B$.

The conditional probability that the animal belongs to category $A$ given that the animal belongs to category 1 is

	$P(\text{Belongs to category }A|\text{Belongs to category 1})$	$=$	${\displaystyle \frac{P(\text{Belongs to category }A\text{ and to category 1})}{P(\text{Belongs to category 1})}}$
		$=$	${\displaystyle \frac{P(\text{Belongs to category }A)}{P(\text{Belongs to category 1})}}$
		$=$	${\displaystyle \frac{\theta /4}{1/2+\theta /4}}={\displaystyle \frac{\theta }{2+\theta }}.$

Thus $y\sim Bin(125,\theta /(2+\theta ))$.

Therefore the Gibbs sampler iterates between the following two equations:

	$\theta |𝐱,y$	$\sim$	$Beta(y+x_4+1,x_2+x_3+1)=\mathrm{Beta}(y+35,39)$
	$y|\theta ,𝐱$	$\sim$	$\mathrm{Bin}(x_1,{\displaystyle \frac{\theta }{\theta +2}})=\mathrm{Bin}(125,{\displaystyle \frac{\theta }{\theta +2}})$

(Note that the above is for general $𝐱$, not just the data given.)

A sample of size 100 from the posterior distribution was obtained using the above Gibbs sampler and the output is presented below. Initial values $\theta =0.5$ and $y=25$ were chosen.

Unnumbered Figure: Link

The algorithm appears to converge immediately, but just to be safe I shall disregard the first 20 iterations as burn-in. The following summary statistics are then obtained for $\theta$:-
Estimate of posterior mean, 0.6258.
Estimate of posterior variance, $0.00281(={0.053}^2)$.
Posterior density plots (obtained by kernel smoothing) suggest a modal value of about 0.645.

It is important to note that the sequence of realisations of $\theta$ are not independent. This is clear from their construction, via the Markov chain. One useful summary of dependence is an ACF (autocorrelation function) plot. The acf plot (in R) gives the estimated value of $\text{Corr}[{\theta }_0,{\theta }_k]$ for $k=0,1,\mathrm{\dots }$, the correlation between samples from the posterior distribution which differ by a lag of $k$ iterations. For independent samples, $\text{Corr}[{\theta }_0,{\theta }_k]=0$ for all $k\ge 1$, although an estimate of the correlation will not be exactly 0. Generally for MCMC, the faster the $\text{Corr}[{\theta }_0,{\theta }_k]$ approaches 0 as $k$ increases the better. (With very few exceptions, MCMC constructs Markov chains with positive correlation, $\text{Corr}[{\theta }_0,{\theta }_k]>0$.)

An acf plot for $\theta$ is given based upon a sample of 1000 iterations following a burn-in of 100 iterations. As we can see there is very little dependence with only $\text{Corr}[{\theta }_0,{\theta }_1]$ significantly different to 0. (Note that $\text{Corr}[{\theta }_0,{\theta }_0]$ is always equal to 1.) This is very good. In practice the dependence is usually a lot more marked. We will say more about dependence and the cost of dependent versus independent samples in later weeks.

Unnumbered Figure: Link

The performance of students in a test is thought to depend upon individual and school variability. In particular the following model is suggested:-

	$X_i$	$\sim$	$N(\theta ,1/\tau )\mathit{\hspace{1em}\hspace{1em}}(i=1,2,\mathrm{\dots },m)$
	$Y_{ij}$	$\sim$	$N(X_i,1/\lambda )\mathit{\hspace{1em}\hspace{1em}}(i=1,2,\mathrm{\dots },m;j=1,2,\mathrm{\dots },n_i),$

where $Y_{ij}$ is the performance of the $j^{th}$ individual from school $i$ and $X_i$ denotes the $i^{th}$ school effect. The parameters $\theta ,\lambda$ and $\tau$ are assumed to be unknown parameters to be estimated. Whilst $\{Y_{ij}\}$ are observed, the $\{X_i\}$ are unobserved. Therefore we use data augmentation (augment with $𝐱=(x_1,x_2,\mathrm{\dots },x_m)$) to estimate $(\theta ,\lambda ,\tau )$ based upon observations $𝐲=({𝐲}_1,{𝐲}_2,\mathrm{\dots },{𝐲}_m)$, where ${𝐲}_i=(y_{i1},y_{i2},\mathrm{\dots },y_{i{n}_i})$.

This is a hierarchical model with students based within schools and the school having an effect on student performance. We could have more hierarchical levels, for examples, the schools could be in different cities and we could have a city effect included.

Note that given $x_i$, ${𝐲}_i$ is independent of $\theta$ and $\tau$. Moreover, given $x_i$, $y_{ij}$ and $y_{ik}$ are independent, for $j\ne k$. Therefore

	$f({𝐲}_i|x_i,\theta ,\tau ,\lambda )$	$=$	${\displaystyle \prod _{j=1}^{n_i}}f(y_{ij}|x_i,\lambda )$		(3.1)
		$=$	${\displaystyle \prod _{j=1}^{n_i}}\sqrt{{\displaystyle \frac{\lambda }{2\pi }}}\exp \left(-{\displaystyle \frac{\lambda }{2}}{(y_{ij}-x_i)}^2\right).$

Now

$f(x_i|\theta ,\tau )$

$=$

$\sqrt{{\displaystyle \frac{\tau }{2\pi }}}\exp \left(-{\displaystyle \frac{\tau }{2}}{(x_i-\theta )}^2\right).$

(3.2)

Since the performance of students from different schools are independent, we have that

$L(\theta ,\tau ,\lambda ;𝐱,𝐲)$

$=$

${\displaystyle \prod _{i=1}^m}\left\{\sqrt{{\displaystyle \frac{\tau }{2\pi }}}\exp \left(-{\displaystyle \frac{\tau }{2}}{(x_i-\theta )}^2\right)\times {\displaystyle \prod _{j=1}^{n_i}}\sqrt{{\displaystyle \frac{\lambda }{2\pi }}}\exp \left(-{\displaystyle \frac{\lambda }{2}}{(y_{ij}-x_i)}^2\right)\right\}.$

(3.3)

To find the conditional distributions of $\theta ,\tau ,\lambda ,x_1,x_2,\mathrm{\dots },x_m$, we combine the likelihood (3.3) with the priors for $\theta$, $\tau$ and $\lambda$. For simplicity, we shall assume that the priors on $\theta$, $\tau$ and $\lambda$ are improper uniform priors. That is, $\pi (\theta )\propto 1$ $(\theta \in ℝ)$ and $\pi (\lambda ),\pi (\tau )\propto 1$ $(\lambda ,\tau \in {ℝ}^{+})$. Therefore, letting $\overline{x}=\frac{1}{m}{\sum }_{i=1}^m{x}_i$ and $N={\sum }_{i=1}^m{n}_i$, we have that

	$f(\tau |\lambda ,\theta ,𝐱,𝐲)$	$\propto$	${\tau }^{m/2}\exp \left(-\tau \times {\displaystyle \frac{1}{2}}{\displaystyle \sum _{i=1}^m}{(x_i-\theta )}^2\right)$
	$f(\theta |\lambda ,\tau ,𝐱,𝐲)$	$\propto$	$\exp \left(-\tau \times {\displaystyle \frac{1}{2}}{\displaystyle \sum _{i=1}^m}{(x_i-\theta )}^2\right)$
		$\propto$	$\exp \left(-\tau \times {\displaystyle \frac{1}{2}}(m{\theta }^2-2\theta m\overline{x})\right)$
	$f(\lambda |\tau ,\theta ,𝐱,𝐲)$	$\propto$	${\lambda }^{N/2}\exp \left(-\lambda \times {\displaystyle \frac{1}{2}}{\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{j=1}^{n_i}}{(y_{ij}-x_i)}^2\right).$

Therefore

	$\tau |\lambda ,\theta ,𝐱,𝐲$	$\sim$	$Gamma({\displaystyle \frac{m}{2}}+1,{\displaystyle \frac{1}{2}}{\displaystyle \sum _{i=1}^m}{(x_i-\theta )}^2)$
	$\theta |\lambda ,\tau ,𝐱,𝐲$	$\sim$	$N(\overline{x},{\displaystyle \frac{1}{m\tau }})$
	$\lambda |\tau ,\theta ,𝐱,𝐲$	$\sim$	$Gamma({\displaystyle \frac{N}{2}}+1,{\displaystyle \frac{1}{2}}{\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{j=1}^{n_i}}{(y_{ij}-x_i)}^2).$

All that is now required is the conditional distributions of $x_1,x_2,\mathrm{\dots },x_m$. Note that the $x_i$’s are independent with

	$f(x_i|\lambda ,\theta ,\tau ,{𝐱}_{i-},𝐲)$	$=$	$f(x_i|\lambda ,\theta ,\tau ,{𝐲}_i)$
		$\propto$	$\exp \left(-\tau \times {\displaystyle \frac{1}{2}}{(x_i-\theta )}^2\right)\times \exp \left(-\lambda \times {\displaystyle \frac{1}{2}}{\displaystyle \sum _{j=1}^{n_i}}{(y_{ij}-x_i)}^2\right)$
		$\propto$	$\exp \left(-{\displaystyle \frac{1}{2}}\left((\tau +n_i\lambda )x_i^2-2x_i\left(\tau \theta +\lambda {\displaystyle \sum _{j=1}^{n_i}}{y}_{ij}\right)\right)\right).$

Therefore, we have that

$$x_i|\lambda ,\theta ,\tau ,{𝐱}_{i-},𝐲\sim N(\frac{\tau \theta +\lambda {\sum }_{j=1}^{n_i}{y}_{ij}}{\tau +n_i\lambda },\frac{1}{\tau +n_i\lambda }).$$

Note that the mean of $x_i$ is a compromise between the mean of the data and the prior mean with

$$\frac{\tau \theta +\lambda {\sum }_{j=1}^{n_i}{y}_{ij}}{\tau +n_i\lambda }=\frac{\tau }{\tau +n_i\lambda }\theta +\frac{n_i\lambda }{\tau +n_i\lambda }{\overline{y}}_{i\cdot }.$$

Furthermore as more data are collected it moves closer to the mean of ${𝐲}_i$.

Example.

Suppose that a study involves 4 schools. Students from each school take a test and score a mark out of 100. The results are presented in the table below.

Note: The data were generated from; $X_i\sim N(60,4^2)$ and $Y_{ij}\sim N(X_i,2^2)$ with the $Y_{ij}$ rounded to give integer values. Therefore the ’true’ values are $\theta =60$, $\tau =0.05625$ and $\lambda =0.25$.

A Gibbs sampler algorithm was run on the above data set to get a sample of 1100 iterations. After discarding the first 100 iterations, we can use the remaining 1000 iterations to estimate the posterior means and standard deviations of the parameters.

The mean of $\tau$ is almost 5 times the true value of $\tau$. The reason for this is that there is only four schools and the means for the 4 schools are fairly similar. ($\tau$ large implies greater precision, ie. smaller variance.) A time series plot of the 1100 iterations of $\tau$ are given in the figure below. We also see that $\tau$ has a right-skewed distribution so that the mean is (quite a bit) larger than the posterior mode.

Unnumbered Figure: Link

Unnumbered Figure: Link

We consider a probit regression example. Let $y_i\in \{0,1\},i=1,\mathrm{\dots },n$ be binary response variables for a collection of $n$ objects with associated covariate measurement ${𝐱}_i$. We assume that

	$y_i$	$\sim$	$\text{Bernoulli }\left(\mathrm{\Phi }({\eta }_i)\right)$
	${\eta }_i$	$=$	${𝐱}_i^T𝜷$

where $\mathrm{\Phi }(\cdot )$ is the standard normal distribution function, ${\eta }_i$ is the linear predictor and $𝜷$ represents a $(p\times 1)$ column vector of regression coefficients. We are interested in the posterior distribution of $𝜷$ and we assume a $N({𝝁}_0,{𝐂}_0)$ prior for $𝜷$.

In the probit model, the mean (the probability of a $1$) is given by ${\mu }_i=\mathrm{\Phi }({\eta }_i)$, hence, it corresponds to the probit link function: $g({\mu }_i)={\mathrm{\Phi }}^{-1}({\mu }_i)$.

Note that

$$L(𝜷;𝐲)=\prod _{i=1}^n\mathrm{\Phi }{({𝐱}_i^T𝜷)}^{y_i}{\{1-\mathrm{\Phi }({𝐱}_i^T𝜷)\}}^{1-y_i}.$$

Therefore we do not have a nice analytical expression for $\pi (𝜷|𝐲)\propto L(𝜷;𝐲)$ and resort again to data augmentation.

For the data augmentation it is useful to think how we could simulate data $𝐲$ from the probit model with coefficients $𝜷$ and covariates ${𝐱}_1,\mathrm{\dots },{𝐱}_n$. To simulate $y_i$ such that $P(y_i=1)=\mathrm{\Phi }({𝐱}_i^T𝜷)$, we could take the following steps:

1. 1.

   Simulate $z_i\sim N({𝐱}_i^T𝜷,1)$.

2. 2.

   If $z_i>0$, set $y_i=1$ and otherwise set $y_i=0$.

Note that

	$P(y_i=1)$	$=$	$P(z_i>0)=P(N({𝐱}_i^T𝜷,1)>0)$
		$=$	$P(N(0,1)>-{𝐱}_i^T𝜷)=\mathrm{\Phi }({𝐱}_i^T𝜷),$

as required. Given the above construction, the $z_i$’s are simple Normal random variables and the $y_i$’s are a simple deterministic functions of the $z_i$’s.

Gibbs sampler
This representation lends itself to efficient simulation using Gibbs.
The joint posterior of $(𝜷,z_1,\mathrm{\dots },z_n)$ is given by

	$\pi (𝜷,z_1,\mathrm{\dots },z_n|𝐲)$	$\propto$	${\displaystyle \prod _{i:y_i=1}}\exp [-{\displaystyle \frac{1}{2}}{(z_i-𝐱_i{}^T𝜷)}^2]I[z_i>0]$
		$\times$
		$\times$	$\exp \left[-{\displaystyle \frac{1}{2}}{(𝜷-{𝝁}_0)}^T{𝐂}_0^{-1}(𝜷-{𝝁}_0)\right].$

The conditional posterior distribution of $𝜷$ is therefore multivariate normal:

$$𝜷|𝐳,𝐲\sim N({({𝐂}_0^{-1}+{𝐗}^T𝐗)}^{-1}({𝐂}_0^{-1}{𝝁}_0+{𝐗}^T𝐳),{({𝐂}_0^{-1}+{𝐗}^T𝐗)}^{-1}),$$

where $𝐳={(z_1,\mathrm{\dots },z_n)}^T$. This is similar to the block update of $𝜷$ in the linear regression in the Week 7 lab session. On the other the conditional posterior distribution for each $z_i$ is truncated normal,

$z_i|𝜷,𝐲$

$\propto$

(3.4)

We can therefore create a Gibbs sampler which, at each step

• •

  Samples from $𝜷|𝐳,𝐲$.

• •

  Samples from $z_i|𝜷,𝐲$, for $i=1,\mathrm{\dots },n$.

Example.

Response: Occurrence or non-occurrence of infection following birth by Caesarian section.

Covariates:

$\bullet$ $x_1=1$ if Caesarian section was not planned and $x_1=0$ otherwise;

$\bullet$ $x_2=1$ if there is presence of one or more risk factors, such as diabetes or excessive weight, and $x_2=0$ otherwise;

$\bullet$ $x_3=1$ if antibiotics were given as prophylaxis and $x_3=0$ otherwise.

In total: $251$ births.

The program for implementing the probit model is available in probit.r. I ran the code for 5500 iterations discarding the first 500 iterations as burn-in. The tables below give the estimates of the parameter means and standard deviations along with the probability that a parameter is positive (based on the sample). We can see from the table that the Caesarian section not being planned (emergency) and having risk factors increase the chance of infection, whilst antibiotics (fortunately) reduce the risk of infection.

Coefficients:	mean	Std	probability
(Intercept)	-0.9854	0.2049	0
noplan	0.5063	0.2313	0.9868
factor	1.0607	0.2386	1
antib	-1.7498	0.2447	0

Finally, in the table below, we give the (estimated) posterior probability of becoming infected for each set of covariates. Note that

$P(\text{Infected}|𝐱)$

$=$

${\displaystyle \int \mathrm{\Phi }({𝐱}^T𝜷)\pi (𝜷|𝐲)𝑑𝜷}.$

This can be estimated using $N$ samples from $𝜷|𝐲$ by

${\displaystyle \frac{1}{N}}{\displaystyle \sum _{j=1}^N}\mathrm{\Phi }({𝐱}^T{𝜷}_j).$

In the table below we compare the posterior probability of infections $(E)$ versus the observed proportion infected $(O)$ for each set of covariates. The table shows generally good agreement between the observed proportions and those given by the model.

We have introduced hierarchical models above and here we consider an alternative parameterisation of the hierarchical Gaussian model. Reparameterisation can be a useful in statistics to improve estimation of the model parameters. The reparameterisations presented below are based upon ideas from [2]. Let $Y=(Y_{11},\mathrm{\dots },Y_{m{n}_m})$ be observed and $X=(X_1,\mathrm{\dots },X_m)$ be unobserved, with

	$Y_{ij}$	$=$	$X_i+{\sigma }_y{ϵ}_{ij},j=1,2,\mathrm{\dots },n_i$
	$X_i$	$=$	$\theta +{\sigma }_x{z}_i,i=1,2,\mathrm{\dots },m,$

where ${\sigma }_x$ and ${\sigma }_y$ are assumed known and ${ϵ}_{ij}$ and $z_i$ are independent standard normal random variables ($N(0,1)$). The parameter $\theta$ is assumed to be unknown and its posterior distribution is of interest. We shall consider the case where $n_i=1$ for all $i$ and consequently drop the subscript $j$. Therefore the model can be rewritten as

	$Y_i$	$=$	$N(X_i,{\sigma }_y^2)$
	$X_i$	$=$	$N(\theta ,{\sigma }_x^2),i=1,2,\mathrm{\dots },m.$

Therefore this is a simplified version of the hierarchical model studied above.

An alternative parameterisation is the non-centered parameterisation proposed by [2],

	$Y_i$	$=$	$\theta +{\tilde{X}}_i+{\sigma }_y{ϵ}_i,$
	${\tilde{X}}_i$	$=$	${\sigma }_x{z}_i,i=1,2,\mathrm{\dots },m.$

Note that $\tilde{X}=({\tilde{X}}_1,\mathrm{\dots },{\tilde{X}}_m)$ and $\theta$ are a priori independent, but conditional upon the data are dependent.

Both (centered and non-centered) parameterisations permit a (data-augmentation) Gibbs sampler as outlined below. (Throughout we shall assign an improper, uniform prior to $\theta$; $\pi (\theta )\propto 1$.)

Centered algorithm

The joint distribution of $(\theta ,𝐱)$ given $𝐲$, $\pi (\theta ,𝐱|𝐲)$ satisfies

	$\pi (\theta ,𝐱|𝐲)$	$=$	${\displaystyle \frac{\pi (𝐲|𝐱,\theta )\pi (𝐱|\theta )\pi (\theta )}{\pi (𝐲)}}$
		$\propto$	$\pi (𝐲|𝐱)\pi (𝐱|\theta ),$

since $\pi (\theta )=1$ and given $𝐱$, $𝐲$ is independent of $\theta$.

Thus

$\pi (\theta ,𝐱|𝐲)$

$\propto$

${\displaystyle \prod _{i=1}^m}{\displaystyle \frac{1}{\sqrt{2\pi }{\sigma }_Y}}\exp \left(-{\displaystyle \frac{1}{2{\sigma }_Y^2}}{(y_i-x_i)}^2\right)\times {\displaystyle \frac{1}{\sqrt{2\pi }{\sigma }_X}}\exp \left(-{\displaystyle \frac{1}{2{\sigma }_X^2}}{(x_i-\theta )}^2\right).$

Therefore

	$\pi (\theta |𝐱,𝐲)$	$\propto$	${\displaystyle \prod _{i=1}^m}\exp \left(-{\displaystyle \frac{1}{2{\sigma }_X^2}}{(x_i-\theta )}^2\right)$
		$=$	$\exp \left(-{\displaystyle \frac{1}{2{\sigma }_X^2}}{\displaystyle \sum _{i=1}^m}{(x_i-\theta )}^2\right)$
		$=$	$\exp \left(-{\displaystyle \frac{1}{2{\sigma }_X^2}}\left\{{\displaystyle \sum _{i=1}^m}{x}_i^2-2{\displaystyle \sum _{i=1}^m}{x}_i\theta +m{\theta }^2\right\}\right)$
		$\propto$	$\exp \left(-{\displaystyle \frac{m}{2{\sigma }_X^2}}\left\{{\theta }^2-2\overline{x}\theta \right\}\right).$

This gives

$$\theta |𝐱,𝐲\sim N(\overline{x},\frac{{\sigma }_X^2}{m}).$$

(Note that if $W$ has pdf $f(w)\propto \exp \left(-\frac{1}{2\beta }(w^2-2\alpha w)\right)$, then $W\sim N(\alpha ,\beta )$.)

For $i=1,2,\mathrm{\dots },n$, $x_i$ depends only upon $\theta$ and $y_i$. Therefore

	$\pi (x_i|\theta ,y_i)$	$\propto$	$\exp \left(-{\displaystyle \frac{1}{2{\sigma }_Y^2}}{(y_i-x_i)}^2\right)\times \exp \left(-{\displaystyle \frac{1}{2{\sigma }_X^2}}{(x_i-\theta )}^2\right)$
		$\propto$	$\exp \left(-{\displaystyle \frac{1}{2{\sigma }_Y^2}}(x_i^2-2y_i{x}_i)\right)\times \exp \left(-{\displaystyle \frac{1}{2{\sigma }_X^2}}(x_i^2-2\theta x_i)\right)$
		$\propto$	$\exp \left(-{\displaystyle \frac{1}{2}}\left\{\left({\displaystyle \frac{1}{{\sigma }_Y^2}}+{\displaystyle \frac{1}{{\sigma }_X^2}}\right)x_i^2-2x_i\left({\displaystyle \frac{y_i}{{\sigma }_Y^2}}+{\displaystyle \frac{\theta }{{\sigma }_X^2}}\right)\right\}\right).$

This gives

$$x_i|\theta ,y_i\sim N(\left\{\frac{y_i}{{\sigma }_Y^2}+\frac{\theta }{{\sigma }_X^2}\right\}/\left\{\frac{1}{{\sigma }_Y^2}+\frac{1}{{\sigma }_X^2}\right\},1/\left\{\frac{1}{{\sigma }_Y^2}+\frac{1}{{\sigma }_X^2}\right\}).$$

It is straightforward using the above conditional distributions to construct a Gibbs sampling algorithm which alternates between updating $\theta |𝐱,𝐲$ and updating $𝐱|\theta ,𝐲$.

Non-centered algorithm

The joint distribution of $(\theta ,\widetilde{𝐱})$ given $𝐲$, $\pi (\theta ,\widetilde{𝐱}|𝐲)$ satisfies

	$\pi (\theta ,\widetilde{𝐱}|𝐲)$	$=$	${\displaystyle \frac{\pi (𝐲|\widetilde{𝐱},\theta )\pi (\widetilde{𝐱}|\theta )\pi (\theta )}{\pi (𝐲)}}$
		$\propto$	$\pi (𝐲|\widetilde{𝐱},\theta )\pi (\widetilde{𝐱}),$

since $\pi (\theta )=1$ and $\widetilde{𝐱}$ is a priori independent of $\theta$.

Thus

$\pi (\theta ,\widetilde{𝐱}|𝐲)$

$\propto$

${\displaystyle \prod _{i=1}^m}{\displaystyle \frac{1}{\sqrt{2\pi }{\sigma }_Y}}\exp \left(-{\displaystyle \frac{1}{2{\sigma }_Y^2}}{(y_i-\theta -{\tilde{x}}_i)}^2\right)\times {\displaystyle \frac{1}{\sqrt{2\pi }{\sigma }_X}}\exp \left(-{\displaystyle \frac{1}{2{\sigma }_X^2}}{\tilde{x}}_i^2\right).$

Therefore, letting $\mu =\frac{1}{n}{\sum }_{i=1}^n{\tilde{x}}_i$, we have that

	$\pi (\theta |\widetilde{𝐱},𝐲)$	$\propto$	${\displaystyle \prod _{i=1}^m}\exp \left(-{\displaystyle \frac{1}{2{\sigma }_Y^2}}{(y_i-{\tilde{x}}_i-\theta )}^2\right)$
		$=$	$\exp \left(-{\displaystyle \frac{1}{2{\sigma }_Y^2}}{\displaystyle \sum _{i=1}^m}{(y_i-{\tilde{x}}_i-\theta )}^2\right)$
		$=$	$\exp \left(-{\displaystyle \frac{1}{2{\sigma }_Y^2}}\left\{{\displaystyle \sum _{i=1}^m}{(y_i-{\tilde{x}}_i)}^2-2{\displaystyle \sum _{i=1}^m}(y_i-{\tilde{x}}_i)\theta +m{\theta }^2\right\}\right)$
		$\propto$	$\exp \left(-{\displaystyle \frac{m}{2{\sigma }_Y^2}}\left\{{\theta }^2-2(\overline{y}-\mu )\theta \right\}\right).$

This gives

$$\theta |\widetilde{𝐱},𝐲\sim N(\overline{y}-\mu ,\frac{{\sigma }_Y^2}{m}).$$

For $i=1,2,\mathrm{\dots },n$, ${\tilde{x}}_i$ depends only upon $\theta$ and $y_i$. Therefore

	$\pi ({\tilde{x}}_i|\theta ,y_i)$	$\propto$	$\exp \left(-{\displaystyle \frac{1}{2{\sigma }_Y^2}}{(y_i-{\tilde{x}}_i-\theta )}^2\right)\times \exp \left(-{\displaystyle \frac{1}{2{\sigma }_X^2}}{\tilde{x}}_i^2\right)$
		$\propto$	$\exp \left(-{\displaystyle \frac{1}{2{\sigma }_Y^2}}({\tilde{x}}_i^2-2(y_i-\theta ){\tilde{x}}_i+{(y_i-\theta )}^2)\right)\times \exp \left(-{\displaystyle \frac{1}{2{\sigma }_X^2}}{\tilde{x}}_i^2\right)$
		$\propto$	$\exp \left(-{\displaystyle \frac{1}{2}}\left\{\left({\displaystyle \frac{1}{{\sigma }_Y^2}}+{\displaystyle \frac{1}{{\sigma }_X^2}}\right){\tilde{x}}_i^2-2{\tilde{x}}_i\left({\displaystyle \frac{y_i-\theta }{{\sigma }_Y^2}}\right)\right\}\right).$

This gives

$${\tilde{x}}_i|\theta ,y_i\sim N(\left\{\frac{y_i-\theta }{{\sigma }_Y^2}\right\}/\left\{\frac{1}{{\sigma }_Y^2}+\frac{1}{{\sigma }_X^2}\right\},1/\left\{\frac{1}{{\sigma }_Y^2}+\frac{1}{{\sigma }_X^2}\right\}).$$