R/pvalue_functions.R

In joshuamschmidt/multiPermr: What the Package Does in One 'Title Case' Line

calculate_empP <- function(permDT) {
  tp <- permDT[,.N,by=.(id,set_pr)][order(id,-set_pr)]
  tp[,cumN:=cumsum(N),by=id]
  tp[,enriched_p:=cumN/max(cumN)]
  tp <- tp[order(id,set_pr)]
  tp[,cumN:=cumsum(N),by=id]
  tp[,depleted_p:=cumN/max(cumN)]
  tp[,two_tailed_p:=2*min(enriched_p,depleted_p),by = row.names(tp)]
  tp[two_tailed_p > 1,two_tailed_p:=1]
  permDT <- permDT[tp,on=.(id,set_pr)]
  permDT[,`:=`(N=NULL,cumN=NULL)]
  return(permDT[])
}

# two tailed is wrong: see here: https://stats.stackexchange.com/questions/34052/two-sided-permutation-test-vs-two-one-sided
correct_fdr <- function(permDT,p_col,n_perm) {
  # fdr correction
  tfdr_obs <- permDT[perm_n==0][order(get(p_col))][,.(p=get(p_col),obs=frank(get(p_col),ties.method="max"))]
  tfdr_exp <- unique(permDT[order(get(p_col))][,.(p=get(p_col),exp=frank(get(p_col),ties.method="max")/n_perm)])
  fdrDT <- tfdr_exp[tfdr_obs,on=.(p)][,.(p,fdr=exp/obs)]
  
}

test_f <- function(x) {
  names(x)
}

correct_fdr_col <- function(p_values,n_set,n_perm) {
  # fdr correction
  tfdr_obs <- data.table(p=p_values[1:n_set],obs=data.table::frank(p_values[1:n_set],ties.method="max"))
  tfdr_exp <- unique(data.table(p=p_values,exp=data.table::frank(p_values,ties.method="max")/n_perm))
  fdrDT <- tfdr_exp[tfdr_obs,on=.(p)][,.(p,fdr=exp/obs,order=.I)]
  # now fdr needs to be made monotonic
  #setkey(fdrDT,fdr)
  fdrDT[,fdr_m:=fdr_monotonic(fdr)]
  setkey(fdrDT,order)
  return(fdrDT$fdr_m)
}

correct_fdr_permDT <- function(permDT) {
  setkey(permDT,perm_n,enriched_p)
  # emp fdr then BH
  permDT[perm_n==0,enriched_fdr:=correct_fdr_col(p_values = permDT$enriched_p,n_set = permDT[,uniqueN(id)],n_perm = permDT[,max(perm_n)])]
  permDT[perm_n==0,enriched_BH:=p.adjust(enriched_p,method = 'BH')]
  setkey(permDT,perm_n,depleted_p)
  permDT[perm_n==0,depleted_fdr:=correct_fdr_col(p_values = permDT$depleted_p,n_set = permDT[,uniqueN(id)],n_perm = permDT[,max(perm_n)])]
  permDT[perm_n==0,depleted_BH:=p.adjust(depleted_p,method = 'BH')]
  setkey(permDT,perm_n,two_tailed_p)
  permDT[perm_n==0,two_tailed_fdr:=correct_fdr_col(p_values = permDT$two_tailed_p,n_set = permDT[,uniqueN(id)],n_perm = permDT[,max(perm_n)])]
  permDT[perm_n==0,two_tailed_BH:=p.adjust(two_tailed_p,method = 'BH')]
  return(permDT[])
}