> 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') > > ##data for two traces > check1.df <-read.csv("check.txt", header=T, fill = TRUE, sep = "") > check2.df <-read.csv("check190.txt", header=T, fill = TRUE, sep = "") > #Set up a model for analysis of multiple profiles h1 and h2 of Table 3 Puch-Solis with unknown > ########## > ##Assume that the contributor to the two traces is the same## > > # UAF version > #H_D: U1 > > ##Lower detection thereshold C=0.001 > mixD<-DNAmixtureUAF(list(check1.df, check2.df), k=1, K = character(0), C=list(0.001,0.001), database = USCaucasian2, dyes = list(SGMplusDyes, SGMplusDyes)) > > > ##Starting up parameters > > startpar <- mixpar( rho = list(2,5), eta= list(200,200), xi = list(0.1, 0.1), + phi= list(c(U1 = 0.9), c(U1 = 0.9))) > > ## Here we set stutter equal for both traces > #restricting xis to be equal for the two profiles > > eqxis <- function(x){ diff(unlist(x[,"xi"])) } > ## Note that for these two traces, we do not expect phi to be equal. > > > #now maximizing the likelihood > > > mlD <- mixML(mixD, startpar, eqxis, val = 0, trace = F, phi.eq = FALSE) > > mlD$lik/log(10) [1] -11.18351 > mlD$mle rho eta xi phi.U1 1 2.033 121.02 0.000000000000000000004063 1 2 2.714 37.52 0.000000000000000000004063 1 > > #H_P:suspect > mixP<-DNAmixtureUAF(list(check1.df, check2.df), k=1, K = "suspect", C=list(0.001,0.001), database = USCaucasian2, dyes = list(SGMplusDyes, SGMplusDyes)) > > ##Starting up parameters > > startpar <- mixpar( rho = list(2,5), eta= list(200,200), xi = list(0.1, 0.1), + phi= list(c(suspect = 0.9), c(suspect = 0.9))) > > ## Here we set stutter equal for both traces > #restricting xis to be equal for the two profiles > > eqxis <- function(x){ diff(unlist(x[,"xi"])) } > ## Note that for these two traces, we do not expect phi to be equal. > > > mlP<-mixML(mixP, startpar, eqxis, val = 0, trace = F, phi.eq = FALSE) > > mlP$lik/log(10) [1] -8.421626 > mlP$mle rho eta xi phi.suspect 1 2.484 100.64 0.00000000000000000001943 1 2 90.370 1.05 0.00000000000000000001943 1 > > #WoE > (mlP$lik/log(10))-(mlD$lik/log(10)) [1] 2.761884 > > # vanilla version > #H_D: U1 > > > ##Lower detection thereshold C=0.001 > mixD<-DNAmixture(list(check1.df, check2.df), k=1, K = character(0), C=list(0.001,0.001), database = USCaucasian2, dyes = list(SGMplusDyes, SGMplusDyes)) > > > ##Starting up parameters > > startpar <- mixpar( rho = list(2,5), eta= list(200,200), xi = list(0.1, 0.1), + phi= list(c(U1 = 0.9), c(U1 = 0.9))) > > ## Here we set stutter equal for both traces > #restricting xis to be equal for the two profiles > > eqxis <- function(x){ diff(unlist(x[,"xi"])) } > ## Note that for these two traces, we do not expect phi to be equal. > > > #now maximizing the likelihood > > > mlD <- mixML(mixD, startpar, eqxis, val = 0, trace = F, phi.eq = FALSE) > > mlD$lik/log(10) [1] -10.04862 > mlD$mle rho eta xi phi.U1 1 2.563 100.030 0.000000000000000000005903 1 2 37.269 2.561 0.000000000000000000005903 1 > > #H_P:suspect > mixP<-DNAmixture(list(check1.df, check2.df), k=1, K = "suspect", C=list(0.001,0.001), database = USCaucasian2, dyes = list(SGMplusDyes, SGMplusDyes)) > > ##Starting up parameters > > startpar <- mixpar( rho = list(2,5), eta= list(200,200), xi = list(0.1, 0.1), + phi= list(c(suspect = 0.9), c(suspect = 0.9))) > > ## Here we set stutter equal for both traces > #restricting xis to be equal for the two profiles > > eqxis <- function(x){ diff(unlist(x[,"xi"])) } > ## Note that for these two traces, we do not expect phi to be equal. > > > mlP<-mixML(mixP, startpar, eqxis, val = 0, trace = F, phi.eq = FALSE) > > mlP$lik/log(10) [1] -8.421626 > mlP$mle rho eta xi phi.suspect 1 2.484 100.64 0.00000000000000000001943 1 2 90.370 1.05 0.00000000000000000001943 1 > > #WoE > (mlP$lik/log(10))-(mlD$lik/log(10)) [1] 1.626997 >