> data(USCaucasian,SGMplusDyes)
> USCaucasian2<-rbind(USCaucasian,data.frame(marker="D18S51",allele=9,frequency=0.001659950))
> ff<-c(paste(Sys.getenv('HOME'),'/Research/Perturb/Sep2013',sep=''),
+     paste(Sys.getenv('HOME'),'/Peter-Julia/mixtures/UAF',sep=''))
> for(f in ff) if(file.exists(f)) {setwd(f); break}
> source('DNAmixtureUAF.R')
> 
> 
> tab3h1.df <-read.csv("Tab3_h1.txt", header=T, fill = TRUE, sep = "")
> tab3h2.df <-read.csv("Tab3_h2.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 contributors assumed identical
> 
> #Under the defence Hypothesis U1&U2
> 
> ##Lower detection thereshold C=30
> mix30<-DNAmixtureUAF(list(tab3h1.df, tab3h2.df), M=49, k=2,  K = character(0), C=list(30,30), database = USCaucasian2, dyes = list(SGMplusDyes, SGMplusDyes))
> 
> 
> ##Starting up parameters
> 
> startpar <- mixpar( rho = list(15,15), eta= list(60,60), xi = list(0.1,0.1),
+         phi= list(c(U1 = 0.8, U2 = 0.2), c(U1 = 0.6, U2 = 0.4)))
> 
> 
> 
> ## 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.
> 
> 
>  ## Note that for these two traces, we do not expect phi to be equal.
> 
> # eqparams <- function(x){ c(diff(unlist(x[,"xi"])), diff(unlist(x[,"eta"])))}
> 
> #now maximizing the likelihood
> 
> 
> ml30 <- mixML(mix30, startpar,  eqxis, val=0, trace = F, phi.eq = FALSE)
> 
> #likelihood
> ml30$lik/log(10)
[1] -236.8617
> 
> #MLE parameters
> ml30$mle
      rho     eta      xi   phi.U1   phi.U2
1   16.69   54.01   0.061   0.8836   0.1164
2   26.09   43.41   0.061   0.5446   0.4554
> 
> #\mu= rho*eta
> ml30$mle[[1]]*ml30$mle[[2]]
[1] 435.6157
> 
> # Coefficient of variation \sigma
> 1/sqrt(ml30$mle[[1]])
[1] 0.2447521
> 
> 
> 
> 
> ######## Prosecution hypothesis:   donor 1 and 1 unknown
> 
> mix30P<-DNAmixtureUAF(list(tab3h1.df, tab3h2.df), M=49, k=2,  K = "donor1", C=list(30,30), database = USCaucasian2, dyes = list(SGMplusDyes, SGMplusDyes))
> 
> 
> ##Starting up parameters
> 
> startpar <- mixpar( rho = list(15,15), eta= list(60,60), xi = list(0.1,0.1),
+         phi= list(c(donor1 = 0.8, U1 = 0.2), c(donor1 = 0.6, U1 = 0.4)))
> 
> 
> 
> ## 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.
> 
> # eqparams <- function(x){ c(diff(unlist(x[,"xi"])), diff(unlist(x[,"eta"])))}
> 
> #now maximizing the likelihood
> 
> 
> #par(mfrow = c(3, 4))
> #plot(mix30P, threshold = TRUE)
> 
> 
> ml30P <- mixML(mix30P, startpar,  eqxis, val=0, trace = F, phi.eq = FALSE)
> 
> 
> ml30P$lik/log(10)
[1] -223.756
> 
> ml30P$mle
      rho     eta      xi   phi.U1   phi.donor1
1   16.76   53.80   0.061   0.1164       0.8836
2   26.07   43.46   0.061   0.4554       0.5446
> 
> 
> #\mu= rho*eta
> ml30P$mle[[1]]*ml30P$mle[[2]]
[1] 436.7349
> 
> 
> # Coefficient of variation \sigma
> 1/sqrt(ml30P$mle[[1]])
[1] 0.2443005
> 
> 
>  
> ####### Prosecution hypothesis:   donor 2 and 1 unknown##########
> 
> mix30P2<-DNAmixtureUAF(list(tab3h1.df, tab3h2.df), M=49, k=2,  K = "donor2", C=list(30,30), database = USCaucasian2, dyes = list(SGMplusDyes, SGMplusDyes))
> 
> 
> ##Starting up parameters
> 
> startpar <- mixpar( rho = list(15,15), eta= list(60,60), xi = list(0.1,0.1),
+         phi= list(c(donor2 = 0.8, U1 = 0.2), c(donor2 = 0.6, U1 = 0.4)))
> 
> 
> 
> ## 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.
> 
> # eqparams <- function(x){ c(diff(unlist(x[,"xi"])), diff(unlist(x[,"eta"])))}
> 
> #now maximizing the likelihood
> 
> 
> 
> ml30P2 <- mixML(mix30P2, startpar,  eqxis, val=0, trace = F, phi.eq = FALSE)
> 
> 
> ml30P2$lik/log(10)
[1] -224.4913
> 
> ml30P2$mle
      rho     eta        xi   phi.U1   phi.donor2
1   16.77   53.74   0.06095   0.8835       0.1165
2   26.18   43.27   0.06095   0.5441       0.4559
> 
> #\mu= rho*eta
> ml30P2$mle[[1]]*ml30P2$mle[[2]]
[1] 439.1532
> 
> # Coefficient of variation \sigma
> 1/sqrt(ml30P2$mle[[1]])
[1] 0.2441612
> 
> 
> #######Prosecution hypothesis:   donor 1 and donor 2 ##########
> 
> mix30P12<-DNAmixtureUAF(list(tab3h1.df, tab3h2.df),  M=49, k=2,  K = c("donor1", "donor2"), C=list(30,30), database = USCaucasian2, dyes = list(SGMplusDyes, SGMplusDyes))
> 
> 
> #par(mfrow = c(3, 4))
> #plot(mix30P12, threshold = TRUE)
> 
> startpar <- mixpar( rho = list(15,15), eta= list(60,60), xi = list(0.1,0.1),
+         phi= list(c(donor1 = 0.8, donor2 = 0.2), c(donor1 = 0.6, donor2 = 0.4)))
> 
> 
> 
> ml30P12 <- mixML(mix30P12, startpar,  eqxis, val=0, trace = F, phi.eq = FALSE)
> 
> 
> 
> ml30P12$lik/log(10)
[1] -211.2854
> 
> ml30P12$mle
      rho    eta        xi   phi.donor1   phi.donor2
1   16.70   54.0   0.06095       0.8835       0.1165
2   26.22   43.2   0.06095       0.5441       0.4559
> 
> 
> #\mu= rho*eta
> ml30P12$mle[[1]]*ml30P12$mle[[2]]
[1] 437.7024
> 
> 
> # Coefficient of variation \sigma
> 1/sqrt(ml30P12$mle[[1]])
[1] 0.2447346
> 
> 
> 
> 
> ###########Likelihood ratios log10(LR)#######
> 
> ###WoE for H_p1: donor 1 & U1 vs H_d: U1 &U2###
> 
> ml30P$lik/log(10) - ml30$lik/log(10)
[1] 13.10579
> 
> 
> ###WoE for H_p1: donor 1 & U1 vs H_d: U1 &U2###
> 
> ml30P$lik/log(10) - ml30$lik/log(10)
[1] 13.10579
> 
> 
> ###WoE for H_p2: donor 2 & U1 vs H_d: U1 &U2###
> 
> ml30P2$lik/log(10) - ml30$lik/log(10)
[1] 12.37039
> 
> 
> ###WoE for H_p12:  donor 1 & donor 2  vs H_d: U1 & U2###
> 
> ml30P12$lik/log(10) - ml30$lik/log(10)
[1] 25.57633
>