A hierarchical t-distribution with unknown degrees of freedom

It is possible to treat the degrees of freedom parameter v as an additional unknown random variable in the model. This is most easily accomplished by specifying v to be a discrete variable. For example, let v take values 2,4,6,8,10,12,15,20,30 or 50. We then assume a uniform prior over the possible categories. That is, each value of v is given equal prior probability = 1/number of categories. The BUGS code is given in blockht.bug, and is shown on the next page.

Analysis

2000 iterations took 67 seconds, after a burn-in of 2000 iterations (note that we have used a longer burn-in than usual, since the degrees of freedom parameter typically converges quite slowly). This produced the following estimates: d = -.256 (SD .062); = -.260 (SD .149); = .095 (SD .054), and v = 10.300 (SD 4.746). A t distribution on v=10.3 degrees of freedom has a slightly heavier tail than Normal. However, the estimates of d and are very similar to those obtained by assuming a Normal population prior for the true treatment effects in each study, suggesting that allowing the different shaped distribution has little influence.

The number of categories and choice of values for v is somewhat arbitary. This can influence the resulting estimates since there is generally little information in the data concerning the value of v. We have found that a prior which places greater weight on low degrees of freedom, but also includes values large enough to give an approximately Normal distribution (e.g. v= 30, 50) works best. Some fine-tuning may be necessary to ensure that `jumps' between successive values of v are small enough to ensure that the sampler does not get `stuck' on a single value for hundreds of iterations. In addition, it is not necessary to assume equal prior probabilities for each category of v (see verdinelli:wasserman:91 for example). Note that assuming a continuous distribution for v would lead to a non log-concave full conditional distribution which BUGS is currently unable to sample from -- see the Dugongs example for further details of this type of problem.

Model specification for blockht example

model blockht;

const
   Num=22,    # Number of studies
   Nbins=10;   # Number of categories for v

var
  rt[Num], nt[Num], rc[Num], nc[Num], pc[Num], 
  pt[Num], mu[Num], delta[Num], d, tau, sigma, 
  delta.new, v, eta[Nbins], k, prior[Nbins];

data rt, nt, rc, nc in "blocker.dat";
inits in "blocker.in";

{

  for (i in 1:Num) {
       rt[i]         ~ dbin(pt[i],nt[i]);
       rc[i]         ~ dbin(pc[i],nc[i]);
       logit(pc[i]) <- mu[i];
       logit(pt[i]) <- mu[i] + delta[i];
       delta[i]      ~ dt(d,tau,v);
       mu[i]         ~ dnorm(0.0, 1.0E-5);
  }
  delta.new      ~ dt(d,tau,v);
  d      ~ dnorm(0.0,1.0E-6);
  tau    ~ dgamma(1.0E-3,1.0E-3);
  sigma <- 1/sqrt(tau);

  for (n in 1:Nbins) {
       prior[n] <- 1/Nbins;   # Uniform prior on v
  }
  k ~ dcat(prior[]);               
  v <- eta[k];           # degrees of freedom for t

  # Specify values taken by v: note that these
  # could be included in the data file instead

  eta[1] <- 2;   eta[2] <- 4;   eta[3] <- 6;
  eta[4] <- 8;   eta[5] <- 10   eta[6] <- 12;
  eta[7] <- 15;  eta[8] <- 20;  eta[9] <- 30;
  eta[10] <- 50;

}