Hearts: a mixture model for count data

The table below presents data given by berry:87 on the effect of a drug used to treat patients with frequent premature ventricular contractions (PVCs) of the heart.

farewell:sprott:88 model this data as a mixture distribution of Poisson counts in which some patients are ``cured'' by the drug, whilst others experience varying levels of response but remain abnormal. A zero count for the post-drug PVC may indicate a ``cure'', or may represent a sampling zero from a patient with a mildly abnormal PVC count. The following model thus is assumed:

To eliminate the nuisance parameters , Farewell and Sprott use the conditional distribution of given . This is equivalent to a binomial likelihood for with denominator and probability (see cox:hinkley:74 pp. 136-137 for further details of the conditional distribution for Poisson variables). Hence the final mixture model may be expressed as follows:

The graph for the hearts model is shown in Figure 4 and the BUGS code follows.

 
Figure 4:   Graphical model for hearts example

model hearts;
const N = 12;
var 
  x[N],y[N],t[N],     # pre-drug, post-drug and total PVC count
  state[N],state1[N], # binary indicator of whether patient is cured
  theta,              # probability of cure (prob of state = 1)
  p, beta,            # p = binomial probability = beta/(1+beta)
  P[2],               # `pick' variable used to select appropriate
                      # value for the binomial probability depending
                      # on whether state1 = 1 or 2 (not cured or cured)
  alpha, delta;       # p and theta transformed to logit scale for normality

data x, y in "hearts.dat";
inits in "hearts.in";

{
# TRANSFORMATIONS
   for (i in 1:N) {
      t[i] <- x[i] + y[i];
   }
# MODEL
   for (i in 1:N) {
      y[i] ~ dbin(P[state1[i]], t[i]);
      state[i] ~ dbern(theta);
      state1[i] <- state[i]+1;   # state[i] takes values 0 or 1, so need to
                                 # add 1 to get values for use as index on P
   }
   P[1] <- p; P[2] <- 0;
   logit(p) <- alpha; alpha ~ dnorm(0,1.0E-4); 
   beta <- exp(alpha);  # beta measures change in rate of PVCs after treatment
   logit(theta) <- delta; delta ~ dnorm(0,1.0E-4)
}

Analysis

10000 iterations took 32 seconds after a 1000 iteration burn-in. The posterior means (95% C.I.) are given below, together with Farewell and Sprott's maximum likelihood estimates.

The BUGS results are in close agreement with those of Farewell and Sprott, and suggest that there is just over a 50% chance of a patient being ``cured'' ( = 0.572); if not, the number of PVCs per minute is likely to fall to about 65% ( = 0.646) of their pre-drug count.