require(DNAmixtures)
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<-DNAmixture(list(tab3h1.df, tab3h2.df),   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)

#MLE parameters
ml30$mle

#\mu= rho*eta
ml30$mle[[1]]*ml30$mle[[2]]

# Coefficient of variation \sigma
1/sqrt(ml30$mle[[1]])




######## Prosecution hypothesis:   donor 1 and 1 unknown

mix30P<-DNAmixture(list(tab3h1.df, tab3h2.df),   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)

ml30P$mle


#\mu= rho*eta
ml30P$mle[[1]]*ml30P$mle[[2]]


# Coefficient of variation \sigma
1/sqrt(ml30P$mle[[1]])


 
####### Prosecution hypothesis:   donor 2 and 1 unknown##########

mix30P2<-DNAmixture(list(tab3h1.df, tab3h2.df),   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)

ml30P2$mle

#\mu= rho*eta
ml30P2$mle[[1]]*ml30P2$mle[[2]]

# Coefficient of variation \sigma
1/sqrt(ml30P2$mle[[1]])


#######Prosecution hypothesis:   donor 1 and donor 2 ##########

mix30P12<-DNAmixture(list(tab3h1.df, tab3h2.df),    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)

ml30P12$mle


#\mu= rho*eta
ml30P12$mle[[1]]*ml30P12$mle[[2]]


# Coefficient of variation \sigma
1/sqrt(ml30P12$mle[[1]])




###########Likelihood ratios log10(LR)#######

###WoE for H_p1: donor 1 & U1 vs H_d: U1 &U2###

ml30P$lik/log(10) - ml30$lik/log(10)


###WoE for H_p1: donor 1 & U1 vs H_d: U1 &U2###

ml30P$lik/log(10) - ml30$lik/log(10)


###WoE for H_p2: donor 2 & U1 vs H_d: U1 &U2###

ml30P2$lik/log(10) - ml30$lik/log(10)


###WoE for H_p12:  donor 1 & donor 2  vs H_d: U1 & U2###

ml30P12$lik/log(10) - ml30$lik/log(10)