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') source('markerwiseLikelihood.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))) ##Defence hypothesis two unknownsH_d:U1 &U2### ################################################ ##detection thereshold C=30## ##UAF## ################################ 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) ###overall likelihood## ml30UAF$lik/log(10) ###################### #markerwise analysis## ###################### logL(mix30UAF)(ml30UAF$mle) source('markerwiselikelihood.R') logL_m(mix30UAF)(ml30UAF$mle) ###################### # NO UAF ## ###################### mix30<-DNAmixture(list(tab3.df), k=2, K = character(0), C=list(30), database = USCaucasian2, dyes = list(SGMplusDyes)) ml30 <- mixML(mix30, startpar, trace = F) #likelihood ml30$lik/log(10) ###################### #markerwise analysis## ###################### logL(mix30)(ml30$mle) source('markerwiselikelihood.R') logL_m(mix30)(ml30$mle) ######## Prosecution hypothesis: donor 1 and 1 unknown ##UAF## mix30UAFP<-DNAmixtureUAF(list(tab3.df), k=2, K = "donor1", C=list(30), database = USCaucasian2, dyes = list(SGMplusDyes)) 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) ###################### #markerwise analysis## ###################### logL(mix30UAFP)(ml30UAFP$mle) source('markerwiselikelihood.R') logL_m(mix30UAFP)(ml30UAFP$mle) ######## Prosecution hypothesis: donor 1 and 1 unknown ## NO UAF## mix30P<-DNAmixture(list(tab3.df), k=2, K = "donor1", C=list(30), database = USCaucasian2, dyes = list(SGMplusDyes)) startpar <- mixpar( rho = list(15), eta= list(60), xi = list(0.07), phi= list(c(donor1 = 0.8, U1 = 0.2))) ml30P <- mixML(mix30P, startpar, trace = F) ml30P$lik/log(10) ###################### #markerwise analysis## ###################### logL(mix30P)(ml30P$mle) source('markerwiselikelihood.R') logL_m(mix30P)(ml30P$mle) ###########Likelihood ratios log10(LR)####### ###WoE for H_p1: donor 1 & U1 vs H_d: U1 &U2### ml30UAFP$lik/log(10) - ml30UAF$lik/log(10) ##LRmUAF## LRmUAF=logL_m(mix30UAFP)(ml30UAFP$mle)/log(10)- logL_m(mix30UAF)(ml30UAF$mle)/log(10) LRmUAF ##LRm## LRm=logL_m(mix30P)(ml30P$mle)/log(10)- logL_m(mix30)(ml30$mle)/log(10) LRm LRmUAF/LRm ###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)) # DNA mixtures 1.488128e+13 whereas PS give 1.7e +9 ####### Prosecution hypothesis: donor 2 and 1 unknown########## ###UAF###### mix30UAFP2<-DNAmixtureUAF(list(tab3.df), 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) ###################### #markerwise analysis## ###################### logL(mix30UAFP2)(ml30UAFP2$mle) source('markerwiselikelihood.R') logL_m(mix30UAFP2)(ml30UAFP2$mle) ####### Prosecution hypothesis: donor 2 and 1 unknown########## ### NO UAF###### mix30P2<-DNAmixture(list(tab3.df), 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))) ml30P2 <- mixML(mix30P2, startpar, trace = F) ml30P2$lik/log(10) ###################### #markerwise analysis## ###################### logL(mix30P2)(ml30P2$mle) source('markerwiselikelihood.R') logL_m(mix30P2)(ml30P2$mle) ###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.470102 ##LRmUAF## LR2mUAF=logL_m(mix30UAFP2)(ml30UAFP2$mle)/log(10)- logL_m(mix30UAF)(ml30UAF$mle)/log(10) LR2mUAF ##LR2m## LR2m=logL_m(mix30P2)(ml30P2$mle)/log(10)- logL_m(mix30)(ml30$mle)/log(10) LR2m LR2mUAF/LR2m # DNA mixtures 29.52e +6 whereas PS give 4.8e +6