diff --git a/R/cTMed-boot-med-dot.R b/R/cTMed-boot-med-dot.R
index bc76c47..5ba7204 100644
--- a/R/cTMed-boot-med-dot.R
+++ b/R/cTMed-boot-med-dot.R
@@ -9,58 +9,116 @@
   par <- FALSE
   if (!is.null(ncores)) {
     ncores <- as.integer(ncores)
+    R <- length(phi)
+    if (ncores > R) {
+      ncores <- R
+    }
     if (ncores > 1) {
       par <- TRUE
     }
   }
   if (par) {
-    cl <- parallel::makeCluster(ncores)
-    on.exit(
-      parallel::stopCluster(cl = cl)
-    )
-    output <- lapply(
-      X = delta_t,
-      FUN = function(i) {
-        thetahatstar <- parallel::parLapply(
-          cl = cl,
-          X = phi,
-          fun = .Med,
-          delta_t = i,
-          from = from,
-          to = to,
-          med = med
-        )
-        thetahatstar <- do.call(
-          what = "rbind",
-          args = thetahatstar
-        )
-        colnames(thetahatstar) <- c(
-          "total",
-          "direct",
-          "indirect",
-          "interval"
-        )
-        est <- .Med(
-          phi = phi_hat,
-          delta_t = i,
-          from = from,
-          to = to,
-          med = med
-        )
-        names(est) <- c(
-          "total",
-          "direct",
-          "indirect",
-          "interval"
-        )
-        out <- list(
-          delta_t = i,
-          est = est,
-          thetahatstar = thetahatstar
-        )
-        return(out)
-      }
-    )
+    os_type <- Sys.info()["sysname"]
+    if (os_type == "Darwin") {
+      fork <- TRUE
+    } else if (os_type == "Linux") {
+      fork <- TRUE
+    } else {
+      fork <- FALSE
+    }
+    if (fork) {
+      output <- lapply(
+        X = delta_t,
+        FUN = function(i) {
+          thetahatstar <- parallel::mclapply(
+            X = phi,
+            FUN = .Med,
+            delta_t = i,
+            from = from,
+            to = to,
+            med = med,
+            mc.cores = ncores
+          )
+          thetahatstar <- do.call(
+            what = "rbind",
+            args = thetahatstar
+          )
+          colnames(thetahatstar) <- c(
+            "total",
+            "direct",
+            "indirect",
+            "interval"
+          )
+          est <- .Med(
+            phi = phi_hat,
+            delta_t = i,
+            from = from,
+            to = to,
+            med = med
+          )
+          names(est) <- c(
+            "total",
+            "direct",
+            "indirect",
+            "interval"
+          )
+          out <- list(
+            delta_t = i,
+            est = est,
+            thetahatstar = thetahatstar
+          )
+          return(out)
+        }
+      )
+    } else {
+      cl <- parallel::makeCluster(ncores)
+      on.exit(
+        parallel::stopCluster(cl = cl)
+      )
+      output <- lapply(
+        X = delta_t,
+        FUN = function(i) {
+          thetahatstar <- parallel::parLapply(
+            cl = cl,
+            X = phi,
+            fun = .Med,
+            delta_t = i,
+            from = from,
+            to = to,
+            med = med
+          )
+          thetahatstar <- do.call(
+            what = "rbind",
+            args = thetahatstar
+          )
+          colnames(thetahatstar) <- c(
+            "total",
+            "direct",
+            "indirect",
+            "interval"
+          )
+          est <- .Med(
+            phi = phi_hat,
+            delta_t = i,
+            from = from,
+            to = to,
+            med = med
+          )
+          names(est) <- c(
+            "total",
+            "direct",
+            "indirect",
+            "interval"
+          )
+          out <- list(
+            delta_t = i,
+            est = est,
+            thetahatstar = thetahatstar
+          )
+          return(out)
+        }
+      )
+    }
     # nocov end
   } else {
     output <- lapply(
diff --git a/R/cTMed-delta-beta-std.R b/R/cTMed-delta-beta-std.R
index 9eba052..e63acd9 100644
--- a/R/cTMed-delta-beta-std.R
+++ b/R/cTMed-delta-beta-std.R
@@ -260,6 +260,10 @@ DeltaBetaStd <- function(phi,
   par <- FALSE
   if (!is.null(ncores)) {
     ncores <- as.integer(ncores)
+    R <- length(delta_t)
+    if (ncores > R) {
+      ncores <- R
+    }
     if (ncores > 1) {
       par <- TRUE
     }
diff --git a/R/cTMed-delta-beta.R b/R/cTMed-delta-beta.R
index 85124ab..a7a7a10 100644
--- a/R/cTMed-delta-beta.R
+++ b/R/cTMed-delta-beta.R
@@ -229,22 +229,44 @@ DeltaBeta <- function(phi,
   par <- FALSE
   if (!is.null(ncores)) {
     ncores <- as.integer(ncores)
+    R <- length(delta_t)
+    if (ncores > R) {
+      ncores <- R
+    }
     if (ncores > 1) {
       par <- TRUE
     }
   }
   if (par) {
-    cl <- parallel::makeCluster(ncores)
-    on.exit(
-      parallel::stopCluster(cl = cl)
-    )
-    output <- parallel::parLapply(
-      cl = cl,
-      X = delta_t,
-      fun = .DeltaBeta,
-      phi = phi,
-      vcov_phi_vec = vcov_phi_vec
-    )
+    os_type <- Sys.info()["sysname"]
+    if (os_type == "Darwin") {
+      fork <- TRUE
+    } else if (os_type == "Linux") {
+      fork <- TRUE
+    } else {
+      fork <- FALSE
+    }
+    if (fork) {
+      output <- parallel::mclapply(
+        X = delta_t,
+        FUN = .DeltaBeta,
+        phi = phi,
+        vcov_phi_vec = vcov_phi_vec,
+        mc.cores = ncores
+      )
+    } else {
+      cl <- parallel::makeCluster(ncores)
+      on.exit(
+        parallel::stopCluster(cl = cl)
+      )
+      output <- parallel::parLapply(
+        cl = cl,
+        X = delta_t,
+        fun = .DeltaBeta,
+        phi = phi,
+        vcov_phi_vec = vcov_phi_vec
+      )
+    }
     # nocov end
   } else {
     output <- lapply(
diff --git a/R/cTMed-delta-indirect-central.R b/R/cTMed-delta-indirect-central.R
index e2e4de5..5082305 100644
--- a/R/cTMed-delta-indirect-central.R
+++ b/R/cTMed-delta-indirect-central.R
@@ -226,23 +226,46 @@ DeltaIndirectCentral <- function(phi,
   par <- FALSE
   if (!is.null(ncores)) {
     ncores <- as.integer(ncores)
+    R <- length(delta_t)
+    if (ncores > R) {
+      ncores <- R
+    }
     if (ncores > 1) {
       par <- TRUE
     }
   }
   if (par) {
-    cl <- parallel::makeCluster(ncores)
-    on.exit(
-      parallel::stopCluster(cl = cl)
-    )
-    output <- parallel::parLapply(
-      cl = cl,
-      X = delta_t,
-      fun = .DeltaCentral,
-      phi = phi,
-      vcov_phi_vec = vcov_phi_vec,
-      total = total
-    )
+    os_type <- Sys.info()["sysname"]
+    if (os_type == "Darwin") {
+      fork <- TRUE
+    } else if (os_type == "Linux") {
+      fork <- TRUE
+    } else {
+      fork <- FALSE
+    }
+    if (fork) {
+      output <- parallel::mclapply(
+        X = delta_t,
+        FUN = .DeltaCentral,
+        phi = phi,
+        vcov_phi_vec = vcov_phi_vec,
+        total = total,
+        mc.cores = ncores
+      )
+    } else {
+      cl <- parallel::makeCluster(ncores)
+      on.exit(
+        parallel::stopCluster(cl = cl)
+      )
+      output <- parallel::parLapply(
+        cl = cl,
+        X = delta_t,
+        fun = .DeltaCentral,
+        phi = phi,
+        vcov_phi_vec = vcov_phi_vec,
+        total = total
+      )
+    }
     # nocov end
   } else {
     output <- lapply(
diff --git a/R/cTMed-delta-med-std.R b/R/cTMed-delta-med-std.R
index 33b9305..45563f9 100644
--- a/R/cTMed-delta-med-std.R
+++ b/R/cTMed-delta-med-std.R
@@ -298,26 +298,52 @@ DeltaMedStd <- function(phi,
   par <- FALSE
   if (!is.null(ncores)) {
     ncores <- as.integer(ncores)
+    R <- length(delta_t)
+    if (ncores > R) {
+      ncores <- R
+    }
     if (ncores > 1) {
       par <- TRUE
     }
   }
   if (par) {
-    cl <- parallel::makeCluster(ncores)
-    on.exit(
-      parallel::stopCluster(cl = cl)
-    )
-    output <- parallel::parLapply(
-      cl = cl,
-      X = delta_t,
-      fun = .DeltaMedStd,
-      phi = phi,
-      sigma = sigma,
-      vcov_theta = vcov_theta,
-      from = from,
-      to = to,
-      med = med
-    )
+    os_type <- Sys.info()["sysname"]
+    if (os_type == "Darwin") {
+      fork <- TRUE
+    } else if (os_type == "Linux") {
+      fork <- TRUE
+    } else {
+      fork <- FALSE
+    }
+    if (fork) {
+      output <- parallel::mclapply(
+        X = delta_t,
+        FUN = .DeltaMedStd,
+        phi = phi,
+        sigma = sigma,
+        vcov_theta = vcov_theta,
+        from = from,
+        to = to,
+        med = med,
+        mc.cores = ncores
+      )
+    } else {
+      cl <- parallel::makeCluster(ncores)
+      on.exit(
+        parallel::stopCluster(cl = cl)
+      )
+      output <- parallel::parLapply(
+        cl = cl,
+        X = delta_t,
+        fun = .DeltaMedStd,
+        phi = phi,
+        sigma = sigma,
+        vcov_theta = vcov_theta,
+        from = from,
+        to = to,
+        med = med
+      )
+    }
     # nocov end
   } else {
     output <- lapply(
diff --git a/R/cTMed-delta-med.R b/R/cTMed-delta-med.R
index 31c0bb4..22ea87a 100644
--- a/R/cTMed-delta-med.R
+++ b/R/cTMed-delta-med.R
@@ -267,25 +267,50 @@ DeltaMed <- function(phi,
   par <- FALSE
   if (!is.null(ncores)) {
     ncores <- as.integer(ncores)
+    R <- length(delta_t)
+    if (ncores > R) {
+      ncores <- R
+    }
     if (ncores > 1) {
       par <- TRUE
     }
   }
   if (par) {
-    cl <- parallel::makeCluster(ncores)
-    on.exit(
-      parallel::stopCluster(cl = cl)
-    )
-    output <- parallel::parLapply(
-      cl = cl,
-      X = delta_t,
-      fun = .DeltaMed,
-      phi = phi,
-      vcov_phi_vec = vcov_phi_vec,
-      from = from,
-      to = to,
-      med = med
-    )
+    os_type <- Sys.info()["sysname"]
+    if (os_type == "Darwin") {
+      fork <- TRUE
+    } else if (os_type == "Linux") {
+      fork <- TRUE
+    } else {
+      fork <- FALSE
+    }
+    if (fork) {
+      output <- parallel::mclapply(
+        X = delta_t,
+        FUN = .DeltaMed,
+        phi = phi,
+        vcov_phi_vec = vcov_phi_vec,
+        from = from,
+        to = to,
+        med = med,
+        mc.cores = ncores
+      )
+    } else {
+      cl <- parallel::makeCluster(ncores)
+      on.exit(
+        parallel::stopCluster(cl = cl)
+      )
+      output <- parallel::parLapply(
+        cl = cl,
+        X = delta_t,
+        fun = .DeltaMed,
+        phi = phi,
+        vcov_phi_vec = vcov_phi_vec,
+        from = from,
+        to = to,
+        med = med
+      )
+    }
     # nocov end
   } else {
     output <- lapply(
diff --git a/R/cTMed-delta-total-central.R b/R/cTMed-delta-total-central.R
index bffa137..e104115 100644
--- a/R/cTMed-delta-total-central.R
+++ b/R/cTMed-delta-total-central.R
@@ -226,23 +226,46 @@ DeltaTotalCentral <- function(phi,
   par <- FALSE
   if (!is.null(ncores)) {
     ncores <- as.integer(ncores)
+    R <- length(delta_t)
+    if (ncores > R) {
+      ncores <- R
+    }
     if (ncores > 1) {
       par <- TRUE
     }
   }
   if (par) {
-    cl <- parallel::makeCluster(ncores)
-    on.exit(
-      parallel::stopCluster(cl = cl)
-    )
-    output <- parallel::parLapply(
-      cl = cl,
-      X = delta_t,
-      fun = .DeltaCentral,
-      phi = phi,
-      vcov_phi_vec = vcov_phi_vec,
-      total = total
-    )
+    os_type <- Sys.info()["sysname"]
+    if (os_type == "Darwin") {
+      fork <- TRUE
+    } else if (os_type == "Linux") {
+      fork <- TRUE
+    } else {
+      fork <- FALSE
+    }
+    if (fork) {
+      output <- parallel::mclapply(
+        X = delta_t,
+        FUN = .DeltaCentral,
+        phi = phi,
+        vcov_phi_vec = vcov_phi_vec,
+        total = total,
+        mc.cores = ncores
+      )
+    } else {
+      cl <- parallel::makeCluster(ncores)
+      on.exit(
+        parallel::stopCluster(cl = cl)
+      )
+      output <- parallel::parLapply(
+        cl = cl,
+        X = delta_t,
+        fun = .DeltaCentral,
+        phi = phi,
+        vcov_phi_vec = vcov_phi_vec,
+        total = total
+      )
+    }
     # nocov end
   } else {
     output <- lapply(
diff --git a/R/cTMed-mc-beta-dot.R b/R/cTMed-mc-beta-dot.R
index cbbae74..54ec137 100644
--- a/R/cTMed-mc-beta-dot.R
+++ b/R/cTMed-mc-beta-dot.R
@@ -23,62 +23,120 @@
   par <- FALSE
   if (!is.null(ncores)) {
     ncores <- as.integer(ncores)
+    if (ncores > R) {
+      ncores <- R
+    }
     if (ncores > 1) {
       par <- TRUE
     }
   }
   if (par) {
-    # generate phi
-    cl <- parallel::makeCluster(ncores)
-    on.exit(
-      parallel::stopCluster(cl = cl)
-    )
-    if (!is.null(seed)) {
-      parallel::clusterSetRNGStream(
-        cl = cl,
-        iseed = seed
-      )
+    os_type <- Sys.info()["sysname"]
+    if (os_type == "Darwin") {
+      fork <- TRUE
+    } else if (os_type == "Linux") {
+      fork <- TRUE
+    } else {
+      fork <- FALSE
     }
-    phis <- parallel::parLapply(
-      cl = cl,
-      X = 1:R,
-      fun = function(i) {
-        return(
-          .MCPhiI(
+    if (fork) {
+      # generate phi
+      if (!is.null(seed)) {
+        set.seed(seed)
+      }
+      phis <- parallel::mclapply(
+        X = seq_len(R),
+        FUN = function(i) {
+          return(
+            .MCPhiI(
+              phi = phi,
+              vcov_phi_vec_l = t(chol(vcov_phi_vec)),
+              test_phi = test_phi
+            )
+          )
+        },
+        mc.cores = ncores
+      )
+      output <- lapply(
+        X = delta_t,
+        FUN = function(i) {
+          thetahatstar <- parallel::mclapply(
+            X = phis,
+            FUN = .TotalDeltaT,
+            delta_t = i,
+            mc.cores = ncores
+          )
+          thetahatstar <- do.call(
+            what = "rbind",
+            args = thetahatstar
+          )
+          colnames(thetahatstar) <- varnames
+          est <- .TotalDeltaT(
             phi = phi,
-            vcov_phi_vec_l = t(chol(vcov_phi_vec)),
-            test_phi = test_phi
+            delta_t = i
           )
-        )
-      }
-    )
-    output <- lapply(
-      X = delta_t,
-      FUN = function(i) {
-        thetahatstar <- parallel::parLapply(
+          names(est) <- varnames
+          out <- list(
+            delta_t = i,
+            est = est,
+            thetahatstar = thetahatstar
+          )
+          return(out)
+        }
+      )
+    } else {
+      # generate phi
+      cl <- parallel::makeCluster(ncores)
+      on.exit(
+        parallel::stopCluster(cl = cl)
+      )
+      if (!is.null(seed)) {
+        parallel::clusterSetRNGStream(
           cl = cl,
-          X = phis,
-          fun = .TotalDeltaT,
-          delta_t = i
+          iseed = seed
         )
-        thetahatstar <- do.call(
-          what = "rbind",
-          args = thetahatstar
-        )
-        colnames(thetahatstar) <- varnames
-        est <- .TotalDeltaT(
-          phi = phi,
-          delta_t = i
-        )
-        names(est) <- varnames
-        out <- list(
-          delta_t = i,
-          est = est,
-          thetahatstar = thetahatstar
-        )
-        return(out)
       }
-    )
+      phis <- parallel::parLapply(
+        cl = cl,
+        X = seq_len(R),
+        fun = function(i) {
+          return(
+            .MCPhiI(
+              phi = phi,
+              vcov_phi_vec_l = t(chol(vcov_phi_vec)),
+              test_phi = test_phi
+            )
+          )
+        }
+      )
+      output <- lapply(
+        X = delta_t,
+        FUN = function(i) {
+          thetahatstar <- parallel::parLapply(
+            cl = cl,
+            X = phis,
+            fun = .TotalDeltaT,
+            delta_t = i
+          )
+          thetahatstar <- do.call(
+            what = "rbind",
+            args = thetahatstar
+          )
+          colnames(thetahatstar) <- varnames
+          est <- .TotalDeltaT(
+            phi = phi,
+            delta_t = i
+          )
+          names(est) <- varnames
+          out <- list(
+            delta_t = i,
+            est = est,
+            thetahatstar = thetahatstar
+          )
+          return(out)
+        }
+      )
+    }
     # nocov end
   } else {
     # generate phi
@@ -86,7 +144,7 @@
       set.seed(seed)
     }
     phis <- lapply(
-      X = 1:R,
+      X = seq_len(R),
       FUN = function(i) {
         return(
           .MCPhiI(
diff --git a/R/cTMed-mc-beta-std-dot.R b/R/cTMed-mc-beta-std-dot.R
index a0156b1..ab1c22d 100644
--- a/R/cTMed-mc-beta-std-dot.R
+++ b/R/cTMed-mc-beta-std-dot.R
@@ -28,72 +28,140 @@
   par <- FALSE
   if (!is.null(ncores)) {
     ncores <- as.integer(ncores)
+    if (ncores > R) {
+      ncores <- R
+    }
     if (ncores > 1) {
       par <- TRUE
     }
   }
   if (par) {
-    # generate phi
-    cl <- parallel::makeCluster(ncores)
-    on.exit(
-      parallel::stopCluster(cl = cl)
-    )
-    if (!is.null(seed)) {
-      parallel::clusterSetRNGStream(
-        cl = cl,
-        iseed = seed
-      )
+    os_type <- Sys.info()["sysname"]
+    if (os_type == "Darwin") {
+      fork <- TRUE
+    } else if (os_type == "Linux") {
+      fork <- TRUE
+    } else {
+      fork <- FALSE
     }
-    phisigmas <- parallel::parLapply(
-      cl = cl,
-      X = 1:R,
-      fun = function(i) {
-        return(
-          .MCPhiSigmaI(
-            theta = theta,
-            vcov_theta = vcov_theta,
-            test_phi = test_phi
-          )
-        )
+    if (fork) {
+      # generate phi
+      if (!is.null(seed)) {
+        set.seed(seed)
       }
-    )
-    output <- lapply(
-      X = delta_t,
-      FUN = function(i) {
-        thetahatstar <- parallel::parLapply(
-          cl = cl,
-          X = phisigmas,
-          fun = function(x,
-                         delta_t) {
-            return(
-              .TotalStdDeltaT(
-                phi = x[[1]],
-                sigma = x[[2]],
-                delta_t = delta_t
-              )
+      phisigmas <- parallel::mclapply(
+        X = seq_len(R),
+        FUN = function(i) {
+          return(
+            .MCPhiSigmaI(
+              theta = theta,
+              vcov_theta = vcov_theta,
+              test_phi = test_phi
             )
-          },
-          delta_t = i
-        )
-        thetahatstar <- do.call(
-          what = "rbind",
-          args = thetahatstar
-        )
-        colnames(thetahatstar) <- varnames
-        est <- .TotalStdDeltaT(
-          phi = phi,
-          sigma = sigma,
-          delta_t = i
-        )
-        names(est) <- varnames
-        out <- list(
-          delta_t = i,
-          est = est,
-          thetahatstar = thetahatstar
+          )
+        },
+        mc.cores = ncores
+      )
+      output <- lapply(
+        X = delta_t,
+        FUN = function(i) {
+          thetahatstar <- parallel::mclapply(
+            X = phisigmas,
+            FUN = function(x,
+                           delta_t) {
+              return(
+                .TotalStdDeltaT(
+                  phi = x[[1]],
+                  sigma = x[[2]],
+                  delta_t = delta_t
+                )
+              )
+            },
+            delta_t = i,
+            mc.cores = ncores
+          )
+          thetahatstar <- do.call(
+            what = "rbind",
+            args = thetahatstar
+          )
+          colnames(thetahatstar) <- varnames
+          est <- .TotalStdDeltaT(
+            phi = phi,
+            sigma = sigma,
+            delta_t = i
+          )
+          names(est) <- varnames
+          out <- list(
+            delta_t = i,
+            est = est,
+            thetahatstar = thetahatstar
+          )
+          return(out)
+        }
+      )
+    } else {
+      # generate phi
+      cl <- parallel::makeCluster(ncores)
+      on.exit(
+        parallel::stopCluster(cl = cl)
+      )
+      if (!is.null(seed)) {
+        parallel::clusterSetRNGStream(
+          cl = cl,
+          iseed = seed
         )
-        return(out)
       }
-    )
+      phisigmas <- parallel::parLapply(
+        cl = cl,
+        X = seq_len(R),
+        fun = function(i) {
+          return(
+            .MCPhiSigmaI(
+              theta = theta,
+              vcov_theta = vcov_theta,
+              test_phi = test_phi
+            )
+          )
+        }
+      )
+      output <- lapply(
+        X = delta_t,
+        FUN = function(i) {
+          thetahatstar <- parallel::parLapply(
+            cl = cl,
+            X = phisigmas,
+            fun = function(x,
+                           delta_t) {
+              return(
+                .TotalStdDeltaT(
+                  phi = x[[1]],
+                  sigma = x[[2]],
+                  delta_t = delta_t
+                )
+              )
+            },
+            delta_t = i
+          )
+          thetahatstar <- do.call(
+            what = "rbind",
+            args = thetahatstar
+          )
+          colnames(thetahatstar) <- varnames
+          est <- .TotalStdDeltaT(
+            phi = phi,
+            sigma = sigma,
+            delta_t = i
+          )
+          names(est) <- varnames
+          out <- list(
+            delta_t = i,
+            est = est,
+            thetahatstar = thetahatstar
+          )
+          return(out)
+        }
+      )
+    }
     # nocov end
   } else {
     # generate phi
@@ -101,7 +169,7 @@
       set.seed(seed)
     }
     phisigmas <- lapply(
-      X = 1:R,
+      X = seq_len(R),
       FUN = function(i) {
         return(
           .MCPhiSigmaI(
diff --git a/R/cTMed-mc-central-dot.R b/R/cTMed-mc-central-dot.R
index c92e480..1de6c06 100644
--- a/R/cTMed-mc-central-dot.R
+++ b/R/cTMed-mc-central-dot.R
@@ -15,72 +15,140 @@
   par <- FALSE
   if (!is.null(ncores)) {
     ncores <- as.integer(ncores)
+    if (ncores > R) {
+      ncores <- R
+    }
     if (ncores > 1) {
       par <- TRUE
     }
   }
   if (par) {
-    # generate phi
-    cl <- parallel::makeCluster(ncores)
-    on.exit(
-      parallel::stopCluster(cl = cl)
-    )
-    if (!is.null(seed)) {
-      parallel::clusterSetRNGStream(
-        cl = cl,
-        iseed = seed
-      )
+    os_type <- Sys.info()["sysname"]
+    if (os_type == "Darwin") {
+      fork <- TRUE
+    } else if (os_type == "Linux") {
+      fork <- TRUE
+    } else {
+      fork <- FALSE
     }
-    phis <- parallel::parLapply(
-      cl = cl,
-      X = 1:R,
-      fun = function(i) {
-        return(
-          .MCPhiI(
-            phi = phi,
-            vcov_phi_vec_l = t(chol(vcov_phi_vec)),
-            test_phi = test_phi
-          )
-        )
+    if (fork) {
+      # generate phi
+      if (!is.null(seed)) {
+        set.seed(seed)
       }
-    )
-    output <- lapply(
-      X = delta_t,
-      FUN = function(i) {
-        thetahatstar <- parallel::parLapply(
+      phis <- parallel::mclapply(
+        X = seq_len(R),
+        FUN = function(i) {
+          return(
+            .MCPhiI(
+              phi = phi,
+              vcov_phi_vec_l = t(chol(vcov_phi_vec)),
+              test_phi = test_phi
+            )
+          )
+        },
+        mc.cores = ncores
+      )
+      output <- lapply(
+        X = delta_t,
+        FUN = function(i) {
+          thetahatstar <- parallel::mclapply(
+            X = phis,
+            FUN = Fun,
+            delta_t = i,
+            mc.cores = ncores
+          )
+          thetahatstar <- do.call(
+            what = "rbind",
+            args = thetahatstar
+          )
+          colnames(thetahatstar) <- colnames(phi)
+          thetahatstar <- cbind(
+            thetahatstar,
+            interval = i
+          )
+          est <- c(
+            Fun(
+              phi = phi,
+              delta_t = i
+            ),
+            i
+          )
+          names(est) <- c(
+            colnames(phi),
+            "interval"
+          )
+          out <- list(
+            delta_t = i,
+            est = est,
+            thetahatstar = thetahatstar
+          )
+          return(out)
+        }
+      )
+    } else {
+      # generate phi
+      cl <- parallel::makeCluster(ncores)
+      on.exit(
+        parallel::stopCluster(cl = cl)
+      )
+      if (!is.null(seed)) {
+        parallel::clusterSetRNGStream(
           cl = cl,
-          X = phis,
-          fun = Fun,
-          delta_t = i
-        )
-        thetahatstar <- do.call(
-          what = "rbind",
-          args = thetahatstar
-        )
-        colnames(thetahatstar) <- colnames(phi)
-        thetahatstar <- cbind(
-          thetahatstar,
-          interval = i
-        )
-        est <- c(
-          Fun(
-            phi = phi,
-            delta_t = i
-          ),
-          i
-        )
-        names(est) <- c(
-          colnames(phi),
-          "interval"
-        )
-        out <- list(
-          delta_t = i,
-          est = est,
-          thetahatstar = thetahatstar
+          iseed = seed
         )
-        return(out)
       }
-    )
+      phis <- parallel::parLapply(
+        cl = cl,
+        X = seq_len(R),
+        fun = function(i) {
+          return(
+            .MCPhiI(
+              phi = phi,
+              vcov_phi_vec_l = t(chol(vcov_phi_vec)),
+              test_phi = test_phi
+            )
+          )
+        }
+      )
+      output <- lapply(
+        X = delta_t,
+        FUN = function(i) {
+          thetahatstar <- parallel::parLapply(
+            cl = cl,
+            X = phis,
+            fun = Fun,
+            delta_t = i
+          )
+          thetahatstar <- do.call(
+            what = "rbind",
+            args = thetahatstar
+          )
+          colnames(thetahatstar) <- colnames(phi)
+          thetahatstar <- cbind(
+            thetahatstar,
+            interval = i
+          )
+          est <- c(
+            Fun(
+              phi = phi,
+              delta_t = i
+            ),
+            i
+          )
+          names(est) <- c(
+            colnames(phi),
+            "interval"
+          )
+          out <- list(
+            delta_t = i,
+            est = est,
+            thetahatstar = thetahatstar
+          )
+          return(out)
+        }
+      )
+    }
     # nocov end
   } else {
     # generate phi
@@ -88,7 +156,7 @@
       set.seed(seed)
     }
     phis <- lapply(
-      X = 1:R,
+      X = seq_len(R),
       FUN = function(i) {
         return(
           .MCPhiI(
diff --git a/R/cTMed-mc-med-dot.R b/R/cTMed-mc-med-dot.R
index 0669886..3666e16 100644
--- a/R/cTMed-mc-med-dot.R
+++ b/R/cTMed-mc-med-dot.R
@@ -12,78 +12,152 @@
   par <- FALSE
   if (!is.null(ncores)) {
     ncores <- as.integer(ncores)
+    if (ncores > R) {
+      ncores <- R
+    }
     if (ncores > 1) {
       par <- TRUE
     }
   }
   if (par) {
-    # generate phi
-    cl <- parallel::makeCluster(ncores)
-    on.exit(
-      parallel::stopCluster(cl = cl)
-    )
-    if (!is.null(seed)) {
-      parallel::clusterSetRNGStream(
-        cl = cl,
-        iseed = seed
-      )
+    os_type <- Sys.info()["sysname"]
+    if (os_type == "Darwin") {
+      fork <- TRUE
+    } else if (os_type == "Linux") {
+      fork <- TRUE
+    } else {
+      fork <- FALSE
     }
-    phis <- parallel::parLapply(
-      cl = cl,
-      X = 1:R,
-      fun = function(i) {
-        return(
-          .MCPhiI(
+    if (fork) {
+      # generate phi
+      if (!is.null(seed)) {
+        set.seed(seed)
+      }
+      phis <- parallel::mclapply(
+        X = seq_len(R),
+        FUN = function(i) {
+          return(
+            .MCPhiI(
+              phi = phi,
+              vcov_phi_vec_l = t(chol(vcov_phi_vec)),
+              test_phi = test_phi
+            )
+          )
+        },
+        mc.cores = ncores
+      )
+      output <- lapply(
+        X = delta_t,
+        FUN = function(i) {
+          thetahatstar <- parallel::mclapply(
+            X = phis,
+            FUN = .Med,
+            delta_t = i,
+            from = from,
+            to = to,
+            med = med,
+            mc.cores = ncores
+          )
+          thetahatstar <- do.call(
+            what = "rbind",
+            args = thetahatstar
+          )
+          colnames(thetahatstar) <- c(
+            "total",
+            "direct",
+            "indirect",
+            "interval"
+          )
+          est <- .Med(
             phi = phi,
-            vcov_phi_vec_l = t(chol(vcov_phi_vec)),
-            test_phi = test_phi
+            delta_t = i,
+            from = from,
+            to = to,
+            med = med
           )
-        )
-      }
-    )
-    output <- lapply(
-      X = delta_t,
-      FUN = function(i) {
-        thetahatstar <- parallel::parLapply(
+          names(est) <- c(
+            "total",
+            "direct",
+            "indirect",
+            "interval"
+          )
+          out <- list(
+            delta_t = i,
+            est = est,
+            thetahatstar = thetahatstar
+          )
+          return(out)
+        }
+      )
+    } else {
+      # generate phi
+      cl <- parallel::makeCluster(ncores)
+      on.exit(
+        parallel::stopCluster(cl = cl)
+      )
+      if (!is.null(seed)) {
+        parallel::clusterSetRNGStream(
           cl = cl,
-          X = phis,
-          fun = .Med,
-          delta_t = i,
-          from = from,
-          to = to,
-          med = med
+          iseed = seed
         )
-        thetahatstar <- do.call(
-          what = "rbind",
-          args = thetahatstar
-        )
-        colnames(thetahatstar) <- c(
-          "total",
-          "direct",
-          "indirect",
-          "interval"
-        )
-        est <- .Med(
-          phi = phi,
-          delta_t = i,
-          from = from,
-          to = to,
-          med = med
-        )
-        names(est) <- c(
-          "total",
-          "direct",
-          "indirect",
-          "interval"
-        )
-        out <- list(
-          delta_t = i,
-          est = est,
-          thetahatstar = thetahatstar
-        )
-        return(out)
       }
-    )
+      phis <- parallel::parLapply(
+        cl = cl,
+        X = seq_len(R),
+        fun = function(i) {
+          return(
+            .MCPhiI(
+              phi = phi,
+              vcov_phi_vec_l = t(chol(vcov_phi_vec)),
+              test_phi = test_phi
+            )
+          )
+        }
+      )
+      output <- lapply(
+        X = delta_t,
+        FUN = function(i) {
+          thetahatstar <- parallel::parLapply(
+            cl = cl,
+            X = phis,
+            fun = .Med,
+            delta_t = i,
+            from = from,
+            to = to,
+            med = med
+          )
+          thetahatstar <- do.call(
+            what = "rbind",
+            args = thetahatstar
+          )
+          colnames(thetahatstar) <- c(
+            "total",
+            "direct",
+            "indirect",
+            "interval"
+          )
+          est <- .Med(
+            phi = phi,
+            delta_t = i,
+            from = from,
+            to = to,
+            med = med
+          )
+          names(est) <- c(
+            "total",
+            "direct",
+            "indirect",
+            "interval"
+          )
+          out <- list(
+            delta_t = i,
+            est = est,
+            thetahatstar = thetahatstar
+          )
+          return(out)
+        }
+      )
+    }
     # nocov end
   } else {
     # generate phi
@@ -91,7 +165,7 @@
       set.seed(seed)
     }
     phis <- lapply(
-      X = 1:R,
+      X = seq_len(R),
       FUN = function(i) {
         return(
           .MCPhiI(
diff --git a/R/cTMed-mc-med-std-dot.R b/R/cTMed-mc-med-std-dot.R
index 5defab8..dec5438 100644
--- a/R/cTMed-mc-med-std-dot.R
+++ b/R/cTMed-mc-med-std-dot.R
@@ -17,93 +17,182 @@
   par <- FALSE
   if (!is.null(ncores)) {
     ncores <- as.integer(ncores)
+    if (ncores > R) {
+      ncores <- R
+    }
     if (ncores > 1) {
       par <- TRUE
     }
   }
   if (par) {
-    # generate phi
-    cl <- parallel::makeCluster(ncores)
-    on.exit(
-      parallel::stopCluster(cl = cl)
-    )
-    if (!is.null(seed)) {
-      parallel::clusterSetRNGStream(
-        cl = cl,
-        iseed = seed
-      )
+    os_type <- Sys.info()["sysname"]
+    if (os_type == "Darwin") {
+      fork <- TRUE
+    } else if (os_type == "Linux") {
+      fork <- TRUE
+    } else {
+      fork <- FALSE
     }
-    phis <- parallel::parLapply(
-      cl = cl,
-      X = 1:R,
-      fun = function(i) {
-        return(
-          .MCPhiSigmaI(
-            theta = theta,
-            vcov_theta = vcov_theta,
-            test_phi = test_phi
-          )
-        )
+    if (fork) {
+      # generate phi
+      if (!is.null(seed)) {
+        set.seed(seed)
       }
-    )
-    output <- lapply(
-      X = delta_t,
-      FUN = function(i) {
-        thetahatstar <- parallel::parLapply(
-          cl = cl,
-          X = phis,
-          fun = function(x,
-                         delta_t,
-                         from,
-                         to,
-                         med) {
-            return(
-              .MedStd(
-                phi = x[[1]],
-                sigma = x[[2]],
-                delta_t = delta_t,
-                from = from,
-                to = to,
-                med = med
-              )
+      phis <- parallel::mclapply(
+        X = seq_len(R),
+        FUN = function(i) {
+          return(
+            .MCPhiSigmaI(
+              theta = theta,
+              vcov_theta = vcov_theta,
+              test_phi = test_phi
             )
-          },
-          delta_t = i,
-          from = from,
-          to = to,
-          med = med
-        )
-        thetahatstar <- do.call(
-          what = "rbind",
-          args = thetahatstar
-        )
-        colnames(thetahatstar) <- c(
-          "total",
-          "direct",
-          "indirect",
-          "interval"
-        )
-        est <- .Med(
-          phi = phi,
-          delta_t = i,
-          from = from,
-          to = to,
-          med = med
-        )
-        names(est) <- c(
-          "total",
-          "direct",
-          "indirect",
-          "interval"
-        )
-        out <- list(
-          delta_t = i,
-          est = est,
-          thetahatstar = thetahatstar
+          )
+        },
+        mc.cores = ncores
+      )
+      output <- lapply(
+        X = delta_t,
+        FUN = function(i) {
+          thetahatstar <- parallel::mclapply(
+            X = phis,
+            FUN = function(x,
+                           delta_t,
+                           from,
+                           to,
+                           med) {
+              return(
+                .MedStd(
+                  phi = x[[1]],
+                  sigma = x[[2]],
+                  delta_t = delta_t,
+                  from = from,
+                  to = to,
+                  med = med
+                )
+              )
+            },
+            delta_t = i,
+            from = from,
+            to = to,
+            med = med,
+            mc.cores = ncores
+          )
+          thetahatstar <- do.call(
+            what = "rbind",
+            args = thetahatstar
+          )
+          colnames(thetahatstar) <- c(
+            "total",
+            "direct",
+            "indirect",
+            "interval"
+          )
+          est <- .Med(
+            phi = phi,
+            delta_t = i,
+            from = from,
+            to = to,
+            med = med
+          )
+          names(est) <- c(
+            "total",
+            "direct",
+            "indirect",
+            "interval"
+          )
+          out <- list(
+            delta_t = i,
+            est = est,
+            thetahatstar = thetahatstar
+          )
+          return(out)
+        }
+      )
+    } else {
+      # generate phi
+      cl <- parallel::makeCluster(ncores)
+      on.exit(
+        parallel::stopCluster(cl = cl)
+      )
+      if (!is.null(seed)) {
+        parallel::clusterSetRNGStream(
+          cl = cl,
+          iseed = seed
         )
-        return(out)
       }
-    )
+      phis <- parallel::parLapply(
+        cl = cl,
+        X = seq_len(R),
+        fun = function(i) {
+          return(
+            .MCPhiSigmaI(
+              theta = theta,
+              vcov_theta = vcov_theta,
+              test_phi = test_phi
+            )
+          )
+        }
+      )
+      output <- lapply(
+        X = delta_t,
+        FUN = function(i) {
+          thetahatstar <- parallel::parLapply(
+            cl = cl,
+            X = phis,
+            fun = function(x,
+                           delta_t,
+                           from,
+                           to,
+                           med) {
+              return(
+                .MedStd(
+                  phi = x[[1]],
+                  sigma = x[[2]],
+                  delta_t = delta_t,
+                  from = from,
+                  to = to,
+                  med = med
+                )
+              )
+            },
+            delta_t = i,
+            from = from,
+            to = to,
+            med = med
+          )
+          thetahatstar <- do.call(
+            what = "rbind",
+            args = thetahatstar
+          )
+          colnames(thetahatstar) <- c(
+            "total",
+            "direct",
+            "indirect",
+            "interval"
+          )
+          est <- .Med(
+            phi = phi,
+            delta_t = i,
+            from = from,
+            to = to,
+            med = med
+          )
+          names(est) <- c(
+            "total",
+            "direct",
+            "indirect",
+            "interval"
+          )
+          out <- list(
+            delta_t = i,
+            est = est,
+            thetahatstar = thetahatstar
+          )
+          return(out)
+        }
+      )
+    }
     # nocov end
   } else {
     # generate phi
@@ -111,7 +200,7 @@
       set.seed(seed)
     }
     phis <- lapply(
-      X = 1:R,
+      X = seq_len(R),
       FUN = function(i) {
         return(
           .MCPhiSigmaI(
diff --git a/R/cTMed-mc-phi.R b/R/cTMed-mc-phi.R
index 15ab131..8eef8fb 100644
--- a/R/cTMed-mc-phi.R
+++ b/R/cTMed-mc-phi.R
@@ -135,34 +135,64 @@ MCPhi <- function(phi,
   par <- FALSE
   if (!is.null(ncores)) {
     ncores <- as.integer(ncores)
+    if (ncores > R) {
+      ncores <- R
+    }
     if (ncores > 1) {
       par <- TRUE
     }
   }
   if (par) {
-    cl <- parallel::makeCluster(ncores)
-    on.exit(
-      parallel::stopCluster(cl = cl)
-    )
-    if (!is.null(seed)) {
-      parallel::clusterSetRNGStream(
-        cl = cl,
-        iseed = seed
-      )
+    os_type <- Sys.info()["sysname"]
+    if (os_type == "Darwin") {
+      fork <- TRUE
+    } else if (os_type == "Linux") {
+      fork <- TRUE
+    } else {
+      fork <- FALSE
     }
-    output <- parallel::parLapply(
-      cl = cl,
-      X = 1:R,
-      fun = function(i) {
-        return(
-          .MCPhiI(
-            phi = phi,
-            vcov_phi_vec_l = t(chol(vcov_phi_vec)),
-            test_phi = test_phi
+    if (fork) {
+      if (!is.null(seed)) {
+        set.seed(seed)
+      }
+      output <- parallel::mclapply(
+        X = seq_len(R),
+        FUN = function(i) {
+          return(
+            .MCPhiI(
+              phi = phi,
+              vcov_phi_vec_l = t(chol(vcov_phi_vec)),
+              test_phi = test_phi
+            )
           )
+        },
+        mc.cores = ncores
+      )
+    } else {
+      cl <- parallel::makeCluster(ncores)
+      on.exit(
+        parallel::stopCluster(cl = cl)
+      )
+      if (!is.null(seed)) {
+        parallel::clusterSetRNGStream(
+          cl = cl,
+          iseed = seed
         )
       }
-    )
+      output <- parallel::parLapply(
+        cl = cl,
+        X = seq_len(R),
+        fun = function(i) {
+          return(
+            .MCPhiI(
+              phi = phi,
+              vcov_phi_vec_l = t(chol(vcov_phi_vec)),
+              test_phi = test_phi
+            )
+          )
+        }
+      )
+    }
     # nocov end
   } else {
     if (!is.null(seed)) {
diff --git a/R/cTMed-posterior-beta-dot.R b/R/cTMed-posterior-beta-dot.R
index 9044b46..7dd07c5 100644
--- a/R/cTMed-posterior-beta-dot.R
+++ b/R/cTMed-posterior-beta-dot.R
@@ -36,42 +36,84 @@
   par <- FALSE
   if (!is.null(ncores)) {
     ncores <- as.integer(ncores)
+    R <- length(delta_t)
+    if (ncores > R) {
+      ncores <- R
+    }
     if (ncores > 1) {
       par <- TRUE
     }
   }
   if (par) {
-    cl <- parallel::makeCluster(ncores)
-    on.exit(
-      parallel::stopCluster(cl = cl)
-    )
-    output <- lapply(
-      X = delta_t,
-      FUN = function(i) {
-        thetahatstar <- parallel::parLapply(
-          cl = cl,
-          X = phi,
-          fun = .TotalDeltaT,
-          delta_t = i
-        )
-        thetahatstar <- do.call(
-          what = "rbind",
-          args = thetahatstar
-        )
-        colnames(thetahatstar) <- varnames
-        est <- .TotalDeltaT(
-          phi = phi_mean,
-          delta_t = i
-        )
-        names(est) <- varnames
-        out <- list(
-          delta_t = i,
-          est = est,
-          thetahatstar = thetahatstar
-        )
-        return(out)
-      }
-    )
+    os_type <- Sys.info()["sysname"]
+    if (os_type == "Darwin") {
+      fork <- TRUE
+    } else if (os_type == "Linux") {
+      fork <- TRUE
+    } else {
+      fork <- FALSE
+    }
+    if (fork) {
+      output <- lapply(
+        X = delta_t,
+        FUN = function(i) {
+          thetahatstar <- parallel::mclapply(
+            X = phi,
+            FUN = .TotalDeltaT,
+            delta_t = i,
+            mc.cores = ncores
+          )
+          thetahatstar <- do.call(
+            what = "rbind",
+            args = thetahatstar
+          )
+          colnames(thetahatstar) <- varnames
+          est <- .TotalDeltaT(
+            phi = phi_mean,
+            delta_t = i
+          )
+          names(est) <- varnames
+          out <- list(
+            delta_t = i,
+            est = est,
+            thetahatstar = thetahatstar
+          )
+          return(out)
+        }
+      )
+    } else {
+      cl <- parallel::makeCluster(ncores)
+      on.exit(
+        parallel::stopCluster(cl = cl)
+      )
+      output <- lapply(
+        X = delta_t,
+        FUN = function(i) {
+          thetahatstar <- parallel::parLapply(
+            cl = cl,
+            X = phi,
+            fun = .TotalDeltaT,
+            delta_t = i
+          )
+          thetahatstar <- do.call(
+            what = "rbind",
+            args = thetahatstar
+          )
+          colnames(thetahatstar) <- varnames
+          est <- .TotalDeltaT(
+            phi = phi_mean,
+            delta_t = i
+          )
+          names(est) <- varnames
+          out <- list(
+            delta_t = i,
+            est = est,
+            thetahatstar = thetahatstar
+          )
+          return(out)
+        }
+      )
+    }
     # nocov end
   } else {
     output <- lapply(
diff --git a/R/cTMed-posterior-central-dot.R b/R/cTMed-posterior-central-dot.R
index f24a492..7fa9743 100644
--- a/R/cTMed-posterior-central-dot.R
+++ b/R/cTMed-posterior-central-dot.R
@@ -11,51 +11,102 @@
   par <- FALSE
   if (!is.null(ncores)) {
     ncores <- as.integer(ncores)
+    R <- length(delta_t)
+    if (ncores > R) {
+      ncores <- R
+    }
     if (ncores > 1) {
       par <- TRUE
     }
   }
   if (par) {
-    cl <- parallel::makeCluster(ncores)
-    on.exit(
-      parallel::stopCluster(cl = cl)
-    )
-    output <- lapply(
-      X = delta_t,
-      FUN = function(i) {
-        thetahatstar <- parallel::parLapply(
-          cl = cl,
-          X = phi,
-          fun = function(phi) {
-            out <- c(
-              Fun(
-                phi = phi,
-                delta_t = i
-              ),
-              i
-            )
-            names(out) <- c(
-              colnames(phi),
-              "interval"
+    os_type <- Sys.info()["sysname"]
+    if (os_type == "Darwin") {
+      fork <- TRUE
+    } else if (os_type == "Linux") {
+      fork <- TRUE
+    } else {
+      fork <- FALSE
+    }
+    if (fork) {
+      output <- lapply(
+        X = delta_t,
+        FUN = function(i) {
+          thetahatstar <- parallel::mclapply(
+            X = phi,
+            FUN = function(phi) {
+              out <- c(
+                Fun(
+                  phi = phi,
+                  delta_t = i
+                ),
+                i
+              )
+              names(out) <- c(
+                colnames(phi),
+                "interval"
+              )
+              return(
+                out
+              )
+            },
+            mc.cores = ncores
+          )
+          thetahatstar <- do.call(
+            what = "rbind",
+            args = thetahatstar
+          )
+          return(
+            list(
+              delta_t = i,
+              est = colMeans(thetahatstar),
+              thetahatstar = thetahatstar
             )
-            return(
-              out
+          )
+        }
+      )
+    } else {
+      cl <- parallel::makeCluster(ncores)
+      on.exit(
+        parallel::stopCluster(cl = cl)
+      )
+      output <- lapply(
+        X = delta_t,
+        FUN = function(i) {
+          thetahatstar <- parallel::parLapply(
+            cl = cl,
+            X = phi,
+            fun = function(phi) {
+              out <- c(
+                Fun(
+                  phi = phi,
+                  delta_t = i
+                ),
+                i
+              )
+              names(out) <- c(
+                colnames(phi),
+                "interval"
+              )
+              return(
+                out
+              )
+            }
+          )
+          thetahatstar <- do.call(
+            what = "rbind",
+            args = thetahatstar
+          )
+          return(
+            list(
+              delta_t = i,
+              est = colMeans(thetahatstar),
+              thetahatstar = thetahatstar
             )
-          }
-        )
-        thetahatstar <- do.call(
-          what = "rbind",
-          args = thetahatstar
-        )
-        return(
-          list(
-            delta_t = i,
-            est = colMeans(thetahatstar),
-            thetahatstar = thetahatstar
           )
-        )
-      }
-    )
+        }
+      )
+    }
     # nocov end
   } else {
     output <- lapply(
diff --git a/R/cTMed-posterior-med-dot.R b/R/cTMed-posterior-med-dot.R
index ed1f469..fdc493d 100644
--- a/R/cTMed-posterior-med-dot.R
+++ b/R/cTMed-posterior-med-dot.R
@@ -8,53 +8,106 @@
   par <- FALSE
   if (!is.null(ncores)) {
     ncores <- as.integer(ncores)
+    R <- length(delta_t)
+    if (ncores > R) {
+      ncores <- R
+    }
     if (ncores > 1) {
       par <- TRUE
     }
   }
   if (par) {
-    cl <- parallel::makeCluster(ncores)
-    on.exit(
-      parallel::stopCluster(cl = cl)
-    )
-    output <- lapply(
-      X = delta_t,
-      FUN = function(i) {
-        thetahatstar <- parallel::parLapply(
-          cl = cl,
-          X = phi,
-          fun = function(phi) {
-            out <- .Med(
-              phi = phi,
+    os_type <- Sys.info()["sysname"]
+    if (os_type == "Darwin") {
+      fork <- TRUE
+    } else if (os_type == "Linux") {
+      fork <- TRUE
+    } else {
+      fork <- FALSE
+    }
+    if (fork) {
+      output <- lapply(
+        X = delta_t,
+        FUN = function(i) {
+          thetahatstar <- parallel::mclapply(
+            X = phi,
+            FUN = function(phi) {
+              out <- .Med(
+                phi = phi,
+                delta_t = i,
+                from = from,
+                to = to,
+                med = med
+              )
+              names(out) <- c(
+                "total",
+                "direct",
+                "indirect",
+                "interval"
+              )
+              return(
+                out
+              )
+            },
+            mc.cores = ncores
+          )
+          thetahatstar <- do.call(
+            what = "rbind",
+            args = thetahatstar
+          )
+          return(
+            list(
               delta_t = i,
-              from = from,
-              to = to,
-              med = med
-            )
-            names(out) <- c(
-              "total",
-              "direct",
-              "indirect",
-              "interval"
+              est = colMeans(thetahatstar),
+              thetahatstar = thetahatstar
             )
-            return(
-              out
+          )
+        }
+      )
+    } else {
+      cl <- parallel::makeCluster(ncores)
+      on.exit(
+        parallel::stopCluster(cl = cl)
+      )
+      output <- lapply(
+        X = delta_t,
+        FUN = function(i) {
+          thetahatstar <- parallel::parLapply(
+            cl = cl,
+            X = phi,
+            fun = function(phi) {
+              out <- .Med(
+                phi = phi,
+                delta_t = i,
+                from = from,
+                to = to,
+                med = med
+              )
+              names(out) <- c(
+                "total",
+                "direct",
+                "indirect",
+                "interval"
+              )
+              return(
+                out
+              )
+            }
+          )
+          thetahatstar <- do.call(
+            what = "rbind",
+            args = thetahatstar
+          )
+          return(
+            list(
+              delta_t = i,
+              est = colMeans(thetahatstar),
+              thetahatstar = thetahatstar
             )
-          }
-        )
-        thetahatstar <- do.call(
-          what = "rbind",
-          args = thetahatstar
-        )
-        return(
-          list(
-            delta_t = i,
-            est = colMeans(thetahatstar),
-            thetahatstar = thetahatstar
           )
-        )
-      }
-    )
+        }
+      )
+    }
     # nocov end
   } else {
     output <- lapply(