Returning to the ABO example

Here's data from the ABO blood group:²¹

Phenotype	A	AB	B	O	Total
Observed	862	131	365	702	2060

To estimate the underlying allele frequencies, $p_A$, $p_B$, and $p_O$, we have to remember how the allele frequencies map to phenotype frequencies:²²

$$\begin{eqnarray*} \hbox{Freq}(A) &=& p_A^2 + 2p_Ap_O \\ \hbox{Freq}{AB} &=& 2p_... ...req}{B} &=& p_B^2 + 2p_Bp_O \\ \hbox{Freq}{O} &=& p_O^2 \quad . \end{eqnarray*}$$

Hers's the WinBUGS code we use to estimate the allele frequencies:

model {
   # likelihood 
   pi[1] <- p.a*p.a + 2*p.a*p.o
   pi[2] <- 2*p.a*p.b
   pi[3] <- p.b*p.b + 2*p.b*p.o
   pi[4] <- p.o*p.o
   x[1:4] ~ dmulti(pi[],n)

   # priors
   a1 ~ dexp(1)
   b1 ~ dexp(1)
   o1 ~ dexp(1)
   p.a <- a1/(a1 + b1 + o1)
   p.b <- b1/(a1 + b1 + o1)
   p.o <- o1/(a1 + b1 + o1)

   n <- sum(x[])
}

list(x=c(862, 131, 365, 702))

The dmulti() is a multinomial probability, a simple generalization of the binomial probability to samples when there are more than two categories. The priors are some mumbo jumbo necessary to produce the rough equivalent of uniform [0,1] priors with more than two alleles.²³ sum() is a built-in function that saves the trouble of calculating the sample size, and ensures that the n in dmulti() is consistent with the individual sample components. The x=c() produces a vector of counts arranged in the same order as the frequencies in pi[]. The results are in Figure 2 Notice that the posterior means are very close to the maximum-likelihood estimates, but that we also have 95% credible intervals so that we have an assessment of how reliable the Bayesian estimates are. Getting them from a likelihood analysis is possible, but it takes a fair amount of additional work.

Figure 2: Results of a WinBUGS analysis of the ABO data.