---
title: "Population specific betas"
author: "Jérôme Goudet & Bruce Weir"
date: '`r Sys.Date()`'
output: html_document
---

Now some more about population specific betas, and their interpretation.  We will look at slides 176-181 in some details, using equations from slides 150-154, but without the sample size correction since we have allele frequencies.

```{r}
np<-4
#Bruce's example (slides 176-177)
paf1<-c(0.5,0.5)
paf2<-c(0.3,0.7)
pin1<-c(0.28,0.72)
pin2<-c(0.32,0.68)
mp<-cbind(paf1,paf2,pin1,pin2) # all pop all freq
hs<-1-t(mp)%*%mp #all diversities
(hb<-1/np/(np-1)*sum(hs-diag(4)*diag(hs))) # world Hb
(hba<-1/2*sum(hs[1:2,1:2]-diag(2)*hs[1:2,1:2])) # afr Hb
(hbi<-1/2*sum(hs[3:4,3:4]-diag(2)*hs[3:4,3:4])) # In Hb
(1-diag(hs)/hb) # betai world
(1-hs[1,1]/hba) # beta af1 rel africa
(1-hs[2,2]/hba) # beta af2 rel africa
(1-hs[3,3]/hbi) # beta in1 rel in
(1-hs[4,4]/hbi) # beta in2 rel in
```

Now we can write a function to carry out the calculations, basically a wrapper around the code of teh previous section:

```{r}
betas.freq<-function(freq){
np<-dim(freq)[2]
hs<-1-t(freq)%*%freq #all diversities
hb<-1/np/(np-1)*sum(hs-diag(np)*diag(hs)) # world Hb
betas<-1-diag(hs)/hb
betaw<-mean(betas)
return(list(betaw=betaw,betai=betas,hs=hs,hb=hb))
}
betas.freq(mp) #world
betas.freq(mp[,1:2]) #africa only
betas.freq(mp[,3:4]) #Inuit only
```

Now, imagine that we are dealing with microsats (STRs) rather than snps.  Each locus can have many alleles.  For instance, 4 alleles.  We are thinking of the first 2 populations being from Africa,  with a lot of variation, while the other two are for instance inuit, and have been populated from a sample of the african variation.  For instance:

```{r}
paf1<-c(0.23,0.22,0.25,0.30)
paf2<-c(0.24,0.23,0.27,0.26)
pin1<-c(0.82,0.18,0,0)
pin2<-c(0.87,0.13,0,0)

mp<-cbind(paf1,paf2,pin1,pin2) # all pop all freq
```

We can now estimate beta worlwide, or by continent:

```{r}
betas.freq(mp) #world
betas.freq(mp[,1:2]) #africa only
betas.freq(mp[,3:4]) #Inuit only
```