A Bayesian approach

We saw last time how to use WinBUGS to provide allele frequency estimates from phenotypic data at the ABO locus.

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))

As you may recall, this produced the results in Figure 1.

Figure 1: Results from WinBUGS analysis of the ABO data assuming genotypes are in Hardy-Weinberg proportions.

Now that we know about inbreeding coefficients and that the allow us to measure the departure of genotype frequencies from Hardy-Weinberg proportions, we can modify this a bit and estimate allele frequencies without assuming that genotypes are in Hardy-Weinberg proportions.

model {
   # likelihood 
   pi[1] <- p.a*p.a + f*p.a*(1-p.a) + 2*p.a*p.o*(1-f)
   pi[2] <- 2*p.a*p.b*(1-f)
   pi[3] <- p.b*p.b + f*p.b*(1-p.b) + 2*p.b*p.o*(1-f)
   pi[4] <- p.o*p.o + f*p.o*(1-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)

   f ~ dunif(0,1)

   n <- sum(x[])
}

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

This produces the results in Figure 2

Figure 2: Results from WinBUGS analysis of the ABO data relaxing the assumption that genotypes are in Hardy-Weinberg proportions.

Notice that the allele frequency estimates have changed quite a bit and that the posterior mean of $f$ is about 0.41. On first appearance, that would seem to indicate that we have lots of inbreeding in this sample. BUT it's a human population. That doesn't seem very likely. Take a closer look. The 95% credible interval for $f$ is between 0.06 and 0.55. That suggests that we don't have much information at all about $f$ from these data.³ How can we tell if the model with inbreeding is better than the model that assumes genotypes are in Hardy-Weinberg proportions?