> require(DNAmixtures) > data(USCaucasian,SGMplusDyes) > USCaucasian2<-rbind(USCaucasian,data.frame(marker="D18S51",allele=9,frequency=0.001659950)) > ff<-c('C:/Users/PeterGreen/Documents/Research/Perturb/Sep2013', + 'C:/Documents and Settings/mortera/Documenti/Peter-Julia/mixtures/UAF') > for(f in ff) if(file.exists(f)) {setwd(f); break} > source('DNAmixtureUAF.R') > > > tab3.df <-read.csv("Tab3_h1.txt", header=T, fill = TRUE, sep = "") > ##Starting up parameters > > startpar <- mixpar( rho = list(15), eta= list(60), xi = list(0.07), + phi= list(c(U1 = 0.8, U2 = 0.2))) > > ##Lower detection thereshold C=30 > mix30UAF<-DNAmixtureUAF(list(tab3.df), M=49, k=2, K = character(0), C=list(30), database = USCaucasian2, dyes = list(SGMplusDyes)) > > ml30UAF <- mixML(mix30UAF, startpar, trace = F) > > #likelihood > ml30UAF$lik/log(10) [1] -122.7308 > > #MLE parameters > ml30UAF$mle rho eta xi phi.U1 phi.U2 1 13.7 65.83 0.07198 0.9167 0.08334 > > #\mu= rho*eta > ml30UAF$mle[[1]]*ml30UAF$mle[[2]] [1] 901.8654 > # 901.9651 > > # Coefficient of variation \sigma > 1/sqrt(ml30UAF$mle[[1]]) [1] 0.2701663 > #[1] 0.2709708 > > #mixture proportion for two unknowns > ml30UAF$mle[[4]] U1 U2 0.91665723 0.08334277 > # U1 U2 > # 0.91849985 0.08150015 > > ######## Prosecution hypothesis: donor 1 and 1 unknown > > mix30UAFP<-DNAmixtureUAF(list(tab3.df), M=49, k=2, K = "donor1", C=list(30), database = USCaucasian2, dyes = list(SGMplusDyes)) > > par(mfrow = c(3, 4)) > plot(mix30UAFP, threshold = TRUE) > > startpar <- mixpar( rho = list(15), eta= list(60), xi = list(0.07), phi= list(c(donor1 = 0.8, U1 = 0.2))) > > ml30UAFP <- mixML(mix30UAFP, startpar, trace = F) > > ml30UAFP$lik/log(10) [1] -109.5012 > > ml30UAFP$mle rho eta xi phi.U1 phi.donor1 1 13.63 66.18 0.07168 0.08106 0.9189 > > > #\mu= rho*eta > ml30UAFP$mle[[1]]*ml30UAFP$mle[[2]] [1] 901.9493 > > > # Coefficient of variation \sigma > 1/sqrt(ml30UAFP$mle[[1]]) [1] 0.270884 > > > #mixture proportion for H_P: donor 1 and 1 unknown U1 > ml30UAFP$mle[[4]] U1 donor1 0.08106162 0.91893838 > # U1 donor 1 > > > ####### Prosecution hypothesis: donor 2 and 1 unknown########## > > mix30UAFP2<-DNAmixtureUAF(list(tab3.df), M=49, k=2, K = "donor2", C=list(30), database = USCaucasian2, + dyes = list(SGMplusDyes)) > > startpar <- mixpar( rho = list(15), eta= list(60), xi = list(0.07), phi= list(c(donor2 = 0.2, U1 = 0.8))) > > ml30UAFP2 <- mixML(mix30UAFP2, startpar, trace = F) > > ml30UAFP2$lik/log(10) [1] -115.2245 > > ml30UAFP2$mle rho eta xi phi.U1 phi.donor2 1 16.37 55.08 0.07011 0.8839 0.1161 > > #\mu= rho*eta > ml30UAFP2$mle[[1]]*ml30UAFP2$mle[[2]] [1] 901.5648 > > # Coefficient of variation \sigma > 1/sqrt(ml30UAFP2$mle[[1]]) [1] 0.2471769 > > #mixture proportion for H_P: donor 2 and 1 unknown U1 > ml30UAFP2$mle[[4]] U1 donor2 0.8839257 0.1160743 > > #######Prosecution hypothesis: donor 1 and donor 2 ########## > > mix30UAFP12<-DNAmixtureUAF(list(tab3.df), M=49, k=2, K = c("donor1", "donor2"), C=list(30), database = USCaucasian2, + dyes = list(SGMplusDyes)) > > #par(mfrow = c(3, 4)) > #plot(mix30UAFP12, threshold = TRUE) > > startpar <- mixpar( rho = list(15), eta= list(60), xi = list(0.07), + phi= list(c(donor1 = 0.8, donor2 = 0.2))) > > ml30UAFP12 <- mixML(mix30UAFP12, startpar, trace = F + #, phi.eq = FALSE + ) > > ml30UAFP12$lik/log(10) [1] -102.0184 > > ml30UAFP12$mle rho eta xi phi.donor1 phi.donor2 1 16.37 55.08 0.07011 0.8839 0.1161 > > > #\mu= rho*eta > ml30UAFP12$mle[[1]]*ml30UAFP12$mle[[2]] [1] 901.565 > > > # Coefficient of variation \sigma > 1/sqrt(ml30UAFP12$mle[[1]]) [1] 0.2471767 > > > > #mixture proportion for H_P: donor 1 and donor 2 > ml30UAFP12$mle[[4]] donor1 donor2 0.8839258 0.1160742 > > > ###########Likelihood ratios log10(LR)####### > > ###WoE for H_p1: donor 1 & U1 vs H_d: U1 &U2### > > ml30UAFP$lik/log(10) - ml30UAF$lik/log(10) [1] 13.22953 > > > ###COMPARISON with Puch-Solis et al paper Example 5 (PS): > ###WoE for H_p1: donor 1 & U1 vs H_d: U1 &U2### > > 10^(ml30UAFP$lik/log(10) - ml30UAF$lik/log(10)) [1] 1.69639e+13 > # DNA mixtures UAF > > # whereas PS give 1.7e +9 > > ###COMPARISON with Puch-Solis et al paper Example 5 (PS): > ###WoE for H_p2: donor 2 & U1 vs H_d: U1 &U2### > > ml30UAFP2$lik/log(10) - ml30UAF$lik/log(10) [1] 7.506286 > > 10^(ml30UAFP2$lik/log(10) - ml30UAF$lik/log(10)) [1] 32083803 > > # DNA mixtures UAF > > # whereas PS give 4.8e +6 > > ####WoE for H_p12: donor 1 & donor 2 vs H_p1: donor 1 &U1### > > ml30UAFP12$lik/log(10) - ml30UAFP$lik/log(10) [1] 7.482813 > > ###WoE for H_p12: donor 1 & donor 2 vs H_d: U1 & U2### > > ml30UAFP12$lik/log(10) - ml30UAF$lik/log(10) [1] 20.71234 >