diff --git a/R/netThresh.R b/R/netThresh.R
index cdb47be..6974a8d 100644
--- a/R/netThresh.R
+++ b/R/netThresh.R
@@ -81,13 +81,15 @@ netThresh <- function(net, asvTab, asvCols = NULL, clusterCol = NULL, cluster =
   clusterColumns <- colnames(clust_ag)
   clust_ag <- cbind(asvTab[, -which(colnames(asvTab) %in% asvCols)], clust_ag)
   #* `calibrate phenotype by calibratePheno`
-  netThreshOut <- do.call(rbind, lapply(phenoCols, function(phenotype) {
-    clust_ag <- clust_ag[!is.na(clust_ag[[phenotype]]), ]
-    if (!is.null(calibratePheno)) {
+  if (!is.null(calibratePheno)) {
+    for (phenotype in phenoCols) {
       formString <- paste0(phenotype, "~", paste0(calibratePheno, collapse = "+"))
-      clust_ag[[phenotype]] <- residuals(lm(as.formula(formString), data = clust_ag))
+      clust_ag[[phenotype]] <- residuals(lm(as.formula(formString),
+                                            data = clust_ag, na.action = na.exclude))
     }
-    thresh_df <- do.call(rbind, parallel::mclapply(clusterColumns, function(col) {
+  }
+  netThreshOut <- do.call(rbind, parallel::mclapply(clusterColumns, function(col) {
+    thresh_df <- do.call(rbind, lapply(phenoCols, function(phenotype) {
       if (model == "hinge" | model == "M01") {
         model <- "hinge"
         f1 <- as.formula(paste0(phenotype, "~1"))
@@ -125,8 +127,8 @@ netThresh <- function(net, asvTab, asvCols = NULL, clusterCol = NULL, cluster =
         warning = function(war) {},
         error = function(err) {}
       )
-    }, mc.cores = cores))
+    }))
     return(thresh_df)
-  }))
+  }, mc.cores = cores))
   return(netThreshOut)
 }
diff --git a/R/thresh.R b/R/thresh.R
index cfbe5d1..3b23670 100644
--- a/R/thresh.R
+++ b/R/thresh.R
@@ -31,15 +31,15 @@ thresh <- function(asvTab, phenoCols, asvCols = NULL, model = "hinge",
   if (is.null(asvCols)) {
     asvCols <- colnames(asvTab)[grepl("ASV", colnames(asvTab))]
   }
-  threshOut <- do.call(rbind, lapply(phenoCols, function(phenotype) {
-    message(paste0("Running ", phenotype))
-    if (!is.null(calibratePheno)) {
+  if (!is.null(calibratePheno)) {
+    for (phenotype in phenoCols) {
       formString <- paste0(phenotype, "~", paste0(calibratePheno, collapse = "+"))
-      asvTab[[phenotype]] <- residuals(lm(as.formula(formString), data = asvTab, na.action = na.exclude))
-    } else {
-      asvTab[[phenotype]] <- asvTab[[phenotype]]
+      asvTab[[phenotype]] <- residuals(lm(as.formula(formString),
+                                          data = asvTab, na.action = na.exclude))
     }
-    thresh_df <- do.call(rbind, parallel::mclapply(asvCols, function(asv_col) {
+  }
+  threshOut <- do.call(rbind, parallel::mclapply(asvCols, function(asv_col) {
+    thresh_df <- do.call(rbind, lapply(phenoCols, function(phenotype) {
       if (model == "hinge" | model == "M01") {
         model <- "hinge"
         f1 <- as.formula(paste0(phenotype, "~1"))
@@ -53,7 +53,6 @@ thresh <- function(asvTab, phenoCols, asvCols = NULL, model = "hinge",
         f1 <- as.formula(paste0(phenotype, "~1"))
         f2 <- as.formula(paste0("~", asv_col))
       }
-
       sub <- asvTab[, c(phenotype, asv_col)]
       tryCatch(
         {
@@ -77,8 +76,8 @@ thresh <- function(asvTab, phenoCols, asvCols = NULL, model = "hinge",
         warning = function(war) {},
         error = function(err) {}
       )
-    }, mc.cores = cores))
+    }))
     thresh_df
-  }))
+  }, mc.cores = cores))
   return(threshOut)
 }