diff --git a/tests/testthat/test-inference-help.R b/tests/testthat/test-inference-help.R
index 30e3c57..efd4e8b 100644
--- a/tests/testthat/test-inference-help.R
+++ b/tests/testthat/test-inference-help.R
@@ -1,48 +1,118 @@
-# make sure that you are using Mac Os or Linux (mclapply is currently not supported for Windows unless some mclapply hack function for Windows is preloaded in your R session)
-
-#*********************** IMPORTANT******************************************************
-# if a multithreaded backend openBLAS for matrix multiplications
-# is installed on your machine, please force it to use 1 thread explicitly
-# unless ncores in LogrRegr is reasonably small, in the latter case
-# you might want to experiment with the combinations of blas_set_num_threads and ncores
-library(RhpcBLASctl)
-blas_set_num_threads(1)
-omp_set_num_threads(1)
-#*********************** IMPORTANT******************************************************
-
 # simulate Gaussian responses
-
+threads <- 1L
 set.seed(040590)
-X1 <- as.data.frame(array(data = rbinom(n = 50 * 1000, size = 1, prob = runif(n = 50 * 1000, 0, 1)), dim = c(1000, 50)))
-Y1 <- rnorm(n = 1000, mean = 1 + 0.7 * (X1$V1 * X1$V4) + 0.8896846 * (X1$V8 * X1$V11) + 1.434573 * (X1$V5 * X1$V9), sd = 1)
+X1 <- as.data.frame(
+  array(
+    data = rbinom(n = 50 * 1000, size = 1, prob = runif(n = 50 * 1000, 0, 1)),
+    dim = c(1000, 50)
+  )
+)
+Y1 <- rnorm(
+  n = 1000,
+  mean = 1 + 0.7 * (X1$V1 * X1$V4) + 0.8896846 * (X1$V8 * X1$V11) + 1.434573 * (X1$V5 * X1$V9),
+  sd = 1
+)
 X1$Y1 <- Y1
 
 # specify the initial formula
-formula1 <- as.formula(paste(colnames(X1)[51], "~ 1 +", paste0(colnames(X1)[-c(51)], collapse = "+")))
+formula1 <- as.formula(
+  paste(colnames(X1)[51], "~ 1 +", paste0(colnames(X1)[-c(51)], collapse = "+"))
+)
 data.example <- as.data.frame(X1)
 
 # run the inference with robust g prior
-res4G <- LogicRegr(formula = formula1, data = data.example, family = "Gaussian", prior = "G", report.level = 0.5, d = 15, cmax = 2, kmax = 15, p.and = 0.9, p.not = 0.01, p.surv = 0.2, ncores = 32)
-print(res4G$feat.stat)
+res4G <- suppressMessages(
+  LogicRegr(
+    formula = formula1, data = data.example, family = "Gaussian", prior = "G",
+    report.level = 0.5, d = 15, cmax = 2, kmax = 15, p.and = 0.9, p.not = 0.01,
+    p.surv = 0.2, ncores = threads, print.freq = 0L
+  )
+)
 # run the inference with Jeffrey's prior
-res4J <- LogicRegr(formula = formula1, data = data.example, family = "Gaussian", prior = "J", report.level = 0.5, d = 15, cmax = 2, kmax = 15, p.and = 0.9, p.not = 0.01, p.surv = 0.2, ncores = 32)
-print(res4J$feat.stat)
-
+res4J <- suppressMessages(
+  LogicRegr(
+    formula = formula1, data = data.example, family = "Gaussian", prior = "J",
+    report.level = 0.5, d = 15, cmax = 2, kmax = 15, p.and = 0.9, p.not = 0.01,
+    p.surv = 0.2, ncores = threads, print.freq = 0L
+  )
+)
 
 # change to Bernoulli responses
-X1 <- as.data.frame(array(data = rbinom(n = 50 * 1000, size = 1, prob = 0.3), dim = c(1000, 50)))
+X1 <- as.data.frame(
+  array(data = rbinom(n = 50 * 1000, size = 1, prob = 0.3), dim = c(1000, 50))
+)
 Y1 <- -0.7 + 1 * ((1 - X1$V1) * (X1$V4)) + 1 * (X1$V8 * X1$V11) + 1 * (X1$V5 * X1$V9)
 X1$Y1 <- round(1.0 / (1.0 + exp(-Y1)))
 
 # specify the initial formula
-formula1 <- as.formula(paste(colnames(X1)[51], "~ 1 +", paste0(colnames(X1)[-c(51)], collapse = "+")))
+formula1 <- as.formula(
+  paste(colnames(X1)[51], "~ 1 +", paste0(colnames(X1)[-c(51)], collapse = "+"))
+)
 data.example <- as.data.frame(X1)
 
-
 # run the inference with robust g prior
-res1G <- LogicRegr(formula = formula1, data = data.example, family = "Bernoulli", prior = "G", report.level = 0.5, d = 15, cmax = 2, kmax = 15, p.and = 0.9, p.not = 0.2, p.surv = 0.2, ncores = 32)
-print(res1G$feat.stat)
+res1G <- suppressWarnings(
+  suppressMessages(
+    LogicRegr(
+      formula = formula1, data = data.example, family = "Bernoulli", prior = "G",
+      report.level = 0.5, d = 15, cmax = 2, kmax = 15, p.and = 0.9, p.not = 0.2,
+      p.surv = 0.2, ncores = threads, print.freq = 0L
+    )
+  )
+)
 
 # run the inference with Jeffrey's prior
-res1J <- LogicRegr(formula = formula1, data = data.example, family = "Bernoulli", prior = "J", report.level = 0.5, d = 15, cmax = 2, kmax = 15, p.and = 0.9, p.not = 0.2, p.surv = 0.2, ncores = 32)
-print(res1J$feat.stat)
+res1J <- suppressWarnings(
+  suppressMessages(
+    LogicRegr(
+      formula = formula1, data = data.example, family = "Bernoulli", prior = "J",
+      report.level = 0.5, d = 15, cmax = 2, kmax = 15, p.and = 0.9, p.not = 0.2,
+      p.surv = 0.2, ncores = threads, print.freq = 0L
+    )
+  )
+)
+test_that("outputs are correct", {
+  expect_equal(
+    res4G$feat.stat,
+    matrix(
+    c(
+      "I(V9)", "1", "I(V5)", "1", "I(((V5))|((V9)))", "1", "I(V8)",
+      "0.999798166079656", "I(V11)", "0.998980768894478", "I(V1)",
+      "0.996827860563252", "I(V4)", "0.989626671838795"
+    ),
+    nrow = 7, ncol = 2, byrow = TRUE
+    )
+  )
+  expect_equal(
+    res4J$feat.stat,
+    matrix(
+    c(
+      "I(((V9))&((V5)))", "0.999999999999948", "I(((V8))&((V11)))",
+      "0.9999999984096", "I(((V4))&((V1)))", "0.99896404225898"
+    ),
+    nrow = 3, ncol = 2, byrow = TRUE
+    )
+  )
+  expect_equal(
+    res1G$feat.stat,
+    matrix(
+      c(
+        "I(((V9))&((V5)))", "1", "I(V8)", "0.999999999999935",
+        "I(((V1))&(1-(V4)))", "0.999986154466532",
+        "I(V11)", "0.999445448885581", "I(V1)", "0.997930154511171",
+        "I(V4)", "0.997375287784621"
+      ),
+      nrow = 6, ncol = 2, byrow = TRUE
+    )
+  )
+  expect_equal(
+    res1J$feat.stat,
+    matrix(
+      c(
+        "I(V1)", "1", "I(V4)", "1", "I(V11)", "1", "I(((V4))|((V8)))", "1",
+        "I(((V5))&((V9)))", "1"
+      ),
+      nrow = 5, ncol = 2, byrow = TRUE
+    )
+  )
+})