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') ##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), M=100000, 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) mlD$mle #H_P:suspect mixP<-DNAmixtureUAF(list(check1.df, check2.df), M=100000, 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) mlP$mle #WoE (mlP$lik/log(10))-(mlD$lik/log(10)) # 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) mlD$mle #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) mlP$mle #WoE (mlP$lik/log(10))-(mlD$lik/log(10))