diff --git a/.gitattributes b/.gitattributes
deleted file mode 100644
index 00ed779..0000000
--- a/.gitattributes
+++ /dev/null
@@ -1 +0,0 @@
-*.html	    linguist-vendored=true
\ No newline at end of file
diff --git a/DESCRIPTION b/DESCRIPTION
index f8da4f9..09d0c64 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: DAtest
 Title: Comparing Differential Abundance/Expression Methods
-Version: 2.7.14
+Version: 2.7.15
 Authors@R: person("Jakob", "Russel",
     email = "russel2620@gmail.com", 
     role = c("aut", "cre"))
diff --git a/NAMESPACE b/NAMESPACE
index 87645cf..eb9d2b0 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -41,7 +41,6 @@ export(DA.mva)
 export(DA.neb)
 export(DA.pea)
 export(DA.per)
-export(DA.phyloseq)
 export(DA.poi)
 export(DA.qpo)
 export(DA.qua)
@@ -56,17 +55,13 @@ export(DA.zig)
 export(DA.znb)
 export(DA.zpo)
 export(DA.zzz)
-export(add.tax.DA)
 export(allDA)
 export(featurePlot)
-export(gm_mean)
 export(norm_alr)
 export(norm_clr)
 export(powerDA)
 export(preDA)
-export(pruneTests)
 export(runtimeDA)
-export(spikein)
 export(testDA)
 export(vennDA)
 import(MASS)
diff --git a/R/DA.adx.R b/R/DA.adx.R
index 5e21b5f..1010a61 100644
--- a/R/DA.adx.R
+++ b/R/DA.adx.R
@@ -41,7 +41,7 @@ DA.adx <- function(data, predictor, ...){
     x[!is.na(x$effect) & x$effect < 0,"ordering"] <- paste0(levels(as.factor(predictor))[1],">",levels(as.factor(predictor))[2])
     x$Feature <- rownames(x)
     
-    if(class(data) == "phyloseq") x <- add.tax.DA(data, x)
+    if(class(data) == "phyloseq") x <- addTax(data, x)
     
     return(x)
   } else {
diff --git a/R/DA.aoa.R b/R/DA.aoa.R
index 2fc881b..2c89329 100644
--- a/R/DA.aoa.R
+++ b/R/DA.aoa.R
@@ -69,7 +69,7 @@ DA.aoa <- function(data, predictor, covars = NULL, p.adj = "fdr", delta = 1, all
     res$pval.adj <- p.adjust(res$pval, method = p.adj)
     res$Feature <- rownames(res)
     res$Method <- "ANOVA - ALR (aoa)"
-    if(class(data) == "phyloseq") res <- add.tax.DA(data, res)
+    if(class(data) == "phyloseq") res <- addTax(data, res)
     return(res)
   }
 }
diff --git a/R/DA.aoc.R b/R/DA.aoc.R
index 10b5f38..2d04396 100644
--- a/R/DA.aoc.R
+++ b/R/DA.aoc.R
@@ -67,7 +67,7 @@ DA.aoc <- function(data, predictor, covars = NULL, p.adj = "fdr", delta = 1, all
     res$pval.adj <- p.adjust(res$pval, method = p.adj)
     res$Feature <- rownames(res)
     res$Method <- "ANOVA - CLR (aoc)"
-    if(class(data) == "phyloseq") res <- add.tax.DA(data, res)
+    if(class(data) == "phyloseq") res <- addTax(data, res)
     return(res)
   }
 }
diff --git a/R/DA.aov.R b/R/DA.aov.R
index a2ca696..050a8ae 100644
--- a/R/DA.aov.R
+++ b/R/DA.aov.R
@@ -67,7 +67,7 @@ DA.aov <- function(data, predictor, covars = NULL, relative = TRUE, p.adj = "fdr
     res$pval.adj <- p.adjust(res$pval, method = p.adj)
     res$Feature <- rownames(res)
     res$Method <- "ANOVA (aov)"
-    if(class(data) == "phyloseq") res <- add.tax.DA(data, res)
+    if(class(data) == "phyloseq") res <- addTax(data, res)
     return(res)
   }
 }
diff --git a/R/DA.bay.R b/R/DA.bay.R
index a53561a..4dcb406 100644
--- a/R/DA.bay.R
+++ b/R/DA.bay.R
@@ -57,7 +57,7 @@ DA.bay <- function(data, predictor, allResults = FALSE, ...){
     
     output_df$Method <- "baySeq (bay)"
     
-    if(class(data) == "phyloseq") output_df <- add.tax.DA(data, output_df)
+    if(class(data) == "phyloseq") output_df <- addTax(data, output_df)
     
     if(allResults) return(CD) else return(output_df)
   } else {
diff --git a/R/DA.ds2.R b/R/DA.ds2.R
index 3c93143..26a9ccc 100644
--- a/R/DA.ds2.R
+++ b/R/DA.ds2.R
@@ -126,7 +126,7 @@ DA.ds2 <- function(data, predictor, paired = NULL, covars = NULL, out.all = NULL
     res$Feature <- DESeq2::results(x)@rownames
     res$Method <- "DESeq2 man. geoMeans (ds2)"
     
-    if(class(data) == "phyloseq") res <- add.tax.DA(data, res)
+    if(class(data) == "phyloseq") res <- addTax(data, res)
     
     if(allResults) return(x) else return(res)
   } else {
diff --git a/R/DA.ds2x.R b/R/DA.ds2x.R
index 8ce2316..ab1ace3 100644
--- a/R/DA.ds2x.R
+++ b/R/DA.ds2x.R
@@ -118,7 +118,7 @@ DA.ds2x <- function(data, predictor, paired = NULL, covars = NULL, out.all = NUL
     res$Feature <- DESeq2::results(x)@rownames
     res$Method <- "DESeq2 (ds2x)"
     
-    if(class(data) == "phyloseq") res <- add.tax.DA(data, res)
+    if(class(data) == "phyloseq") res <- addTax(data, res)
     
     if(allResults) return(x) else return(res)  
   } else {
diff --git a/R/DA.ere.R b/R/DA.ere.R
index 0e0fe4f..41e6404 100644
--- a/R/DA.ere.R
+++ b/R/DA.ere.R
@@ -59,7 +59,7 @@ DA.ere <- function(data, predictor, p.adj = "fdr", ...){
     ta$Feature <- rownames(ta)
     ta$Method <- "EdgeR exact - TMM (ere)"
     
-    if(class(data) == "phyloseq") ta <- add.tax.DA(data, ta)
+    if(class(data) == "phyloseq") ta <- addTax(data, ta)
     
     return(ta) 
     
diff --git a/R/DA.ere2.R b/R/DA.ere2.R
index 33027c7..310e44b 100644
--- a/R/DA.ere2.R
+++ b/R/DA.ere2.R
@@ -58,7 +58,7 @@ DA.ere2 <- function(data, predictor, p.adj = "fdr", ...){
     ta$Feature <- rownames(ta)
     ta$Method <- "EdgeR exact - RLE (ere2)"
     
-    if(class(data) == "phyloseq") ta <- add.tax.DA(data, ta)
+    if(class(data) == "phyloseq") ta <- addTax(data, ta)
     
     return(ta) 
   } else {
diff --git a/R/DA.erq.R b/R/DA.erq.R
index 608127a..30ee611 100644
--- a/R/DA.erq.R
+++ b/R/DA.erq.R
@@ -109,7 +109,7 @@ DA.erq <- function(data, predictor, paired = NULL, covars = NULL, out.all = NULL
     ta$Feature <- rownames(ta)
     ta$Method <- "EdgeR qll - TMM (erq)"
     
-    if(class(data) == "phyloseq") ta <- add.tax.DA(data, ta)
+    if(class(data) == "phyloseq") ta <- addTax(data, ta)
     
     if(allResults) return(qlf) else return(ta)
     
diff --git a/R/DA.erq2.R b/R/DA.erq2.R
index 0fd4dc7..40b50b9 100644
--- a/R/DA.erq2.R
+++ b/R/DA.erq2.R
@@ -109,7 +109,7 @@ DA.erq2 <- function(data, predictor, paired = NULL, covars = NULL, out.all = NUL
     ta$Feature <- rownames(ta)
     ta$Method <- "EdgeR qll - RLE (erq2)"
     
-    if(class(data) == "phyloseq") ta <- add.tax.DA(data, ta)
+    if(class(data) == "phyloseq") ta <- addTax(data, ta)
     
     if(allResults) return(qlf) else return(ta)
     
diff --git a/R/DA.fri.R b/R/DA.fri.R
index 9c9dd7f..aeebf0c 100644
--- a/R/DA.fri.R
+++ b/R/DA.fri.R
@@ -62,7 +62,7 @@ DA.fri <- function(data, predictor, paired = NULL, relative = TRUE, p.adj = "fdr
     res$Feature <- gsub(".Friedman.*","",rownames(res))
     res$Method <- "Friedman (fri)" 
     
-    if(class(data) == "phyloseq") res <- add.tax.DA(data, res)
+    if(class(data) == "phyloseq") res <- addTax(data, res)
     
     return(res)
   }
diff --git a/R/DA.kru.R b/R/DA.kru.R
index 42201ab..c86bf3e 100644
--- a/R/DA.kru.R
+++ b/R/DA.kru.R
@@ -54,7 +54,7 @@ DA.kru <- function(data, predictor, relative = TRUE, p.adj = "fdr", allResults =
     res$pval.adj <- p.adjust(res$pval, method = p.adj)
     res$Feature <- rownames(res)
     res$Method <- "Kruskal-Wallis (kru)" 
-    if(class(data) == "phyloseq") res <- add.tax.DA(data, res)
+    if(class(data) == "phyloseq") res <- addTax(data, res)
     return(res)
   }
  
diff --git a/R/DA.lao.R b/R/DA.lao.R
index 8e60251..36604ff 100644
--- a/R/DA.lao.R
+++ b/R/DA.lao.R
@@ -65,7 +65,7 @@ DA.lao <- function(data, predictor, covars = NULL, relative = TRUE, p.adj = "fdr
     res$pval.adj <- p.adjust(res$pval, method = p.adj)
     res$Feature <- rownames(res)
     res$Method <- "Log ANOVA (lao)"
-    if(class(data) == "phyloseq") res <- add.tax.DA(data, res)
+    if(class(data) == "phyloseq") res <- addTax(data, res)
     return(res)
   }
 }
diff --git a/R/DA.lao2.R b/R/DA.lao2.R
index 698e282..a19ae3e 100644
--- a/R/DA.lao2.R
+++ b/R/DA.lao2.R
@@ -64,7 +64,7 @@ DA.lao2 <- function(data, predictor, covars = NULL, p.adj = "fdr", delta = 0.001
     res$pval.adj <- p.adjust(res$pval, method = p.adj)
     res$Feature <- rownames(res)
     res$Method <- "Log ANOVA 2 (lao2)"
-    if(class(data) == "phyloseq") res <- add.tax.DA(data, res)
+    if(class(data) == "phyloseq") res <- addTax(data, res)
     return(res)
   }
 }
diff --git a/R/DA.lia.R b/R/DA.lia.R
index cd08260..b27355d 100644
--- a/R/DA.lia.R
+++ b/R/DA.lia.R
@@ -104,7 +104,7 @@ DA.lia <- function(data, predictor, paired = NULL, covars = NULL, out.all = NULL
     res$Feature <- rownames(res)
     res$Method <- "LIMMA - ALR (lia)"
     
-    if(class(data) == "phyloseq") res <- add.tax.DA(data, res)
+    if(class(data) == "phyloseq") res <- addTax(data, res)
     
     if(allResults) return(fit.eb) else return(res) 
     
diff --git a/R/DA.lic.R b/R/DA.lic.R
index 0cb1447..68c45db 100644
--- a/R/DA.lic.R
+++ b/R/DA.lic.R
@@ -104,7 +104,7 @@ DA.lic <- function(data, predictor, paired = NULL, covars = NULL, out.all = NULL
     res$Feature <- rownames(res)
     res$Method <- "LIMMA - CLR (lic)"
     
-    if(class(data) == "phyloseq") res <- add.tax.DA(data, res)
+    if(class(data) == "phyloseq") res <- addTax(data, res)
     
     if(allResults) return(fit.eb) else return(res) 
     
diff --git a/R/DA.lim.R b/R/DA.lim.R
index 213d2b6..667528a 100644
--- a/R/DA.lim.R
+++ b/R/DA.lim.R
@@ -105,7 +105,7 @@ DA.lim <- function(data, predictor, paired = NULL, covars = NULL, relative = TRU
     res$Feature <- rownames(res)
     res$Method <- "LIMMA (lim)"
     
-    if(class(data) == "phyloseq") res <- add.tax.DA(data, res)
+    if(class(data) == "phyloseq") res <- addTax(data, res)
     
     if(allResults) return(fit.eb) else return(res)
     
diff --git a/R/DA.lli.R b/R/DA.lli.R
index 67f363f..f03414d 100644
--- a/R/DA.lli.R
+++ b/R/DA.lli.R
@@ -103,7 +103,7 @@ DA.lli <- function(data, predictor, paired = NULL, covars = NULL, relative = TRU
     res$Feature <- rownames(res)
     res$Method <- "Log LIMMA (lli)"
     
-    if(class(data) == "phyloseq") res <- add.tax.DA(data, res)
+    if(class(data) == "phyloseq") res <- addTax(data, res)
     
     if(allResults) return(fit.eb) else return(res) 
     
diff --git a/R/DA.lli2.R b/R/DA.lli2.R
index b9d5d64..7ab7acd 100644
--- a/R/DA.lli2.R
+++ b/R/DA.lli2.R
@@ -103,7 +103,7 @@ DA.lli2 <- function(data, predictor, paired = NULL, covars = NULL, out.all = NUL
     res$Feature <- rownames(res)
     res$Method <- "Log LIMMA 2 (lli2)"
     
-    if(class(data) == "phyloseq") res <- add.tax.DA(data, res)
+    if(class(data) == "phyloseq") res <- addTax(data, res)
     
     if(allResults) return(fit.eb) else return(res) 
     
diff --git a/R/DA.llm.R b/R/DA.llm.R
index ddc0948..a5429a0 100644
--- a/R/DA.llm.R
+++ b/R/DA.llm.R
@@ -156,7 +156,7 @@ DA.llm <- function(data, predictor, paired = NULL, covars = NULL, relative = TRU
     res$pval.adj <- p.adjust(res$pval, method = p.adj)
     res$Feature <- rownames(res)
     res$Method <- "Log Linear reg. (llm)"
-    if(class(data) == "phyloseq") res <- add.tax.DA(data, res)
+    if(class(data) == "phyloseq") res <- addTax(data, res)
     return(res)
   } 
   
diff --git a/R/DA.llm2.R b/R/DA.llm2.R
index 04ac672..46d37ef 100644
--- a/R/DA.llm2.R
+++ b/R/DA.llm2.R
@@ -155,7 +155,7 @@ DA.llm2 <- function(data, predictor, paired = NULL, covars = NULL, out.all = NUL
     res$pval.adj <- p.adjust(res$pval, method = p.adj)
     res$Feature <- rownames(res)
     res$Method <- "Log Linear reg. 2 (llm2)"
-    if(class(data) == "phyloseq") res <- add.tax.DA(data, res)
+    if(class(data) == "phyloseq") res <- addTax(data, res)
     return(res)
   }  
 }
diff --git a/R/DA.lma.R b/R/DA.lma.R
index 7d3ccc6..9c532f0 100644
--- a/R/DA.lma.R
+++ b/R/DA.lma.R
@@ -159,7 +159,7 @@ DA.lma <- function(data, predictor, paired = NULL, covars = NULL, out.all = NULL
     res$pval.adj <- p.adjust(res$pval, method = p.adj)
     res$Feature <- rownames(res)
     res$Method <- "Linear model - ALR (lma)"
-    if(class(data) == "phyloseq") res <- add.tax.DA(data, res)
+    if(class(data) == "phyloseq") res <- addTax(data, res)
     return(res)
   } 
   
diff --git a/R/DA.lmc.R b/R/DA.lmc.R
index 48d19db..5fa797b 100644
--- a/R/DA.lmc.R
+++ b/R/DA.lmc.R
@@ -158,7 +158,7 @@ DA.lmc <- function(data, predictor, paired = NULL, covars = NULL, out.all = NULL
     res$pval.adj <- p.adjust(res$pval, method = p.adj)
     res$Feature <- rownames(res)
     res$Method <- "Linear model - CLR (lmc)"
-    if(class(data) == "phyloseq") res <- add.tax.DA(data, res)
+    if(class(data) == "phyloseq") res <- addTax(data, res)
     return(res)
   } 
   
diff --git a/R/DA.lrm.R b/R/DA.lrm.R
index d751a49..a36bf01 100644
--- a/R/DA.lrm.R
+++ b/R/DA.lrm.R
@@ -154,7 +154,7 @@ DA.lrm <- function(data, predictor, paired = NULL, covars = NULL, relative = TRU
     res$pval.adj <- p.adjust(res$pval, method = p.adj)
     res$Feature <- rownames(res)
     res$Method <- "Linear regression (lrm)"
-    if(class(data) == "phyloseq") res <- add.tax.DA(data, res)
+    if(class(data) == "phyloseq") res <- addTax(data, res)
     return(res)
   }
   
diff --git a/R/DA.ltt.R b/R/DA.ltt.R
index 6f60b26..156e293 100644
--- a/R/DA.ltt.R
+++ b/R/DA.ltt.R
@@ -82,7 +82,7 @@ DA.ltt <- function(data, predictor, paired = NULL, relative = TRUE, p.adj = "fdr
     res[!is.na(res$log2FC) & res$log2FC < 0,"ordering"] <- paste0(levels(as.factor(predictor))[1],">",levels(as.factor(predictor))[2])
     res$Feature <- rownames(res)
     res$Method <- "Log t-test (ltt)"
-    if(class(data) == "phyloseq") res <- add.tax.DA(data, res)
+    if(class(data) == "phyloseq") res <- addTax(data, res)
     return(res)
   }
 }
diff --git a/R/DA.ltt2.R b/R/DA.ltt2.R
index 87d72ff..bb8f53d 100644
--- a/R/DA.ltt2.R
+++ b/R/DA.ltt2.R
@@ -81,7 +81,7 @@ DA.ltt2 <- function(data, predictor, paired = NULL, p.adj = "fdr", delta = 0.001
     res[!is.na(res$log2FC) & res$log2FC < 0,"ordering"] <- paste0(levels(as.factor(predictor))[1],">",levels(as.factor(predictor))[2])
     res$Feature <- rownames(res)
     res$Method <- "Log t-test2 (ltt2)"
-    if(class(data) == "phyloseq") res <- add.tax.DA(data, res)
+    if(class(data) == "phyloseq") res <- addTax(data, res)
     return(res)
   }
 }
diff --git a/R/DA.msf.R b/R/DA.msf.R
index 0436536..50b004e 100644
--- a/R/DA.msf.R
+++ b/R/DA.msf.R
@@ -61,7 +61,7 @@ DA.msf <- function(data, predictor, p.adj = "fdr", allResults = FALSE, ...){
     temp_table[!is.na(temp_table$logFC) & temp_table$logFC < 0,"ordering"] <- paste0(levels(as.factor(predictor))[1],">",levels(as.factor(predictor))[2])
     temp_table$Method <- "MgSeq Feature (msf)"  
     
-    if(class(data) == "phyloseq") temp_table <- add.tax.DA(data, temp_table)
+    if(class(data) == "phyloseq") temp_table <- addTax(data, temp_table)
     
     if(allResults) return(mgsfit) else return(temp_table)
     
diff --git a/R/DA.mva.R b/R/DA.mva.R
index 99fb6bd..42a09ac 100644
--- a/R/DA.mva.R
+++ b/R/DA.mva.R
@@ -105,7 +105,7 @@ DA.mva <- function(data, predictor, paired = NULL, covars = NULL, relative = TRU
       res$Feature <- rownames(pvals.unadj$uni.test)
       res$Method <- "mvabund (mva)"
       
-      if(class(data) == "phyloseq") res <- add.tax.DA(data, res)
+      if(class(data) == "phyloseq") res <- addTax(data, res)
       
       return(res)
     }
diff --git a/R/DA.neb.R b/R/DA.neb.R
index 21394df..5f98dde 100644
--- a/R/DA.neb.R
+++ b/R/DA.neb.R
@@ -211,7 +211,7 @@ DA.neb <- function(data, predictor, paired = NULL, covars = NULL, relative = TRU
     res$Method <- "Negbinom GLM (neb)"
     
     if(nrow(res) > 1){
-      if(class(data) == "phyloseq") res <- add.tax.DA(data, res)
+      if(class(data) == "phyloseq") res <- addTax(data, res)
     }
     return(res)
   } else {
diff --git a/R/DA.pea.R b/R/DA.pea.R
index fad496e..56f310c 100644
--- a/R/DA.pea.R
+++ b/R/DA.pea.R
@@ -54,7 +54,7 @@ DA.pea <- function(data, predictor, relative = TRUE, p.adj = "fdr", ...){
   res$Method <- "Pearson (pea)"
 
   # Add tax table if available
-  if(class(data) == "phyloseq") res <- add.tax.DA(data, res)
+  if(class(data) == "phyloseq") res <- addTax(data, res)
 
   return(res)  
 }
diff --git a/R/DA.per.R b/R/DA.per.R
index 1cdb8ac..d3c908b 100644
--- a/R/DA.per.R
+++ b/R/DA.per.R
@@ -125,7 +125,7 @@ DA.per <- function(data, predictor, paired = NULL, relative = TRUE, p.adj = "fdr
   output_df$pval.adj <- p.adjust(output_df$pval, method = p.adj)
   output_df$Method <- "Permutation (per)"
   
-  if(class(data) == "phyloseq") output_df <- add.tax.DA(data, output_df)
+  if(class(data) == "phyloseq") output_df <- addTax(data, output_df)
   
   return(output_df)
 }
diff --git a/R/DA.phyloseq.R b/R/DA.phyloseq.R
index c9d4062..622bb8e 100644
--- a/R/DA.phyloseq.R
+++ b/R/DA.phyloseq.R
@@ -1,14 +1,3 @@
-#' Extract data from a \code{phyloseq} object to be used in \code{DAtest}
-#' 
-#' Internal function
-#' 
-#' @param data A \code{phyloseq} object
-#' @param predictor The \code{predictor} of interest
-#' @param paired Factor for paired/blocked experimental designs.
-#' @param covars A character vector with names of the variables in \code{sample_data(data)}
-#' @return A list with data extracted from the \code{phyloseq} object
-#' @export
-
 DA.phyloseq <- function(data, predictor, paired = NULL, covars = NULL){
  
   ok <- tryCatch({
diff --git a/R/DA.poi.R b/R/DA.poi.R
index 72f4767..d0c078c 100644
--- a/R/DA.poi.R
+++ b/R/DA.poi.R
@@ -245,7 +245,7 @@ DA.poi <- function(data, predictor, paired = NULL, covars = NULL, relative = TRU
     res$Method <- "Poisson GLM (poi)"
     
     if(nrow(res) > 1){
-      if(class(data) == "phyloseq") res <- add.tax.DA(data, res)
+      if(class(data) == "phyloseq") res <- addTax(data, res)
     }
     return(res)
   }
diff --git a/R/DA.qpo.R b/R/DA.qpo.R
index 7054e6e..6f2c145 100644
--- a/R/DA.qpo.R
+++ b/R/DA.qpo.R
@@ -158,7 +158,7 @@ DA.qpo <- function(data, predictor, covars = NULL, relative = TRUE, out.all = NU
     res$Method <- "Quasi-Poisson GLM (qpo)"
     
     if(nrow(res) > 1){
-      if(class(data) == "phyloseq") res <- add.tax.DA(data, res)
+      if(class(data) == "phyloseq") res <- addTax(data, res)
     }
     return(res)
   }
diff --git a/R/DA.qua.R b/R/DA.qua.R
index 944a306..de1c2a6 100644
--- a/R/DA.qua.R
+++ b/R/DA.qua.R
@@ -64,7 +64,7 @@ DA.qua <- function(data, predictor, paired = NULL, relative = TRUE, p.adj = "fdr
     res$Feature <- gsub(".Quade.*","",rownames(res))
     res$Method <- "Quade (qua)" 
     
-    if(class(data) == "phyloseq") res <- add.tax.DA(data, res)
+    if(class(data) == "phyloseq") res <- addTax(data, res)
     return(res)
   }
   
diff --git a/R/DA.sam.R b/R/DA.sam.R
index ba5f7c5..eac5021 100644
--- a/R/DA.sam.R
+++ b/R/DA.sam.R
@@ -99,7 +99,7 @@ DA.sam <- function(data, predictor, paired = NULL, fdr.output = 0.05, allResults
     
     df$Method <- "SAMseq (sam)"
     
-    if(class(data) == "phyloseq") df <- add.tax.DA(data, df)
+    if(class(data) == "phyloseq") df <- addTax(data, df)
     
     if(allResults){
       return(res)
diff --git a/R/DA.spe.R b/R/DA.spe.R
index e1e0137..8ca2bca 100644
--- a/R/DA.spe.R
+++ b/R/DA.spe.R
@@ -51,7 +51,7 @@ DA.spe <- function(data, predictor, relative = TRUE, p.adj = "fdr", ...){
   res$Feature <- rownames(res)
   res$Method <- "Spearman (spe)"
   
-  if(class(data) == "phyloseq") res <- add.tax.DA(data, res)
+  if(class(data) == "phyloseq") res <- addTax(data, res)
 
   return(res)  
 }
diff --git a/R/DA.tta.R b/R/DA.tta.R
index 4a25efc..dca71c4 100644
--- a/R/DA.tta.R
+++ b/R/DA.tta.R
@@ -86,7 +86,7 @@ DA.tta <- function(data, predictor, paired = NULL, p.adj = "fdr", delta = 1, tes
     res[!is.na(res$log2FC) & res$log2FC < 0,"ordering"] <- paste0(levels(as.factor(predictor))[1],">",levels(as.factor(predictor))[2])
     res$Feature <- rownames(res)
     res$Method <- "t-test - ALR (tta)"
-    if(class(data) == "phyloseq") res <- add.tax.DA(data, res)
+    if(class(data) == "phyloseq") res <- addTax(data, res)
     return(res)
   }
 }
diff --git a/R/DA.ttc.R b/R/DA.ttc.R
index 39d6ca5..8114265 100644
--- a/R/DA.ttc.R
+++ b/R/DA.ttc.R
@@ -84,7 +84,7 @@ DA.ttc <- function(data, predictor, paired = NULL, p.adj = "fdr", delta = 1, tes
     res[!is.na(res$log2FC) & res$log2FC < 0,"ordering"] <- paste0(levels(as.factor(predictor))[1],">",levels(as.factor(predictor))[2])
     res$Feature <- rownames(res)
     res$Method <- "t-test - CLR (ttc)"
-    if(class(data) == "phyloseq") res <- add.tax.DA(data, res)
+    if(class(data) == "phyloseq") res <- addTax(data, res)
     return(res)
   }
 }
diff --git a/R/DA.ttt.R b/R/DA.ttt.R
index e925091..1b8a45b 100644
--- a/R/DA.ttt.R
+++ b/R/DA.ttt.R
@@ -84,7 +84,7 @@ DA.ttt <- function(data, predictor,paired = NULL, relative = TRUE, p.adj = "fdr"
     res[!is.na(res$log2FC) & res$log2FC < 0,"ordering"] <- paste0(levels(as.factor(predictor))[1],">",levels(as.factor(predictor))[2])
     res$Feature <- rownames(res)
     res$Method <- "t-test (ttt)"
-    if(class(data) == "phyloseq") res <- add.tax.DA(data, res)
+    if(class(data) == "phyloseq") res <- addTax(data, res)
     return(res)
   }
 }
\ No newline at end of file
diff --git a/R/DA.vli.R b/R/DA.vli.R
index 72f1e30..0c021a1 100644
--- a/R/DA.vli.R
+++ b/R/DA.vli.R
@@ -108,7 +108,7 @@ DA.vli <- function(data, predictor, paired = NULL, covars = NULL, out.all = NULL
     res$Feature <- rownames(res)
     res$Method <- "LIMMA voom (vli)"
     
-    if(class(data) == "phyloseq") res <- add.tax.DA(data, res)
+    if(class(data) == "phyloseq") res <- addTax(data, res)
     
     if(allResults) return(fit.eb) else return(res)  
     
diff --git a/R/DA.wil.R b/R/DA.wil.R
index 9832f58..2335230 100644
--- a/R/DA.wil.R
+++ b/R/DA.wil.R
@@ -84,7 +84,7 @@ DA.wil <- function(data, predictor, paired = NULL, relative = TRUE, p.adj = "fdr
     res[!is.na(res$log2FC) & res$log2FC < 0,"ordering"] <- paste0(levels(as.factor(predictor))[1],">",levels(as.factor(predictor))[2])
     res$Feature <- rownames(res)
     res$Method <- "Wilcox (wil)" 
-    if(class(data) == "phyloseq") res <- add.tax.DA(data, res)
+    if(class(data) == "phyloseq") res <- addTax(data, res)
     return(res)
   }
  
diff --git a/R/DA.zig.R b/R/DA.zig.R
index 4da9a50..c3a3bc4 100644
--- a/R/DA.zig.R
+++ b/R/DA.zig.R
@@ -94,7 +94,7 @@ DA.zig <- function(data, predictor, paired = NULL, covars = NULL, p.adj = "fdr",
     temp_table$Feature <- rownames(temp_table)
     temp_table$Method <- "MgSeq ZIG (zig)"
     
-    if(class(data) == "phyloseq") temp_table <- add.tax.DA(data, temp_table)
+    if(class(data) == "phyloseq") temp_table <- addTax(data, temp_table)
     
     if(allResults) return(mgsfit) else return(temp_table)
     
diff --git a/R/DA.znb.R b/R/DA.znb.R
index a1d5758..d5cf197 100644
--- a/R/DA.znb.R
+++ b/R/DA.znb.R
@@ -167,7 +167,7 @@ DA.znb <- function(data, predictor, covars = NULL, relative = TRUE, out.all = NU
       res$Method <- "ZI-NegBin GLM (znb)"
       
       if(nrow(res) > 1){
-        if(class(data) == "phyloseq") res <- add.tax.DA(data, res)
+        if(class(data) == "phyloseq") res <- addTax(data, res)
       }
       return(res)
     }
diff --git a/R/DA.zpo.R b/R/DA.zpo.R
index 146dd24..e8c21bf 100644
--- a/R/DA.zpo.R
+++ b/R/DA.zpo.R
@@ -166,7 +166,7 @@ DA.zpo <- function(data, predictor, covars = NULL, relative = TRUE, out.all = NU
       res$Method <- "ZI-Poisson GLM (zpo)"
       
       if(nrow(res) > 1){
-        if(class(data) == "phyloseq") res <- add.tax.DA(data, res)
+        if(class(data) == "phyloseq") res <- addTax(data, res)
       }
       return(res)
     }
diff --git a/R/DA.zzz.R b/R/DA.zzz.R
index 6927e1b..55d78c1 100644
--- a/R/DA.zzz.R
+++ b/R/DA.zzz.R
@@ -69,7 +69,7 @@ DA.zzz <- function(data, predictor, paired = NULL, covars = NULL, p.adj = "fdr",
   if(!all(c("pval","Feature","Method") %in% colnames(res))) stop("The following columns has to be present in output: pval, Feature and Method")
   res$pval.adj <- p.adjust(res$pval, method = p.adj)
     
-  if(class(data) == "phyloseq") res <- add.tax.DA(data, res)
+  if(class(data) == "phyloseq") res <- addTax(data, res)
   
   return(res)
 }
diff --git a/R/add.tax.DA.R b/R/add.tax.DA.R
index 320a5d5..909f6c5 100644
--- a/R/add.tax.DA.R
+++ b/R/add.tax.DA.R
@@ -1,11 +1,4 @@
-#' Add tax_table from phyloseq object
-#'
-#' Internal function
-#' @param data phyloseq object
-#' @param res data.frame with results
-#' @return \code{res} merged with the tax.table from the phyloseq object
-#' @export
-add.tax.DA <- function(data, res){
+addTax <- function(data, res){
   
   loadNamespace("phyloseq")
   
diff --git a/R/allDA.R b/R/allDA.R
index fa825fe..18ec2c4 100644
--- a/R/allDA.R
+++ b/R/allDA.R
@@ -284,7 +284,7 @@ allDA <- function(data, predictor, paired = NULL, covars = NULL,
   if(length(Pval.raw) > 0){
     df.raw <- suppressWarnings(Reduce(function(x,y) merge(x, y, by= "Feature", all.x = TRUE, all.y = TRUE), Pval.raw))
     colnames(df.raw)[2:ncol(df.raw)] <- names(Pval.raw)
-    if(class(data) == "phyloseq") df.raw <- add.tax.DA(data, df.raw)
+    if(class(data) == "phyloseq") df.raw <- addTax(data, df.raw)
   } else {
     df.raw <- NULL
   }
@@ -307,7 +307,7 @@ allDA <- function(data, predictor, paired = NULL, covars = NULL,
       }
       colnames(df.adj)[ncol(df.adj)] <- "sam"
     }
-    if(class(data) == "phyloseq") df.adj <- add.tax.DA(data, df.adj)
+    if(class(data) == "phyloseq") df.adj <- addTax(data, df.adj)
   } else {
     df.adj <- NULL
   }
@@ -364,7 +364,7 @@ allDA <- function(data, predictor, paired = NULL, covars = NULL,
     list.est <- list.est[!sapply(list.est, is.null)]
     if(length(list.est) > 0){
       df.est <- Reduce(function(x,y) merge(x, y, by= "Feature", all.x = TRUE, all.y = TRUE), list.est)
-      if(class(data) == "phyloseq") df.est <- add.tax.DA(data, df.est)
+      if(class(data) == "phyloseq") df.est <- addTax(data, df.est)
     } else {
       df.est <- NULL
     } 
diff --git a/R/norm.R b/R/norm.R
index bdbefa5..083ea75 100644
--- a/R/norm.R
+++ b/R/norm.R
@@ -1,9 +1,3 @@
-#' Geometric means
-#'
-#' @param x numeric vector
-#' @return The geometric mean
-#' @examples gm_mean(1:10)
-#' @export
 gm_mean = function(x){
   if(any(x < 0, na.rm = TRUE)){
     stop("Negative values not allowed")
diff --git a/R/pruneTests.R b/R/pruneTests.R
index 16ba65d..2ec784d 100644
--- a/R/pruneTests.R
+++ b/R/pruneTests.R
@@ -1,15 +1,3 @@
-#' Prune tests argument for \code{testDA}
-#'
-#' @param tests A character vector with names of tests
-#' @param predictor Vector with predictor variable
-#' @param paired Factor with paired variable
-#' @param covars Named list with covariables
-#' @param relative Include tests that work with relative abundances (TRUE) or only absolute abundances (FALSE)
-#' @param decimal Exclude tests that do not work with decimals (TRUE)
-#' @param zeroes If FALSE will exclude tests that include zero-inflation
-#' @return A character with the tests argument pruned according to the parameters
-#' @export
-
 pruneTests <- function(tests, predictor, paired, covars, relative, decimal, zeroes){
 
   # Prune test argument if packages are not installed
@@ -19,8 +7,7 @@ pruneTests <- function(tests, predictor, paired, covars, relative, decimal, zero
   if(!"metagenomeSeq" %in% rownames(installed.packages())) tests <- tests[!tests %in% c("msf","zig")]
   if(!"DESeq2" %in% rownames(installed.packages())) tests <- tests[!tests %in% c("ds2","ds2x")]
   if(!"limma" %in% rownames(installed.packages())) tests <- tests[!tests %in% c("lim","lli","lli2","vli","lia","lic")]
-  if(!"statmod" %in% rownames(installed.packages())) tests <- tests[!tests %in% c("lim","lli","lli2","vli","lia","lic")]
-  if(!"RAIDA" %in% rownames(installed.packages())) tests <- tests[tests != "rai"]
+  if(!"statmod" %in% rownames(installed.packages()) && is.null(paired)) tests <- tests[!tests %in% c("lim","lli","lli2","vli","lia","lic")]
   if(!"pscl" %in% rownames(installed.packages())) tests <- tests[!tests %in% c("zpo","znb")]
   if(!"samr" %in% rownames(installed.packages())) tests <- tests[!tests %in% c("sam")]
   if(!"mvabund" %in% rownames(installed.packages())) tests <- tests[!tests %in% c("mva")]
diff --git a/R/spikein.R b/R/spikein.R
index ca9b4a7..c36ac27 100644
--- a/R/spikein.R
+++ b/R/spikein.R
@@ -1,19 +1,3 @@
-#' Spike-in
-#'
-#' Internal function for the \code{testDA} function.
-#' 
-#' Modified version of the one from:
-#' https://microbiomejournal.biomedcentral.com/articles/10.1186/s40168-016-0208-8.
-#' 
-#' @param count_table Matrix or data.frame. Table with taxa/genes/proteins as rows and samples as columns
-#' @param predictor Factor or Numeric. The predictor of interest. E.g. case and control. If the \code{predictor} has more than two levels, only the 2. level will be spiked.
-#' @param effectSize Integer. The effect size for the spike-ins. Default 2
-#' @param k Vector of length 3. Number of Features to spike in each tertile (lower, mid, upper). k=c(5,10,15): 5 features spiked in low abundance tertile, 10 features spiked in mid abundance tertile and 15 features spiked in high abundance tertile. Default c(5,5,5)
-#' @param num.pred Logical. Is the \code{predictor} numeric? Default FALSE
-#' @param relative Logical. Are abundances relative? Default TRUE
-#' @return Returns a list where the first element is the count_table with the spiked features, and the second element contains rownames of the spiked features.
-#' @export
-
 spikein <- function(count_table, predictor, effectSize = 2, k, num.pred = FALSE, relative = TRUE){
   
   if(effectSize < 0) stop("Effect size should be positive")
diff --git a/R/zzz.R b/R/zzz.R
index fe4ff23..f6b9958 100644
--- a/R/zzz.R
+++ b/R/zzz.R
@@ -1,3 +1,3 @@
 .onAttach <- function(libname, pkgname){
-  packageStartupMessage("DAtest version 2.7.14")
+  packageStartupMessage("DAtest version 2.7.15")
 }
\ No newline at end of file
diff --git a/man/DA.phyloseq.Rd b/man/DA.phyloseq.Rd
deleted file mode 100644
index 579398a..0000000
--- a/man/DA.phyloseq.Rd
+++ /dev/null
@@ -1,23 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/DA.phyloseq.R
-\name{DA.phyloseq}
-\alias{DA.phyloseq}
-\title{Extract data from a \code{phyloseq} object to be used in \code{DAtest}}
-\usage{
-DA.phyloseq(data, predictor, paired = NULL, covars = NULL)
-}
-\arguments{
-\item{data}{A \code{phyloseq} object}
-
-\item{predictor}{The \code{predictor} of interest}
-
-\item{paired}{Factor for paired/blocked experimental designs.}
-
-\item{covars}{A character vector with names of the variables in \code{sample_data(data)}}
-}
-\value{
-A list with data extracted from the \code{phyloseq} object
-}
-\description{
-Internal function
-}
diff --git a/man/add.tax.DA.Rd b/man/add.tax.DA.Rd
deleted file mode 100644
index a2ff32b..0000000
--- a/man/add.tax.DA.Rd
+++ /dev/null
@@ -1,19 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/add.tax.DA.R
-\name{add.tax.DA}
-\alias{add.tax.DA}
-\title{Add tax_table from phyloseq object}
-\usage{
-add.tax.DA(data, res)
-}
-\arguments{
-\item{data}{phyloseq object}
-
-\item{res}{data.frame with results}
-}
-\value{
-\code{res} merged with the tax.table from the phyloseq object
-}
-\description{
-Internal function
-}
diff --git a/man/gm_mean.Rd b/man/gm_mean.Rd
deleted file mode 100644
index c9c1c06..0000000
--- a/man/gm_mean.Rd
+++ /dev/null
@@ -1,20 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/norm.R
-\name{gm_mean}
-\alias{gm_mean}
-\title{Geometric means}
-\usage{
-gm_mean(x)
-}
-\arguments{
-\item{x}{numeric vector}
-}
-\value{
-The geometric mean
-}
-\description{
-Geometric means
-}
-\examples{
-gm_mean(1:10)
-}
diff --git a/man/pruneTests.Rd b/man/pruneTests.Rd
deleted file mode 100644
index 5cca0bd..0000000
--- a/man/pruneTests.Rd
+++ /dev/null
@@ -1,29 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/pruneTests.R
-\name{pruneTests}
-\alias{pruneTests}
-\title{Prune tests argument for \code{testDA}}
-\usage{
-pruneTests(tests, predictor, paired, covars, relative, decimal, zeroes)
-}
-\arguments{
-\item{tests}{A character vector with names of tests}
-
-\item{predictor}{Vector with predictor variable}
-
-\item{paired}{Factor with paired variable}
-
-\item{covars}{Named list with covariables}
-
-\item{relative}{Include tests that work with relative abundances (TRUE) or only absolute abundances (FALSE)}
-
-\item{decimal}{Exclude tests that do not work with decimals (TRUE)}
-
-\item{zeroes}{If FALSE will exclude tests that include zero-inflation}
-}
-\value{
-A character with the tests argument pruned according to the parameters
-}
-\description{
-Prune tests argument for \code{testDA}
-}
diff --git a/man/spikein.Rd b/man/spikein.Rd
deleted file mode 100644
index 499b8cb..0000000
--- a/man/spikein.Rd
+++ /dev/null
@@ -1,32 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/spikein.R
-\name{spikein}
-\alias{spikein}
-\title{Spike-in}
-\usage{
-spikein(count_table, predictor, effectSize = 2, k, num.pred = FALSE,
-  relative = TRUE)
-}
-\arguments{
-\item{count_table}{Matrix or data.frame. Table with taxa/genes/proteins as rows and samples as columns}
-
-\item{predictor}{Factor or Numeric. The predictor of interest. E.g. case and control. If the \code{predictor} has more than two levels, only the 2. level will be spiked.}
-
-\item{effectSize}{Integer. The effect size for the spike-ins. Default 2}
-
-\item{k}{Vector of length 3. Number of Features to spike in each tertile (lower, mid, upper). k=c(5,10,15): 5 features spiked in low abundance tertile, 10 features spiked in mid abundance tertile and 15 features spiked in high abundance tertile. Default c(5,5,5)}
-
-\item{num.pred}{Logical. Is the \code{predictor} numeric? Default FALSE}
-
-\item{relative}{Logical. Are abundances relative? Default TRUE}
-}
-\value{
-Returns a list where the first element is the count_table with the spiked features, and the second element contains rownames of the spiked features.
-}
-\description{
-Internal function for the \code{testDA} function.
-}
-\details{
-Modified version of the one from:
-https://microbiomejournal.biomedcentral.com/articles/10.1186/s40168-016-0208-8.
-}