Fix outlier detection

R/combine_bgc_output.R

@@ -60,9 +60,19 @@ combine_bgc_output <- function(results.dir,
                    full.names = TRUE)
   check_if_files(hi)
 
+  ta <- list.files(path = results.dir,
+                   pattern = paste0(prefix, ".*ta_\\d+$"),
+                   full.names = TRUE)
+  check_if_files(ta)
+
+  tb <- list.files(path = results.dir,
+                   pattern = paste0(prefix, ".*tb_\\d+$"),
+                   full.names = TRUE)
+  check_if_files(tb)
+
   gc()
 
-  if (!countBGCruns(list(lnl, a0, b0, qa, qb, hi))){
+  if (!countBGCruns(list(lnl, a0, b0, qa, qb, hi, ta, tb))){
     stop("Error: All parameters must have the same number of BGC runs!")
   }
 
@@ -86,14 +96,18 @@ combine_bgc_output <- function(results.dir,
   gc()
   hi.df <- bind_bgc_runs(hi, thin = thin, discard = discard)
   gc()
+  ta.df <- bind_bgc_runs(ta, thin = thin, discard = discard)
+  gc()
+  tb.df <- bind_bgc_runs(tb, thin = thin, discard = discard)
+  gc()
 
   writeLines(paste0("\n\n # MCMC samples = ", ncol(lnl.df)))
 
 
 writeLines(paste0("\n\nDone! Found ", length(lnl), " BGC runs.\n\n"))
 
 df.list <-
-  list(lnl=lnl.df, a0=a0.df, b0=b0.df, qa=qa.df, qb=qb.df, hi=hi.df)
+  list(lnl=lnl.df, a0=a0.df, b0=b0.df, qa=qa.df, qb=qb.df, hi=hi.df, ta=ta.df, tb=tb.df)
 
 return(df.list)
 }

R/get_bgc_outliers.R

@@ -40,7 +40,7 @@ get_cis <- function(df){
 #' @param qn Upper quantile interval value
 #' @return data.frame containing SNP outlier info
 #' @noRd
-get_ab_outliers <- function(a.out, b.out, qa.out, qb.out, loci, qn){
+get_ab_outliers <- function(a.out, b.out, qa.out, qb.out, loci, qn, ta, tb){
   # Get excess ancestry and outliers for alpha and beta parameters.
   names(a.out)[3:4] <- c('lb','ub')
   names(b.out)[3:4] <- c('lb','ub')
@@ -62,8 +62,8 @@ get_ab_outliers <- function(a.out, b.out, qa.out, qb.out, loci, qn){
   # Add in quantiles
   a.out$q <- qa.out$mean
   b.out$q <- qb.out$mean
-
-  # qnorm takes SD, so take the square root of the quantile estimate.
+  
+  # qnorm takes SD, so take square root of the reciprocal of tau.
   # This basically generates the interval to check if alpha and beta are
   # outliers.
   # The interval is defined by qn = n/2 and qn = 1-n/2.
@@ -72,12 +72,13 @@ get_ab_outliers <- function(a.out, b.out, qa.out, qb.out, loci, qn){
   # By default, qn is set to 0.025 and 0.975
   qn_lower <- 1 - qn
 
-  a.out$qlb <- qnorm(qn_lower,0,sqrt(a.out$q))
-  a.out$qub <- qnorm(qn,0,sqrt(a.out$q))
-
-  b.out$qlb <- qnorm(qn_lower,0,sqrt(b.out$q))
-  b.out$qub <- qnorm(qn,0,sqrt(b.out$q))
-
+  a.out$qlb <- qnorm(qn_lower,0,sqrt(1/ta))
+  a.out$qub <- qnorm(qn,0,sqrt(1/ta))
+
+  b.out$qlb <- qnorm(qn_lower,0,sqrt(1/tb))
+  b.out$qub <- qnorm(qn,0,sqrt(1/tb))
+
+
   # Check if median falls outside qn interval.
   a.out$outlier <- NA
   a.out$outlier[a.out$median > a.out$qub] <- 'pos'
@@ -170,6 +171,8 @@ get_bgc_outliers <- function(df.list,
   gc()
   qb <- get_cis(df.list[[5]])
   gc()
+  ta <- mean(unlist(bgc.genes[["ta"]])) # get estimate of Tau-alpha
+  tb <- mean(unlist(bgc.genes[["tb"]])) # get estimate of Tau-beta
 
   clean=FALSE
   #if no loci file provided, make a spoof one here
@@ -189,7 +192,7 @@ get_bgc_outliers <- function(df.list,
   }
 
   # Find SNPs with excess ancestry and outliers SNPs.
-  snps <- get_ab_outliers(a, b, qa, qb, loci.file, qn)
+  snps <- get_ab_outliers(a, b, qa, qb, loci.file, qn, ta, tb)
   gc()
 
   # Get hybrid indexes.

R/plot_traces.R

@@ -50,6 +50,18 @@ plot_traces <- function(df.list, prefix, plotDIR = "./plots", showPLOTS=FALSE){
                 prefix = prefix,
                 plotDIR = plotDIR,
                 showPLOTS=showPLOTS)
+  plotBGCparams(df = df.list[[7]],
+                ylab = "Tau-alpha",
+                bgc_param = "ta",
+                prefix = prefix,
+                plotDIR = plotDIR,
+                showPLOTS=showPLOTS)
+  plotBGCparams(df = df.list[[8]],
+                ylab = "Tau-beta",
+                bgc_param = "tb",
+                prefix = prefix,
+                plotDIR = plotDIR,
+                showPLOTS=showPLOTS)
 
   writeLines(paste("Saved trace plots to", plotDIR))
 

scripts/run_bgc.sh

@@ -199,6 +199,14 @@
 $ESTPOST_PATH/estpost -i ${RESULTS_DIR}/${prefix}_mcmcout_${run}.hdf5 \
 	-p zeta-quantile -o ${RESULTS_DIR}/${prefix}_stat_qb_${run} -s 2 -w 0
 
+# Get tau-alpha (ta)
+$ESTPOST_PATH/estpost -i ${RESULTS_DIR}/${prefix}_mcmcout_${run}.hdf5 \
+	-p tau-alpha -o ${RESULTS_DIR}/${prefix}_stat_ta_${run} -s 2 -w 0
+
+# Get tau-beta (tb)
+$ESTPOST_PATH/estpost -i ${RESULTS_DIR}/${prefix}_mcmcout_${run}.hdf5 \
+	-p tau-beta -o ${RESULTS_DIR}/${prefix}_stat_tb_${run} -s 2 -w 0
+
 echo "Done with estpost!"
 
 exit 0