#Example 1: Motivating example Egeland et al. (2013)
require(paramlink)
y1=swapSex(nuclearPed(3),c(3,4))
p=c(0.1,0.2,0.3,0.4)
alleles=1:length(p)
T1=c(1,1)
T2=c(2,2)
R=1:2
known=list(c(3,T1),c(4,T2))
l1=paraMix(y1,R,id.U=5,alleles=alleles,afreq=p,known_genotypes=known)
y2=swapSex(nuclearPed(1),3)
y2=addOffspring(y2,mother=2,noff=1,sex=2)
y2=relabel(y2,c(1:3,6,4),1:5)
l2=paraMix(y2,R,id.U=6,alleles=alleles,afreq=p,known_genotypes=known)
LR1=l1$lik/l2$lik
exact=1/(2*(p[1]+p[2]))
stopifnot(abs(LR1-exact)<10^(-6))
#Example 2. Example 1 in Egeland et al. (2013) based on Fung and Hu (2008)
#Data:
#Mixture 1/2/3
#Suspect=4, genotype 3/3
#Victim=10, genotype 1/2
#H1: Contributors were the suspect and victim (unrelated)
#H2: Contributors were the father of suspect and victim (unrelated)
#H3: Contributors were the brother of suspect and victim (unrelated)
afreq=c(0.044,0.166,0.110,0.680)
alleles=1:length(afreq)
R=1:3 #Mixture
man_ped=nuclearPed(2)
victim = singleton(id=10, sex=2)
known = list(c(4,3,3),c(10,1,2)) #individual 4 is 3/3, and 10 (the victim) is 1/2.
#The likelihoods corresponding to H1,H2 and H3
l1=paraMix(list(man_ped, victim), R, id.U=NULL, id.V=NULL,
alleles=alleles, afreq=afreq, known_genotypes=known)$lik
l2=paraMix(list(man_ped, victim), R, id.U=1, id.V=4,
alleles=alleles, afreq=afreq, known_genotypes=known)$lik
l3=paraMix(list(man_ped, victim), R, id.U=3, id.V=4,
alleles=alleles, afreq=afreq, known_genotypes=known)$lik
LR12=l1/l2
stopifnot(abs(LR12-3.125)<10^(-6))
LR13=l1/l3
stopifnot(abs(LR13- 2.355296)<10^(-6))
Run the code above in your browser using DataLab