options(java.parameters = "- Xmx1024m")
pkgTest <- function(x) {
  if (!require(x, character.only = TRUE))
  {
    install.packages(x, dep = TRUE)
    if(!require(x, character.only = TRUE)) stop("Package not found")
  }
}
packages <- c("readxl","data.table","stats","fitdistrplus","openxlsx")
lapply(packages, pkgTest)

grill.wd<-"Destination folder for new files" # set here the name of the directory where you would like to save output files
online.wd<-"Source folder for the Excel Data file" # set here the name of the working directory with the Excel file with data from the literature
setwd(grill.wd)

#################
## functions
###########

#################
## end functions
###########

current.wd<-getwd()
setwd(online.wd)
grill.lit.dt<-data.table(read_excel("Excel Data file name here!"))[,Excel.line:=.I+1] # set the name of the file with literature results here
setwd(current.wd)
## if alternative.inclusion==TRUE, log-fire grilling is excluded throughout the analysis and the files are generated accordingly
alternative.inclusion<-FALSE

names(grill.lit.dt)

chem.names<-gsub("^([^ ]+) .+$","\\1",grep("pooled mean",names(grill.lit.dt),value = TRUE))

for(cn in chem.names){
  for (cn.cn in grep(cn,names(grill.lit.dt),value = TRUE)){
    if (grepl("LOD",cn.cn)){
      setnames(grill.lit.dt,cn.cn, paste0(cn,".LOD"))
    }
    if (grepl("LOQ",cn.cn)){
      setnames(grill.lit.dt,cn.cn, paste0(cn,".LOQ"))
    }
    if (grepl("pooled mean$",cn.cn)){
      setnames(grill.lit.dt,cn.cn, paste0(cn,".Mean"))
    }
    if (grepl(" sd$",cn.cn)){
      setnames(grill.lit.dt,cn.cn, paste0(cn,".SD"))
    }
    if (grepl(" min$",cn.cn)){
      setnames(grill.lit.dt,cn.cn, paste0(cn,".min"))
    }
    if (grepl(" maks$",cn.cn)){
      setnames(grill.lit.dt,cn.cn, paste0(cn,".max"))
    }
  }
  cat(grep(cn,names(grill.lit.dt),value = TRUE),sep = "\t")
  cat("\n")
}

setnames(grill.lit.dt,"Grilling method","grilling.method")
included.methods<-c(
  "Briquettes",
  "Charcoal",
  "Charcoal + wood chips",
  "Disposable grill",
  "Electric",
  "Gas",
  "Grilled (unspecified method)",
  "Home or restaurant grilled",
  "Wood/log fire")

setdiff(grill.lit.dt[,grilling.method],included.methods)
setdiff(included.methods,grill.lit.dt[,grilling.method])

setnames(grill.lit.dt,"Food item new","food.name")

grill.lit.dt[food.name=="lamb",food.name:="Lamb"]
grill.lit.dt[food.name=="Other food items",food.name:="Other food"]

table(grill.lit.dt[grilling.method %in% included.methods,.(food.name,grilling.method)],useNA = "always")
table(grill.lit.dt[grilling.method %in% included.methods,.(food.name)],useNA = "always")

if (length(names(grill.lit.dt))!=length(unique(names(grill.lit.dt)))){
  stop("Duplicate names in the table\n")
}

setnames(grill.lit.dt,"n","n.obs")

grill.lit.dt[,
  c("Excel.line","n.obs","grilling.method","food.name",grep(paste("(",chem.names,")",collapse = "|",sep=""),names(grill.lit.dt),value = TRUE)),
  with = FALSE]

grill.dt<-rbindlist(lapply(chem.names, function(cn) {
    dt<-grill.lit.dt[,
      c("Excel.line","n.obs","grilling.method","food.name",grep(cn,names(grill.lit.dt),value = TRUE)),
      with = FALSE][,chem.name:=cn]
    setnames(dt,names(dt),gsub(paste0("^",cn,"[ .]"),"",names(dt)))
  }),use.names = TRUE,fill = TRUE)[!is.na(food.name)]

# rule for filling in unknown n for tests with multiple samples
grill.dt[is.na(n.obs) & !is.na(Mean) & !is.na(SD),n.obs:=3]
grill.dt[is.na(n.obs) & !is.na(nd),n.obs:=3]

if (nrow(grill.dt[!is.na(na) & (!is.na(LOQ)| !is.na(SD)| !is.na(min)| !is.na(max))])>0){
  fwrite(grill.dt[!is.na(na) & (!is.na(LOQ)| !is.na(SD)| !is.na(min)| !is.na(max))],"irregular.na.lines.txt",sep = "\t")
  stop("There are some not fixed na lines")
}

### remove lines without analysis
if (length(grep("\\.NA$",names(unlist(lapply(grill.dt[!is.na(na),.(LOQ,`pooled median`,Mean,nd,SD,min,max)],table,useNA="a"))),invert = TRUE))>0){
  stop("some unexpected values in non-analyzed lines\n")
}
grill.dt<-grill.dt[!(!is.na(na) & is.na(LOQ) & is.na(Mean) & is.na(min) & is.na(max)),]
if (nrow(grill.dt[!is.na(na)])>0){
  stop("some remaining non-analyzed lines\n")
}

## check for LOD/LOQ relative values
if (nrow(grill.dt[LOQ/LOD<2])>0){
  fwrite(grill.dt[LOQ/LOD<2],"relatively.small.LOQs.txt",sep="\t")
  warning("Still some wrong relations between LOQ and LOD")  
}

# LOQ that are coupled with too low means
grill.dt[Mean<LOQ/n.obs]
fwrite(grill.dt[Mean<LOQ/n.obs],"relatively.small.LOQs.v2.txt",sep="\t")

## check for LOQ outliers
grill.dt[,out.up:=LOQ>mean(LOQ[grilling.method %in% included.methods],na.rm=TRUE)*3,by=c("chem.name","food.name")]
grill.dt[,out.down:=LOQ<mean(LOQ[grilling.method %in% included.methods],na.rm=TRUE)/10,by=c("chem.name","food.name")]
fwrite(grill.dt[out.down==TRUE],"absolutely.small.LOQs.txt",sep="\t")
fwrite(grill.dt[out.up==TRUE],"absolutely.large.LOQs.txt",sep="\t")
grill.dt[,out.up:=LOD>mean(LOD[grilling.method %in% included.methods],na.rm=TRUE)*3,by=c("chem.name","food.name")]
grill.dt[,out.down:=LOD<mean(LOD[grilling.method %in% included.methods],na.rm=TRUE)/10,by=c("chem.name","food.name")]
fwrite(grill.dt[out.down==TRUE],"absolutely.small.LODs.txt",sep="\t")
fwrite(grill.dt[out.up==TRUE],"absolutely.large.LODs.txt",sep="\t")

# the above checks are used to clear up all errors that are possible to clear up in the literature data (update Excel input)

## after these are cleared up will fill in the missing values

table(grill.dt[is.na(as.numeric(LOD)),LOD],useNA = "a")
table(grill.dt[is.na(as.numeric(LOQ)),LOQ],useNA = "a")
grill.dt[is.na(LOQ) & !is.na(LOD),LOQ:=3.3*LOD]
grill.dt[!is.na(LOQ) & is.na(LOD),LOD:=LOQ/3.3]
table(grill.dt[is.na(as.numeric(LOD)),LOD],useNA = "a")
table(grill.dt[is.na(as.numeric(LOQ)),LOQ],useNA = "a")
table(grill.dt[is.na(as.numeric(min)),min])
grill.dt[is.na(as.numeric(min)) & !is.na(min),min:=0]

table(grill.dt[is.na(as.numeric(max)),max],useNA = "a")
grill.dt[is.na(as.numeric(max)) & !is.na(max) & !is.na(LOQ),max:=LOQ]
table(grill.dt[is.na(as.numeric(max)),max],useNA = "a")
grill.dt[is.na(as.numeric(max)) & !is.na(max)]

LOQ.table<-merge(
  setNames(merge(
    grill.dt[!is.na(LOQ) & grilling.method %in% included.methods][,sum(n.obs,na.rm=TRUE),by=c("food.name","chem.name")],
    grill.dt[!is.na(LOQ) & !grilling.method %in% included.methods][,sum(n.obs,na.rm=TRUE),by=c("food.name","chem.name")],
    by=c("food.name","chem.name"),all = TRUE,suffixes = c(".included",".excluded")),c("food.name","chem.name","N.included","N.excluded")),
  setNames(merge(
    grill.dt[!is.na(LOQ) & grilling.method %in% included.methods][,mean(LOQ),by=c("food.name","chem.name")],
    grill.dt[!is.na(LOQ) & !grilling.method %in% included.methods][,mean(LOQ),by=c("food.name","chem.name")],
    by=c("food.name","chem.name"),all = TRUE,suffixes = c(".included",".excluded")),c("food.name","chem.name","Mean.LOQ.included","Mean.LOQ.excluded")),
  by=c("chem.name","food.name"),all = TRUE)
fwrite(LOQ.table,"LOQ.summary.table.txt",sep="\t")

## populate missing LOQs with chem/food averages
grill.dt<-merge(grill.dt,LOQ.table[,.(chem.name,food.name,Mean.LOQ.included,Mean.LOQ.excluded)],by=c("chem.name","food.name"),all.x = TRUE)
grill.dt[,(c("out.up","out.down")):=NULL]
grill.dt[,LOQ:=ifelse(is.na(LOQ),ifelse(grilling.method %in% included.methods,Mean.LOQ.included,Mean.LOQ.excluded),LOQ)]
grill.dt[,LOQ:=ifelse(is.na(LOQ),Mean.LOQ.included,LOQ)]
grill.dt[,Mean.LOQ.included:=NULL]
grill.dt[,Mean.LOQ.excluded:=NULL]

## fill in still missing max's
table(grill.dt[is.na(as.numeric(max)),max],useNA = "a")
grill.dt[is.na(as.numeric(max)) & !is.na(max) & !is.na(LOQ),max:=LOQ]
table(grill.dt[is.na(as.numeric(max)),max],useNA = "a")

str(grill.dt)
## remove non-numeric SD
grill.dt[!is.na(SD) & is.na(as.numeric(SD)),SD:=ifelse(grepl("nd",SD,ignore.case = TRUE),NA,SD)]
if (nrow(grill.dt[!is.na(SD) & is.na(as.numeric(SD))])>0){
  stop("Still some non-numeric SDs")
} else {
  grill.dt[,SD:=as.numeric(SD)]
}
if (nrow(grill.dt[!is.na(min) & is.na(as.numeric(min))])>0){
  stop("Still some non-numeric min")
} else {
  grill.dt[,min:=as.numeric(min)]
}
if (nrow(grill.dt[!is.na(max) & is.na(as.numeric(max))])>0){
  stop("Still some non-numeric max")
} else {
  grill.dt[,max:=as.numeric(max)]
}
str(grill.dt)

if (nrow(grill.dt[!is.na(single) & is.na(as.numeric(single)) & single!=nd])>0){
  stop("some unexpected non-numeric values in single")
} else {
  grill.dt[!is.na(single) & is.na(as.numeric(single)) & single==nd,single:=NA]
  grill.dt[,single:=as.numeric(single)]
}
if (nrow(grill.dt[Mean!=single])>0){
  stop("some conflicting values in mean and single")
} else {
  grill.dt[,single:=NULL]
}

unlist(lapply(grill.dt,class))

## fix inconsistent min/max inclusion
grill.dt[xor(is.na(max),is.na(min)),min:=ifelse(!is.na(nd) & is.na(min),0,min)]
grill.dt[xor(is.na(max),is.na(min)),max:=ifelse(!is.na(nd) & is.na(max),LOQ,max)]
if (nrow(grill.dt[xor(is.na(max),is.na(min)),])>0){
  stop("still some unfixed min/max pairs")
}

## lines with not enough info
grill.dt[n.obs!=1 & is.na(max) & is.na(min) & (is.na(Mean) | is.na(SD)) & is.na(nd) & is.na(nd),no.info:=TRUE]
grill.dt[is.na(no.info),no.info:=FALSE]
fwrite(grill.dt[no.info==TRUE],
  "lines.with.missing.info.txt",sep = "\t")

## lines with inconsistent min/max vs Mean/median
grill.dt[Mean>max]
grill.dt[Mean<min]
grill.dt[`pooled median`>max]
grill.dt[`pooled median`<min]
fwrite(setorder(grill.dt[((Mean>=max | Mean<=min) & min!=max) | min>max],Excel.line),"Mean.min.max.problems.txt",sep="\t")

# subset for uniform simulation
## non-detected mean
grill.dt[is.na(Mean) & !is.na(nd),min:=ifelse(is.na(min),0,min)]
grill.dt[is.na(Mean) & !is.na(nd),max:=ifelse(is.na(max),LOQ,max)]
grill.dt[is.na(Mean) & !is.na(min),uniform:=TRUE]
## missing SD and non-missing min/max
grill.dt[n.obs>3 & !is.na(min) & is.na(SD),uniform:=TRUE]
grill.dt[is.na(uniform),uniform:=FALSE]

## subset with individual results
# all n
grill.dt[!is.na(Mean) & is.na(SD) & is.na(min)]
grill.dt[!is.na(Mean) & is.na(SD) & is.na(min),individual:=TRUE]
n.max<-max(grill.dt[!is.na(Mean) & is.na(SD) & is.na(min),n.obs])
grill.dt[!is.na(Mean) & is.na(SD) & is.na(min)][,(paste0("x",1:max(n.max))):=NA]
for (i in 1:n.max){
  grill.dt[!is.na(Mean) & is.na(SD) & is.na(min) & n.obs>=i,(paste0("x",i)):=Mean]
}
# n=1
grill.dt[n.obs==1 & !is.na(Mean)]
grill.dt[n.obs==1 & !is.na(Mean),individual:=TRUE]
grill.dt[n.obs==1 & !is.na(Mean),x1:=Mean]
# n=2
grill.dt[n.obs==2 & !is.na(min) & uniform==FALSE]
grill.dt[n.obs==2 & !is.na(min) & uniform==FALSE,individual:=TRUE]
grill.dt[n.obs==2 & !is.na(min) & uniform==FALSE,x1:=n.obs*Mean-max]
grill.dt[n.obs==2 & !is.na(min) & uniform==FALSE,x2:=max]
grill.dt[n.obs==2 & !is.na(Mean) & !is.na(SD) & is.na(min) & uniform==FALSE]
grill.dt[n.obs==2 & !is.na(Mean) & !is.na(SD) & is.na(min) & uniform==FALSE,individual:=TRUE]
grill.dt[n.obs==2 & !is.na(Mean) & !is.na(SD) & is.na(min) & uniform==FALSE,max:=Mean+SD]
grill.dt[n.obs==2 & !is.na(Mean) & !is.na(SD) & is.na(min) & uniform==FALSE,min:=Mean-SD]
grill.dt[n.obs==2 & !is.na(Mean) & !is.na(SD) & uniform==FALSE,x1:=min]
grill.dt[n.obs==2 & !is.na(Mean) & !is.na(SD) & uniform==FALSE,x2:=max]
grill.dt[n.obs==2 & is.na(x2) & !is.na(max)]
grill.dt[n.obs==2 & is.na(x2) & !is.na(SD)]
grill.dt[n.obs==2 & is.na(x2) & !is.na(`pooled median`)]
if (nrow(grill.dt[min<0])>0){
  stop("generated lines with negative min")
}
#n=3
grill.dt[n.obs==3 & min!=0 & uniform==FALSE]
grill.dt[n.obs==3 & min!=0 & uniform==FALSE,individual:=TRUE]
grill.dt[n.obs==3 & min!=0 & uniform==FALSE,x1:=min]
grill.dt[n.obs==3 & min!=0 & uniform==FALSE,x2:=max]
grill.dt[n.obs==3 & min!=0 & uniform==FALSE,x3:=Mean*n.obs-min-max]
grill.dt[n.obs==3 & min==0 & !is.na(`pooled median`) & uniform==FALSE & is.na(individual)]
grill.dt[n.obs==3 & min==0 & !is.na(`pooled median`) & uniform==FALSE & is.na(individual),x1:=max]
grill.dt[n.obs==3 & min==0 & !is.na(`pooled median`) & uniform==FALSE & is.na(individual),x2:=`pooled median`]
grill.dt[n.obs==3 & min==0 & !is.na(`pooled median`) & uniform==FALSE & is.na(individual),x3:=Mean*n.obs-x1-x2]
grill.dt[n.obs==3 & min==0 & !is.na(`pooled median`) & uniform==FALSE & is.na(individual) & x3<0,(c(paste0("x",1:3),"individual")):=NA]
grill.dt[n.obs==3 & min==0 & !is.na(`pooled median`) & uniform==FALSE & is.na(individual) & x3>=0,individual:=TRUE]
grill.dt[n.obs==3 & !is.na(max) & !is.na(SD) & uniform==FALSE & is.na(individual)]
grill.dt[n.obs==3 & !is.na(max) & !is.na(SD) & uniform==FALSE & is.na(individual),x1:=max]
grill.dt[n.obs==3 & !is.na(max) & !is.na(SD) & uniform==FALSE & is.na(individual),x2:=(3*Mean-x1)/2+sqrt(((SD^2+Mean^2)*3-x1^2)/2-((3*Mean-x1)/2)^2)]
grill.dt[n.obs==3 & !is.na(max) & !is.na(SD) & uniform==FALSE & is.na(individual),x3:=(3*Mean-x1)/2-sqrt(((SD^2+Mean^2)*3-x1^2)/2-((3*Mean-x1)/2)^2)]
grill.dt[n.obs==3 & !is.na(max) & !is.na(SD) & uniform==FALSE & is.na(individual) & (x2<0 | x3<0 | is.nan(x2) | is.nan(x3)),(c(paste0("x",1:3),"individual")):=NA]
grill.dt[n.obs==3 & !is.na(max) & !is.na(SD) & uniform==FALSE & is.na(individual) & x2>0 & x3>0,individual:=TRUE]
grill.dt[n.obs==3 & min==0 & uniform==FALSE & is.na(individual)]
grill.dt[n.obs==3 & min==0 & uniform==FALSE & is.na(individual),x1:=max]
grill.dt[n.obs==3 & min==0 & uniform==FALSE & is.na(individual),x2:=ifelse(3*Mean-max< -Mean*.05,NA,3*Mean-max)]
grill.dt[n.obs==3 & min==0 & uniform==FALSE & is.na(individual) & x2<0,x2:=0]
grill.dt[n.obs==3 & min==0 & uniform==FALSE & is.na(individual),x3:=0]
grill.dt[n.obs==3 & min==0 & uniform==FALSE & is.na(individual) & is.na(x2),(c(paste0("x",1:3),"individual")):=NA]
grill.dt[n.obs==3 & min==0 & uniform==FALSE & is.na(individual) & !is.na(x2),individual:=TRUE]
#n=4
if (nrow(grill.dt[n.obs==4 & !is.na(min) & !is.na(SD) & uniform==FALSE])+
    nrow(grill.dt[n.obs==4 & !is.na(min) & !is.na(`pooled median`) & uniform==FALSE])>0){
  stop("there are some lines with n tests that can be split")
}
# otherwise
grill.dt[is.na(x1),individual:=FALSE]

## lognormal subset
grill.dt[individual==FALSE & uniform==FALSE & !is.na(Mean) & !is.na(SD)]
grill.dt[individual==FALSE & uniform==FALSE & !is.na(Mean) & !is.na(SD),lognormal:=TRUE]
grill.dt[is.na(lognormal),lognormal:=FALSE]

if (nrow(grill.dt[individual+uniform+lognormal>1])>0){
  stop("conflicted classification")
}
### remaining unclassifieds
grill.dt[individual+uniform+lognormal<1]

#########
## transform Mean/SD to logN equivalent
lognorm_E_SD <- function(mean,sd){
  list(log(mean^2/sqrt(sd^2+mean^2)),
    sqrt(log(1 + sd^2/mean^2)))
}

grill.dt[lognormal==TRUE,(c("logMean","logSD")):=lognorm_E_SD(Mean,SD)]

fwrite(grill.dt,"grill.dt.for.simulations.txt",sep="\t")
grill.dt<-fread("grill.dt.for.simulations.txt",sep="\t")

### Excluding "Wood/log fire"
if (alternative.inclusion==TRUE){
  fwrite(grill.dt[!grilling.method %in% "Wood/log fire"],"grill.dt.for.simulations.no.wood.log.txt",sep="\t")
  grill.dt<-fread("grill.dt.for.simulations.no.wood.log.txt",sep="\t")
}

########################################
# simulations
########################################
common.multiple<-100L # for each simulation run
sim.reps<-100L # number of simulation runs

set.seed(1234)
## all na's have to cleared out from n.obs: same rule as for lognormal case is applied
uniform.set<-grill.dt[uniform==TRUE][is.na(n.obs),n.obs:=3][rep(1:.N,n.obs*common.multiple)][rep(1:.N,times=sim.reps)][,
  Sim.Rep:=rep(1:sim.reps,each=.N/sim.reps)][,
    Sim.Val:=runif(.N,0,1)*(max-min)+min][,(grep("^x[0-9]+$",names(grill.dt),value = TRUE)):=NULL]
lognormal.set<-grill.dt[lognormal==TRUE][rep(1:.N,n.obs*common.multiple)][rep(1:.N,times=sim.reps)][,
  Sim.Rep:=rep(1:sim.reps,each=.N/sim.reps)][,
    Sim.Val:=exp(log(rlnorm(.N,0,1))*logSD+logMean)][,(grep("^x[0-9]+$",names(grill.dt),value = TRUE)):=NULL]
individual.set<-melt(grill.dt[individual==TRUE][rep(1:.N,each=common.multiple,times=sim.reps)][,Sim.Rep:=rep(1:sim.reps,each=.N/sim.reps)],
  value.name = "Sim.Val",measure.vars = grep("^x[0-9]+$",names(grill.dt),value = TRUE),na.rm = TRUE)[,
  variable:=NULL]

simulation.dataset<-rbindlist(list(uniform.set,lognormal.set,individual.set))
setkeyv(simulation.dataset,c("Sim.Rep","chem.name","food.name","Sim.Val"))
simulation.dataset[,Sim.Val:=Sim.Val+10^(-7)] # to remove zeros for lognormality

simulation.dataset[,included.method:=ifelse(grilling.method %in% included.methods,TRUE,FALSE)]
simulation.dataset[,Sim.Val.LB:=ifelse(Sim.Val<LOQ,0,Sim.Val)]
simulation.dataset[,Sim.Val.UB:=ifelse(Sim.Val<LOQ,LOQ,Sim.Val)]
simulation.dataset[,Sim.Val.rank:=frank(Sim.Val, ties.method="random"),by=c("chem.name","food.name","Sim.Rep","included.method")]
setkeyv(simulation.dataset,c("chem.name","food.name","Sim.Rep","included.method","Sim.Val.rank"))
if (alternative.inclusion==TRUE){
  fwrite(simulation.dataset,"simulation.dataset.dt.alternative.inclusion.txt",sep="\t")
  # simulation.dataset<-fread("simulation.dataset.dt.alternative.inclusion.txt",sep="\t")
} else {
  fwrite(simulation.dataset,"simulation.dataset.dt.txt",sep="\t")
  # simulation.dataset<-fread("simulation.dataset.dt.txt",sep="\t")
}


## theoretical (lognormal) percentiles
lnorm.parameter.dt<-simulation.dataset[,{fd<-try(fitdist(Sim.Val,"lnorm"),silent = TRUE); if (class(fd)=="try-error"){as.list(as.double(c(NA,NA)))} else {as.list(fd$estimate)}},
  by=c("chem.name","food.name","Sim.Rep","included.method")]
if (alternative.inclusion==TRUE){
  fwrite(lnorm.parameter.dt,"lnorm.parameter.dt.alternative.inclusion.txt",sep="\t")
} else {
  fwrite(lnorm.parameter.dt,"lnorm.parameter.dt.txt",sep="\t")
}

percentile.dt<-copy(lnorm.parameter.dt)[,
  (sprintf("p%01.2f",c(.05,.1,.25,.50,.75,.90,.95))):=
    transpose(lapply(seq_along(meanlog),function(r) (qlnorm(c(.05,.1,.25,.50,.75,.90,.95),meanlog = meanlog[r],sdlog=sdlog[r]))))]
percentile.dt[,Mean:=exp(meanlog+sdlog^2/2)]
percentile.dt<-percentile.dt[,as.list(unlist(lapply(.SD,function(c) unlist(quantile(c,c(.05,.5,.95),na.rm = TRUE))))),
  by=c("chem.name","food.name","included.method"),.SDcols = c("Mean",sprintf("p%01.2f",c(.05,.1,.25,.50,.75,.90,.95)))]

## empirical percentiles
setkeyv(simulation.dataset,c("chem.name","food.name","Sim.Rep","included.method","Sim.Val"))
setorder(simulation.dataset,chem.name,food.name,Sim.Rep,included.method,Sim.Val)
emp.percentile.dt<-simulation.dataset[,c(Mean=mean(Sim.Val),as.list(quantile(Sim.Val,c(.05,.1,.25,.50,.75,.90,.95)))),
  by=c("chem.name","food.name","Sim.Rep","included.method")]
emp.percentile.dt[,as.list(unlist(lapply(.SD,function(c) c(min=min(c),unlist(quantile(c,c(.05,.5,.95))),max=max(c))))),
  by=c("chem.name","food.name","included.method"),.SDcols = c("Mean",sprintf("%d%%",c(.05,.1,.25,.50,.75,.90,.95)*100))]
emp.percentile.dt<-emp.percentile.dt[,as.list(unlist(lapply(.SD,function(c) c(unlist(quantile(c,c(.05,.5,.95))))))),
  by=c("chem.name","food.name","included.method"),.SDcols = c("Mean",sprintf("%d%%",c(.05,.1,.25,.50,.75,.90,.95)*100))]
# LB
emp.percentile.LB.dt<-simulation.dataset[,c(Mean=mean(Sim.Val.LB),as.list(quantile(Sim.Val,c(.05,.1,.25,.50,.75,.90,.95)))),
  by=c("chem.name","food.name","Sim.Rep","included.method")]
emp.percentile.LB.dt[,as.list(unlist(lapply(.SD,function(c) c(min=min(c),unlist(quantile(c,c(.05,.5,.95))),max=max(c))))),
  by=c("chem.name","food.name","included.method"),.SDcols = c("Mean",sprintf("%d%%",c(.05,.1,.25,.50,.75,.90,.95)*100))]
emp.percentile.LB.dt<-emp.percentile.LB.dt[,as.list(unlist(lapply(.SD,function(c) c(unlist(quantile(c,c(.05,.5,.95))))))),
  by=c("chem.name","food.name","included.method"),.SDcols = c("Mean",sprintf("%d%%",c(.05,.1,.25,.50,.75,.90,.95)*100))]
# UB
emp.percentile.UB.dt<-simulation.dataset[,c(Mean=mean(Sim.Val.UB),as.list(quantile(Sim.Val,c(.05,.1,.25,.50,.75,.90,.95)))),
  by=c("chem.name","food.name","Sim.Rep","included.method")]
emp.percentile.UB.dt[,as.list(unlist(lapply(.SD,function(c) c(min=min(c),unlist(quantile(c,c(.05,.5,.95))),max=max(c))))),
  by=c("chem.name","food.name","included.method"),.SDcols = c("Mean",sprintf("%d%%",c(.05,.1,.25,.50,.75,.90,.95)*100))]
emp.percentile.UB.dt<-emp.percentile.UB.dt[,as.list(unlist(lapply(.SD,function(c) c(unlist(quantile(c,c(.05,.5,.95))))))),
  by=c("chem.name","food.name","included.method"),.SDcols = c("Mean",sprintf("%d%%",c(.05,.1,.25,.50,.75,.90,.95)*100))]

### compute 4PAH distribution
sum.dt<-copy(percentile.dt)[,
  lapply(.SD,function(c) ifelse(.N==4,sum(c),as.numeric(NA))),
  by=c("food.name","included.method"),.SDcols=setdiff(names(percentile.dt),c("chem.name","food.name","included.method"))][,
    chem.name:="4PAH"]
emp.sum.dt<-copy(emp.percentile.dt)[,
  lapply(.SD,function(c) ifelse(.N==4,sum(c),as.numeric(NA))),
  by=c("food.name","included.method"),.SDcols=setdiff(names(emp.percentile.dt),c("chem.name","food.name","included.method"))][,
    chem.name:="4PAH"]
emp.sum.LB.dt<-copy(emp.percentile.LB.dt)[,
  lapply(.SD,function(c) ifelse(.N==4,sum(c),as.numeric(NA))),
  by=c("food.name","included.method"),.SDcols=setdiff(names(emp.percentile.LB.dt),c("chem.name","food.name","included.method"))][,
    chem.name:="4PAH"]
emp.sum.UB.dt<-copy(emp.percentile.UB.dt)[,
  lapply(.SD,function(c) ifelse(.N==4,sum(c),as.numeric(NA))),
  by=c("food.name","included.method"),.SDcols=setdiff(names(emp.percentile.UB.dt),c("chem.name","food.name","included.method"))][,
    chem.name:="4PAH"]

## gathering percentiles and turning tables wide\
perc.to.wide<-function(percentile.dt,sum.dt=data.table(),dcast.str="food.name+included.method~...",name.str="chem.name"){
  DT<-rbindlist(list(percentile.dt,sum.dt),use.names = TRUE,fill = TRUE)
  ex.cols<-c(name.str,strsplit(gsub("~.+","",dcast.str),"\\+")[[1]])
  percentile.dt.wide<-dcast(DT,as.formula(dcast.str),
    value.var = setdiff(names(DT),ex.cols))
  setnames(percentile.dt.wide,grep("_",names(percentile.dt.wide),value = TRUE),
    gsub("(.+)_(.+)","\\2_\\1",grep("_",names(percentile.dt.wide),value = TRUE)))
  type.cols=setdiff(names(DT),ex.cols)
  setcolorder(percentile.dt.wide,paste(rep(sort(unique(DT[,get(name.str)])),each=length(type.cols)),type.cols,sep="_"))
  setcolorder(percentile.dt.wide,strsplit(gsub("~.+","",dcast.str),"\\+")[[1]])
  return(percentile.dt.wide)
}
percentile.dt.wide<-perc.to.wide(percentile.dt,sum.dt)
emp.percentile.dt.wide<-perc.to.wide(emp.percentile.dt,emp.sum.dt)
emp.percentile.LB.dt.wide<-perc.to.wide(emp.percentile.LB.dt,emp.sum.LB.dt)
emp.percentile.UB.dt.wide<-perc.to.wide(emp.percentile.UB.dt,emp.sum.UB.dt)

# write to Excel
wb <- createWorkbook()
addWorksheet(wb, "logN Percentiles")
addWorksheet(wb, "Percentiles")
addWorksheet(wb, "LB Percentiles")
addWorksheet(wb, "UB Percentiles")
for (i in 1:4){
  freezePane(
    wb,i,
    firstActiveRow = 2,
    firstActiveCol = 3
  )
  setColWidths(wb, i, cols = 1:ncol(percentile.dt.wide), widths = 17)
}
writeData(wb, 1, x = percentile.dt.wide, withFilter = TRUE)
writeData(wb, 2, x = emp.percentile.dt.wide, withFilter = TRUE)
writeData(wb, 3, x = emp.percentile.LB.dt.wide, withFilter = TRUE)
writeData(wb, 4, x = emp.percentile.UB.dt.wide, withFilter = TRUE)

if (alternative.inclusion==TRUE){
  fname<-paste0("Grill.Occurrence.Perc.alternative.inclusion.",gsub("-","",Sys.Date()))
} else {
  fname<-paste0("Grill.Occurrence.Perc.",gsub("-","",Sys.Date()))
}
i=1
while(file.exists(paste0(fname,"_",i,".xlsx"))){
  i=i+1
}
saveWorkbook(wb, file = paste0(fname,"_",i,".xlsx"),overwrite = TRUE)


### exposure

# Checking units 
table(grill.lit.dt[,`PAH conc units`])
table(grill.lit.dt[,`units for LOD`])
table(grill.lit.dt[,`units for LOQ`])

grill.lit.dt[(!grepl("(µg/kg)|(µg/L)|(ng/g)|(ng/mL)|(ppb)",`PAH conc units`) & !`PAH conc units` %in% c("",NA))|
    (!grepl("(µg/kg)|(µg/L)|(ng/g)|(ng/mL)|(ppb)",`units for LOD`) & !`units for LOD` %in% c("",NA))|
    (!grepl("(µg/kg)|(µg/L)|(ng/g)|(ppb)",`units for LOQ`) & !`units for LOQ` %in% c("",NA))]

table(grill.dt[,food.name])
# all is in ng/g (or equivalent)

plates.dt<-rbindlist(
  list(data.table(food.name=c("Pork_lean","CWoS","Beef","Fish_salmon"),weights=50,type="lean"),
    data.table(food.name=c("Sausage","Beef patties","Pork_fatty","CWS"),weights=c(60,60,60,20),type="fat")))

# mean-based exposure
mean.dt<-dcast(
  copy(grill.dt)[,included.method:=grilling.method %in% included.methods][,
    weighted.mean(ifelse(!is.na(Mean),Mean,(min+max)/2),n.obs,na.rm=TRUE),by=c("food.name","chem.name","included.method")],
  food.name+included.method~...,value.var = "V1")[,
    `4PAH`:=BaA+BaP+BbF+Chry]
mean.LB.dt<-dcast(
  copy(grill.dt)[,included.method:=grilling.method %in% included.methods][Mean<LOQ,Mean:=0][min<LOQ,min:=0][max<=LOQ,max:=0][,
    weighted.mean(ifelse(!is.na(Mean) ,Mean,(min+max)/2),n.obs,na.rm=TRUE),by=c("food.name","chem.name","included.method")],
  food.name+included.method~...,value.var = "V1")[,
    `4PAH`:=BaA+BaP+BbF+Chry]
mean.UB.dt<-dcast(
  copy(grill.dt)[,included.method:=grilling.method %in% included.methods][Mean<LOQ,Mean:=LOQ][min<LOQ,min:=LOQ][max<=LOQ,max:=LOQ][,
    weighted.mean(ifelse(!is.na(Mean),Mean,(min+max)/2),n.obs,na.rm=TRUE),by=c("food.name","chem.name","included.method")],
  food.name+included.method~...,value.var = "V1")[,
    `4PAH`:=BaA+BaP+BbF+Chry]
setorder(mean.dt,-included.method,food.name)
setorder(mean.LB.dt,-included.method,food.name)
setorder(mean.UB.dt,-included.method,food.name)

mean.plate.dt<-merge(mean.dt,plates.dt,by="food.name")[,
  (c("BaA","BaP","BbF","Chry","4PAH")):=lapply(.SD,function(col) col*weights),.SDcols=c("BaA","BaP","BbF","Chry","4PAH")]
mean.plate.dt<-rbindlist(list(mean.plate.dt[,lapply(.SD,sum),by=c("included.method","type"),
  .SDcols = setdiff(names(mean.plate.dt),c("included.method","type","food.name"))][,food.name:="Plate Sum"],mean.plate.dt),
  use.names = TRUE,fill = TRUE)
setcolorder(mean.plate.dt,c("food.name","weights","included.method","type"))
setorder(mean.plate.dt,-included.method,type,weights,food.name)
mean.LB.plate.dt<-merge(mean.LB.dt,plates.dt,by="food.name")[,
  (c("BaA","BaP","BbF","Chry","4PAH")):=lapply(.SD,function(col) col*weights),.SDcols=c("BaA","BaP","BbF","Chry","4PAH")]
mean.LB.plate.dt<-rbindlist(list(mean.LB.plate.dt[,lapply(.SD,sum),by=c("included.method","type"),
  .SDcols = setdiff(names(mean.LB.plate.dt),c("included.method","type","food.name"))][,food.name:="Plate Sum"],mean.LB.plate.dt),
  use.names = TRUE,fill = TRUE)
setcolorder(mean.LB.plate.dt,c("food.name","weights","included.method","type"))
setorder(mean.LB.plate.dt,-included.method,type,weights,food.name)
mean.UB.plate.dt<-merge(mean.UB.dt,plates.dt,by="food.name")[,
  (c("BaA","BaP","BbF","Chry","4PAH")):=lapply(.SD,function(col) col*weights),.SDcols=c("BaA","BaP","BbF","Chry","4PAH")]
mean.UB.plate.dt<-rbindlist(list(mean.UB.plate.dt[,lapply(.SD,sum),by=c("included.method","type"),
  .SDcols = setdiff(names(mean.UB.plate.dt),c("included.method","type","food.name"))][,food.name:="Plate Sum"],mean.UB.plate.dt),
  use.names = TRUE,fill = TRUE)
setcolorder(mean.UB.plate.dt,c("food.name","weights","included.method","type"))
setorder(mean.UB.plate.dt,-included.method,type,weights,food.name)

## lognormal exposure
if (alternative.inclusion==TRUE){
  lnorm.parameter.dt<-fread("lnorm.parameter.dt.alternative.inclusion.txt",sep="\t")
} else {
  lnorm.parameter.dt<-fread("lnorm.parameter.dt.txt",sep="\t")
}
set.seed(1234)
logN.dt<-merge(lnorm.parameter.dt,
  setnames(unique(lnorm.parameter.dt[,.(food.name,Sim.Rep)])[,list(1:1000,runif(1000)),by=c("food.name","Sim.Rep")],c("V1","V2"),c("lineid","CDF.probs")),
  by=c("food.name","Sim.Rep"),allow.cartesian = TRUE)
logN.dt[,Sim.Val:=qlnorm(CDF.probs,meanlog = meanlog,sdlog=sdlog)]
logN.dt<-merge(logN.dt,plates.dt,by="food.name")[,Sim.Val:=Sim.Val*weights]
logN.dt<-rbindlist(list(logN.dt,
    setnames(logN.dt[,list(ifelse(.N==4,weights[1],as.numeric(NA)),ifelse(.N==4,sum(Sim.Val),as.numeric(NA))),
      by=c("food.name","Sim.Rep","included.method","CDF.probs","lineid","type")][,chem.name:="4PAH"],
    c("V1","V2"),c("weights","Sim.Val"))
),use.names = TRUE,fill = TRUE)
logN.dt<-rbindlist(list(logN.dt,
  setnames(logN.dt[,list(sum(weights),sum(Sim.Val)),by=c("Sim.Rep","type","chem.name","lineid","included.method")],c("V1","V2"),c("weights","Sim.Val"))[,
    food.name:="Plate Sum"]),
  use.names = TRUE,fill = TRUE)

plate.logN.dt<-logN.dt[,c(Mean=mean(Sim.Val),as.list(quantile(Sim.Val,c(.05,.1,.25,.50,.75,.90,.95)))),
  by=c("chem.name","food.name","Sim.Rep","included.method","type","weights")]
plate.logN.dt[,as.list(unlist(lapply(.SD,function(c) c(min=min(c),unlist(quantile(c,c(.05,.5,.95))),max=max(c))))),
  by=c("chem.name","food.name","included.method"),.SDcols = c("Mean",sprintf("%d%%",c(.05,.1,.25,.50,.75,.90,.95)*100))]
plate.logN.dt<-plate.logN.dt[,as.list(unlist(lapply(.SD,function(c) c(unlist(quantile(c,c(.05,.5,.95))))))),
  by=c("chem.name","food.name","included.method","type","weights"),.SDcols = c("Mean",sprintf("%d%%",c(.05,.1,.25,.50,.75,.90,.95)*100))]

logN.dt.wide<-perc.to.wide(plate.logN.dt,dcast.str="food.name+weights+included.method+type~...")
setcolorder(logN.dt.wide,c("food.name","weights","included.method","type"))
setorder(logN.dt.wide,-included.method,type,weights,food.name)  

### "empirical" exposure from simulations
simulated.plate.dt<-simulation.dataset[food.name %in% plates.dt[,food.name]]
logN.CDFs<-setnames(unique(logN.dt[!is.na(CDF.probs),.(food.name,CDF.probs,weights,type,lineid,Sim.Rep)]),
  c("food.name","Sim.Rep"),c("food.name.CDF","Sim.Rep.CDF"))
simulated.plate.dt<-merge(merge(
  setnames(simulated.plate.dt[,
    quantile(Sim.Val,logN.CDFs[food.name.CDF==food.name & Sim.Rep.CDF==Sim.Rep,CDF.probs]),
    by=c("food.name","chem.name","included.method","Sim.Rep")],"V1","Sim.Val")[,lineid:=1:.N,by=c("food.name","chem.name","included.method","Sim.Rep")],
  setnames(simulated.plate.dt[,
    quantile(Sim.Val.LB,logN.CDFs[food.name.CDF==food.name & Sim.Rep.CDF==Sim.Rep,CDF.probs]),
    by=c("food.name","chem.name","included.method","Sim.Rep")],"V1","Sim.Val.LB")[,lineid:=1:.N,by=c("food.name","chem.name","included.method","Sim.Rep")],
  by=c("food.name","chem.name","included.method","Sim.Rep","lineid")),
  setnames(simulated.plate.dt[,
    quantile(Sim.Val.UB,logN.CDFs[food.name.CDF==food.name & Sim.Rep.CDF==Sim.Rep,CDF.probs]),
    by=c("food.name","chem.name","included.method","Sim.Rep")],"V1","Sim.Val.UB")[,lineid:=1:.N,by=c("food.name","chem.name","included.method","Sim.Rep")],
  by=c("food.name","chem.name","included.method","Sim.Rep","lineid"))
simulated.plate.dt<-merge(simulated.plate.dt,plates.dt,by="food.name")
simulated.plate.dt[,(c("Sim.Val", "Sim.Val.LB", "Sim.Val.UB")):=lapply(.SD,function(x) x*weights),.SDcols = c("Sim.Val", "Sim.Val.LB", "Sim.Val.UB")]
simulated.plate.dt<-merge(simulated.plate.dt,logN.CDFs[,.(food.name.CDF,CDF.probs,Sim.Rep.CDF,lineid)],
  by.x=c("food.name","Sim.Rep","lineid"),by.y=c("food.name.CDF","Sim.Rep.CDF","lineid"))

simulated.plate.dt<-rbindlist(list(simulated.plate.dt,
  setnames(simulated.plate.dt[,
    list(ifelse(.N==4,weights[1],as.numeric(NA)),
      ifelse(.N==4,sum(Sim.Val),as.numeric(NA)),
      ifelse(.N==4,sum(Sim.Val.LB),as.numeric(NA)),
      ifelse(.N==4,sum(Sim.Val.UB),as.numeric(NA))),
    by=c("food.name","Sim.Rep","included.method","CDF.probs","lineid","type")][,chem.name:="4PAH"],
    paste0("V",1:4),c("weights","Sim.Val","Sim.Val.LB","Sim.Val.UB"))
),use.names = TRUE,fill = TRUE)
simulated.plate.dt<-rbindlist(list(simulated.plate.dt,
  setnames(simulated.plate.dt[,list(sum(weights),sum(Sim.Val),sum(Sim.Val.LB),sum(Sim.Val.UB)),
    by=c("Sim.Rep","type","chem.name","lineid","included.method")],paste0("V",1:4),c("weights","Sim.Val","Sim.Val.LB","Sim.Val.UB"))[,
    food.name:="Plate Sum"]),
  use.names = TRUE,fill = TRUE)

simulated.plate.MB.dt<-simulated.plate.dt[,c(Mean=mean(Sim.Val),as.list(quantile(Sim.Val,c(.05,.1,.25,.50,.75,.90,.95)))),
  by=c("chem.name","food.name","Sim.Rep","included.method","type","weights")]
simulated.plate.MB.dt[,as.list(unlist(lapply(.SD,function(c) c(min=min(c),unlist(quantile(c,c(.05,.5,.95))),max=max(c))))),
  by=c("chem.name","food.name","included.method"),.SDcols = c("Mean",sprintf("%d%%",c(.05,.1,.25,.50,.75,.90,.95)*100))]
simulated.plate.MB.dt<-simulated.plate.MB.dt[,as.list(unlist(lapply(.SD,function(c) c(unlist(quantile(c,c(.05,.5,.95))))))),
  by=c("chem.name","food.name","included.method","type","weights"),.SDcols = c("Mean",sprintf("%d%%",c(.05,.1,.25,.50,.75,.90,.95)*100))]

simulated.plate.MB.dt.wide<-perc.to.wide(simulated.plate.MB.dt,dcast.str="food.name+weights+included.method+type~...")
setcolorder(simulated.plate.MB.dt.wide,c("food.name","weights","included.method","type"))
setorder(simulated.plate.MB.dt.wide,-included.method,type,weights,food.name)  

simulated.plate.LB.dt<-simulated.plate.dt[,c(Mean=mean(Sim.Val.LB),as.list(quantile(Sim.Val.LB,c(.05,.1,.25,.50,.75,.90,.95)))),
  by=c("chem.name","food.name","Sim.Rep","included.method","type","weights")]
simulated.plate.LB.dt[,as.list(unlist(lapply(.SD,function(c) c(min=min(c),unlist(quantile(c,c(.05,.5,.95))),max=max(c))))),
  by=c("chem.name","food.name","included.method"),.SDcols = c("Mean",sprintf("%d%%",c(.05,.1,.25,.50,.75,.90,.95)*100))]
simulated.plate.LB.dt<-simulated.plate.LB.dt[,as.list(unlist(lapply(.SD,function(c) c(unlist(quantile(c,c(.05,.5,.95))))))),
  by=c("chem.name","food.name","included.method","type","weights"),.SDcols = c("Mean",sprintf("%d%%",c(.05,.1,.25,.50,.75,.90,.95)*100))]

simulated.plate.LB.dt.wide<-perc.to.wide(simulated.plate.LB.dt,dcast.str="food.name+weights+included.method+type~...")
setcolorder(simulated.plate.LB.dt.wide,c("food.name","weights","included.method","type"))
setorder(simulated.plate.LB.dt.wide,-included.method,type,weights,food.name)  

simulated.plate.UB.dt<-simulated.plate.dt[,c(Mean=mean(Sim.Val.UB),as.list(quantile(Sim.Val.UB,c(.05,.1,.25,.50,.75,.90,.95)))),
  by=c("chem.name","food.name","Sim.Rep","included.method","type","weights")]
simulated.plate.UB.dt[,as.list(unlist(lapply(.SD,function(c) c(min=min(c),unlist(quantile(c,c(.05,.5,.95))),max=max(c))))),
  by=c("chem.name","food.name","included.method"),.SDcols = c("Mean",sprintf("%d%%",c(.05,.1,.25,.50,.75,.90,.95)*100))]
simulated.plate.UB.dt<-simulated.plate.UB.dt[,as.list(unlist(lapply(.SD,function(c) c(unlist(quantile(c,c(.05,.5,.95))))))),
  by=c("chem.name","food.name","included.method","type","weights"),.SDcols = c("Mean",sprintf("%d%%",c(.05,.1,.25,.50,.75,.90,.95)*100))]

simulated.plate.UB.dt.wide<-perc.to.wide(simulated.plate.UB.dt,dcast.str="food.name+weights+included.method+type~...")
setcolorder(simulated.plate.UB.dt.wide,c("food.name","weights","included.method","type"))
setorder(simulated.plate.UB.dt.wide,-included.method,type,weights,food.name)  

# write plate exposure to Excel
wb <- createWorkbook()
addWorksheet(wb, "Deterministic Plate Exposure")
addWorksheet(wb, "Deterministic Plate LB Exposure")
addWorksheet(wb, "Deterministic Plate UB Exposure")
addWorksheet(wb, "logN Plate Exposure")
addWorksheet(wb, "Simulated Plate Exposure")
addWorksheet(wb, "Simulated Plate LB Exposure")
addWorksheet(wb, "Simulated Plate UB Exposure")
addWorksheet(wb, "Deterministic Occ per Gram")
addWorksheet(wb, "Deterministic Occ LB per Gram")
addWorksheet(wb, "Deterministic Occ UB per Gram")
n.ws<-length(sheets(wb))
for (i in 1:7){
  freezePane(
    wb,i,
    firstActiveRow = 2,
    firstActiveCol = 5
  )
}
for (i in 8:n.ws){
  freezePane(
    wb,i,
    firstActiveRow = 2,
    firstActiveCol = 3
  )
}
for (i in 1:3){
  setColWidths(wb, i, cols = 1:ncol(mean.plate.dt), widths = 17)
}
for (i in 4:7){
  setColWidths(wb, i, cols = 1:ncol(logN.dt.wide), widths = 17)
}
for (i in 8:n.ws){
  setColWidths(wb, i, cols = 1:ncol(mean.dt), widths = 17)
}

dt.list<-list(
  mean.plate.dt,mean.LB.plate.dt,mean.UB.plate.dt,
  logN.dt.wide,
  simulated.plate.MB.dt.wide,simulated.plate.LB.dt.wide,simulated.plate.UB.dt.wide,
  mean.dt,mean.LB.dt,mean.UB.dt)
if (length(dt.list)!=n.ws){
  stop("Inconsistent number of data elements")
}
for (i in 1:length(dt.list)){
  writeData(wb, i, x = dt.list[[i]], withFilter = TRUE)
}

if (alternative.inclusion==TRUE){
  fname<-paste0("Grill.Plates.Perc.alternative.inclusion.",gsub("-","",Sys.Date()))
} else {
  fname<-paste0("Grill.Plates.Perc.",gsub("-","",Sys.Date()))
}
i=1
while(file.exists(paste0(fname,"_",i,".xlsx"))){
  i=i+1
}
saveWorkbook(wb, file = paste0(fname,"_",i,".xlsx"),overwrite = TRUE)

## summary table for raw data in requested format
grill.dt[max<=LOQ,n.nd:=n.obs]
grill.dt[min>=LOQ,n.nd:=0]
grill.dt[Mean<LOQ & n.obs==1,n.nd:=1]
grill.dt[Mean>=LOQ & n.obs==1,n.nd:=0]
grill.dt[Mean-3*SD>=LOQ,n.nd:=0]
grill.dt[Mean>=LOQ & is.na(SD) & is.na(n.nd) & is.na(min),n.nd:=0]
grill.dt[Mean<LOQ & is.na(SD) & is.na(n.nd) & is.na(min),n.nd:=n.obs]
grill.dt[min<LOQ & max>LOQ & is.na(n.nd), n.nd:=ceiling((LOQ-min)/(max-min)*n.obs)]
grill.dt[Mean+3*SD<=LOQ,n.nd:=n.obs]
## lose rule for n.nd assignment
grill.dt[Mean<LOQ & is.na(n.nd),n.nd:=1]
grill.dt[is.na(n.nd),n.nd:=0]

smry.table<-
  dcast(grill.dt[grilling.method %in% included.methods][,
  .(
    sum(n.obs,na.rm = TRUE),
    round(sum(n.nd,na.rm = TRUE)/sum(n.obs,na.rm = TRUE)*100),
    median(LOQ),
    weighted.mean(ifelse(!is.na(Mean),Mean,(min+max)/2),n.obs),
    weighted.mean(ifelse(!is.na(Mean),ifelse(Mean<LOQ,0,Mean),(ifelse(min<LOQ,0,min)+ifelse(max<=LOQ,0,max))/2),n.obs),
    weighted.mean(ifelse(!is.na(Mean),ifelse(Mean<LOQ,LOQ,Mean),(ifelse(min<LOQ,LOQ,min)+ifelse(max<LOQ,LOQ,max))/2),n.obs)
  ),by=c("food.name","chem.name")],
  food.name~...,value.var = paste0("V",1:6))
setnames(smry.table,names(smry.table),gsub("(.+)_(.+)","\\2_\\1",names(smry.table)))
setcolorder(smry.table,sort(names(smry.table)))
setcolorder(smry.table,"food.name")
for (i in 1:6){
  setnames(smry.table,
    names(smry.table),
    gsub(paste0("V",i),c("N","% under LOQ","Median LOQ","Mean","Mean LB","Mean UB")[i],names(smry.table)))
}

# run summary table to Excel
wb <- createWorkbook()
addWorksheet(wb, "Summary Table")
for (i in 1){
  freezePane(
    wb,i,
    firstActiveRow = 2,
    firstActiveCol = 4
  )
  setColWidths(wb, i, cols = 1:ncol(smry.table), widths = 17)
}

writeData(wb, i, x = smry.table, withFilter = TRUE)

if (alternative.inclusion==TRUE){
  fname<-paste0("Summary.Table.alternative.inclusion.",gsub("-","",Sys.Date()))
} else {
  fname<-paste0("Summary.Table.",gsub("-","",Sys.Date()))
}
i=1
while(file.exists(paste0(fname,"_",i,".xlsx"))){
  i=i+1
}
saveWorkbook(wb, file = paste0(fname,"_",i,".xlsx"),overwrite = TRUE)