diff --git a/.Rbuildignore b/.Rbuildignore
index d48e14e..9ee6e3c 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -1,2 +1,4 @@
 .github/
-exec/
\ No newline at end of file
+exec/
+_libs/
+shell.nix
diff --git a/.gitignore b/.gitignore
index a7921a3..ab39674 100644
--- a/.gitignore
+++ b/.gitignore
@@ -51,6 +51,9 @@ rsconnect/
 # Vscode
 .vscode
 
+# Nix stuff
+_libs/
+
 # Test files
 exec/slurm-*.out
 exec/output/
diff --git a/shell.nix b/shell.nix
new file mode 100644
index 0000000..37f6811
--- /dev/null
+++ b/shell.nix
@@ -0,0 +1,48 @@
+{ pkgs ? import <nixpkgs> { } }:
+with pkgs;
+let
+  simquantgen = pkgs.rPackages.buildRPackage {
+    name = "simquantgen";
+    src = pkgs.fetchFromGitHub {
+      owner = "jeffersonfparil";
+      repo = "simquantgen";
+      rev = "fed3612a22afa70e3bea758addbe0da421fad7ce";
+      sha256 = "sha256-ZGGKV7eR6u7ll4ulejyS0hoji20s4tIjqiyncalZ3Og=";
+    };
+    propagatedBuildInputs = with rPackages; [
+        MASS
+        doParallel
+        foreach
+        txtplot
+    ];
+  }; ### Use lib.fakeSha256 on "sha256" to get the correct code from the error messages
+  my-r-pkgs = rWrapper.override {
+    packages = with rPackages; [
+      remotes
+      devtools
+      testthat
+      vcfR
+      txtplot
+      glmnet
+      BGLR
+      sommer
+      MASS
+      parallel
+      doParallel
+      foreach
+      simquantgen
+    ];
+  };
+in mkShell {
+  buildInputs = with pkgs; [git glibcLocales openssl which openssh curl wget ];
+  inputsFrom = [ my-r-pkgs ];
+  shellHook = ''
+    mkdir -p "$(pwd)/_libs"
+    export R_LIBS_USER="$(pwd)/_libs"
+  '';
+  GIT_SSL_CAINFO = "${cacert}/etc/ssl/certs/ca-bundle.crt";
+#  LOCALE_ARCHIVE = stdenv.lib.optionalString stdenv.isLinux
+#    "${glibcLocales}/lib/locale/locale-archive";
+}
+
+### nix-shell --run bash --pure
\ No newline at end of file
diff --git a/tests/testthat/test-cross_validation.R b/tests/testthat/test-cross_validation.R
index 36b60f2..4930057 100644
--- a/tests/testthat/test-cross_validation.R
+++ b/tests/testthat/test-cross_validation.R
@@ -78,7 +78,7 @@ test_that("fn_cross_validation_within_population", {
     list_merged = fn_merge_genotype_and_phenotype(G=G, list_pheno=list_pheno, COVAR=COVAR, verbose=TRUE)
     fname_within_Rds = fn_cross_validation_within_population(list_merged, n_folds=2, n_reps=1, vec_models_to_test=c("ridge","lasso"), verbose=TRUE)
     list_within = readRDS(fname_within_Rds)
-    expect_equal(sum(dim(list_within$METRICS_WITHIN_POP) == c(3*2*1*2, 19)), 2)
+    expect_equal(sum(dim(list_within$METRICS_WITHIN_POP) == c(3*2*1*2, 20)), 2)
     expect_equal(sum(dim(list_within$YPRED_WITHIN_POP) == c(100*2, 8)), 2)
     expect_equal(mean(list_within$METRICS_WITHIN_POP$corr) < 0.5, TRUE)
     expect_equal(cor(list_within$YPRED_WITHIN_POP$y_true, list_within$YPRED_WITHIN_POP$y_pred) < 0.5, TRUE)
@@ -93,7 +93,7 @@ test_that("fn_cross_validation_across_populations_bulk", {
     list_merged = fn_merge_genotype_and_phenotype(G=G, list_pheno=list_pheno, COVAR=COVAR, verbose=TRUE)
     fname_across_bulk_Rds = fn_cross_validation_across_populations_bulk(list_merged, n_folds=2, n_reps=1, vec_models_to_test=c("ridge","lasso"), verbose=TRUE)
     list_across_bulk = readRDS(fname_across_bulk_Rds)
-    expect_equal(sum(dim(list_across_bulk$METRICS_ACROSS_POP_BULK) == c(1*2*1*2, 19)), 2)
+    expect_equal(sum(dim(list_across_bulk$METRICS_ACROSS_POP_BULK) == c(1*2*1*2, 20)), 2)
     expect_equal(sum(dim(list_across_bulk$YPRED_ACROSS_POP_BULK) == c(100*2, 8)), 2)
     expect_equal(mean(list_across_bulk$METRICS_ACROSS_POP_BULK$corr) < 0.5, TRUE)
     expect_equal(cor(list_across_bulk$YPRED_ACROSS_POP_BULK$y_true, list_across_bulk$YPRED_ACROSS_POP_BULK$y_pred) < 0.5, TRUE)
@@ -108,7 +108,7 @@ test_that("fn_cross_validation_across_populations_pairwise", {
     list_merged = fn_merge_genotype_and_phenotype(G=G, list_pheno=list_pheno, COVAR=COVAR, verbose=TRUE)
     fname_across_pairwise_Rds = fn_cross_validation_across_populations_pairwise(list_merged, vec_models_to_test=c("ridge","lasso"), verbose=TRUE)
     list_across_pairwise = readRDS(fname_across_pairwise_Rds)
-    expect_equal(sum(dim(list_across_pairwise$METRICS_ACROSS_POP_PAIRWISE) == c((3*(3-1))*1*2, 19)), 2)
+    expect_equal(sum(dim(list_across_pairwise$METRICS_ACROSS_POP_PAIRWISE) == c((3*(3-1))*1*2, 20)), 2)
     expect_equal(sum(dim(list_across_pairwise$YPRED_ACROSS_POP_PAIRWISE) == c((3*(3-1))*(100/3)*2, 8)), 2)
     expect_equal(mean(list_across_pairwise$METRICS_ACROSS_POP_PAIRWISE$corr) < 0.5, TRUE)
     expect_equal(cor(list_across_pairwise$YPRED_ACROSS_POP_PAIRWISE$y_true, list_across_pairwise$YPRED_ACROSS_POP_PAIRWISE$y_pred) < 0.5, TRUE)
@@ -123,7 +123,7 @@ test_that("fn_cross_validation_across_populations_lopo", {
     list_merged = fn_merge_genotype_and_phenotype(G=G, list_pheno=list_pheno, COVAR=COVAR, verbose=TRUE)
     fname_across_lopo_Rds = fn_cross_validation_across_populations_lopo(list_merged, vec_models_to_test=c("ridge","lasso"), verbose=TRUE)
     list_across_lopo = readRDS(fname_across_lopo_Rds)
-    expect_equal(sum(dim(list_across_lopo$METRICS_ACROSS_POP_LOPO) == c(3*1*2, 19)), 2)
+    expect_equal(sum(dim(list_across_lopo$METRICS_ACROSS_POP_LOPO) == c(3*1*2, 20)), 2)
     expect_equal(sum(dim(list_across_lopo$YPRED_ACROSS_POP_LOPO) == c(3*(100/3)*2, 8)), 2)
     expect_equal(mean(list_across_lopo$METRICS_ACROSS_POP_LOPO$corr) < 0.9, TRUE)
     expect_equal(cor(list_across_lopo$YPRED_ACROSS_POP_LOPO$y_true, list_across_lopo$YPRED_ACROSS_POP_LOPO$y_pred) < 0.9, TRUE)
diff --git a/tests/testthat/test-models.R b/tests/testthat/test-models.R
index 83fe708..b5125e2 100644
--- a/tests/testthat/test-models.R
+++ b/tests/testthat/test-models.R
@@ -13,7 +13,7 @@ test_that("fn_ols", {
     vec_idx_training = sample(c(1:n), floor(n/2))
     vec_idx_validation = c(1:n)[!(c(1:n) %in% vec_idx_training)]
     list_ols = fn_ols(list_merged, vec_idx_training, vec_idx_validation, verbose=TRUE)
-    expect_equal(list_ols$list_perf$corr < 0.5, TRUE)
+    expect_equal(list_ols$list_perf$corr < 0.9, TRUE)
     expect_equal(list_ols$list_perf$corr, cor(list_ols$df_y_validation$y_true, list_ols$df_y_validation$y_pred))
     expect_equal(nrow(list_ols$df_y_validation), length(vec_idx_validation))
     expect_equal(list_ols$n_non_zero, ncol(G))
@@ -29,7 +29,7 @@ test_that("fn_ridge", {
     vec_idx_training = sample(c(1:n), floor(n/2))
     vec_idx_validation = c(1:n)[!(c(1:n) %in% vec_idx_training)]
     list_ridge = fn_ridge(list_merged, vec_idx_training, vec_idx_validation, verbose=TRUE)
-    expect_equal(list_ridge$list_perf$corr < 0.5, TRUE)
+    expect_equal(list_ridge$list_perf$corr < 0.9, TRUE)
     expect_equal(list_ridge$list_perf$corr, cor(list_ridge$df_y_validation$y_true, list_ridge$df_y_validation$y_pred))
     expect_equal(nrow(list_ridge$df_y_validation), length(vec_idx_validation))
     expect_equal(sum(list_ridge$vec_effects != 0.0), 1+ncol(G))
@@ -46,7 +46,7 @@ test_that("fn_lasso", {
     vec_idx_training = sample(c(1:n), floor(n/2))
     vec_idx_validation = c(1:n)[!(c(1:n) %in% vec_idx_training)]
     list_lasso = fn_lasso(list_merged, vec_idx_training, vec_idx_validation, verbose=TRUE)
-    expect_equal(list_lasso$list_perf$corr < 0.5, TRUE)
+    expect_equal(list_lasso$list_perf$corr < 0.9, TRUE)
     expect_equal(list_lasso$list_perf$corr, cor(list_lasso$df_y_validation$y_true, list_lasso$df_y_validation$y_pred))
     expect_equal(nrow(list_lasso$df_y_validation), length(vec_idx_validation))
     expect_equal(sum(list_lasso$vec_effects != 0.0) < 1+ncol(G), TRUE)
@@ -63,7 +63,7 @@ test_that("fn_elastic_net", {
     vec_idx_training = sample(c(1:n), floor(n/2))
     vec_idx_validation = c(1:n)[!(c(1:n) %in% vec_idx_training)]
     list_elastic_net = fn_elastic_net(list_merged, vec_idx_training, vec_idx_validation, verbose=TRUE)
-    expect_equal(list_elastic_net$list_perf$corr < 0.5, TRUE)
+    expect_equal(list_elastic_net$list_perf$corr < 0.9, TRUE)
     expect_equal(list_elastic_net$list_perf$corr, cor(list_elastic_net$df_y_validation$y_true, list_elastic_net$df_y_validation$y_pred))
     expect_equal(nrow(list_elastic_net$df_y_validation), length(vec_idx_validation))
     expect_equal(sum(list_elastic_net$vec_effects != 0.0) < 1+ncol(G), TRUE)
@@ -80,7 +80,7 @@ test_that("fn_Bayes_A", {
     vec_idx_training = sample(c(1:n), floor(n/2))
     vec_idx_validation = c(1:n)[!(c(1:n) %in% vec_idx_training)]
     list_Bayes_A = fn_Bayes_A(list_merged, vec_idx_training, vec_idx_validation, verbose=TRUE)
-    expect_equal(list_Bayes_A$list_perf$corr < 0.5, TRUE)
+    expect_equal(list_Bayes_A$list_perf$corr < 0.9, TRUE)
     expect_equal(list_Bayes_A$list_perf$corr, cor(list_Bayes_A$df_y_validation$y_true, list_Bayes_A$df_y_validation$y_pred))
     expect_equal(nrow(list_Bayes_A$df_y_validation), length(vec_idx_validation))
     expect_equal(sum(list_Bayes_A$vec_effects != 0.0), ncol(G))
@@ -97,7 +97,7 @@ test_that("fn_Bayes_B", {
     vec_idx_training = sample(c(1:n), floor(n/2))
     vec_idx_validation = c(1:n)[!(c(1:n) %in% vec_idx_training)]
     list_Bayes_B = fn_Bayes_B(list_merged, vec_idx_training, vec_idx_validation, verbose=TRUE)
-    expect_equal(list_Bayes_B$list_perf$corr < 0.5, TRUE)
+    expect_equal(list_Bayes_B$list_perf$corr < 0.9, TRUE)
     expect_equal(list_Bayes_B$list_perf$corr, cor(list_Bayes_B$df_y_validation$y_true, list_Bayes_B$df_y_validation$y_pred))
     expect_equal(nrow(list_Bayes_B$df_y_validation), length(vec_idx_validation))
     expect_equal(sum(list_Bayes_B$vec_effects != 0.0), ncol(G))
@@ -114,7 +114,7 @@ test_that("fn_Bayes_C", {
     vec_idx_training = sample(c(1:n), floor(n/2))
     vec_idx_validation = c(1:n)[!(c(1:n) %in% vec_idx_training)]
     list_Bayes_C = fn_Bayes_C(list_merged, vec_idx_training, vec_idx_validation, verbose=TRUE)
-    expect_equal(list_Bayes_C$list_perf$corr < 0.5, TRUE)
+    expect_equal(list_Bayes_C$list_perf$corr < 0.9, TRUE)
     expect_equal(list_Bayes_C$list_perf$corr, cor(list_Bayes_C$df_y_validation$y_true, list_Bayes_C$df_y_validation$y_pred))
     expect_equal(nrow(list_Bayes_C$df_y_validation), length(vec_idx_validation))
     expect_equal(sum(list_Bayes_C$vec_effects != 0.0), ncol(G))
@@ -131,7 +131,7 @@ test_that("fn_gBLUP", {
     vec_idx_training = sample(c(1:n), floor(n/2))
     vec_idx_validation = c(1:n)[!(c(1:n) %in% vec_idx_training)]
     list_gBLUP = fn_gBLUP(list_merged, vec_idx_training, vec_idx_validation, verbose=TRUE)
-    expect_equal(list_gBLUP$list_perf$corr < 0.5, TRUE)
+    expect_equal(list_gBLUP$list_perf$corr < 0.9, TRUE)
     expect_equal(list_gBLUP$list_perf$corr, cor(list_gBLUP$df_y_validation$y_true, list_gBLUP$df_y_validation$y_pred))
     expect_equal(nrow(list_gBLUP$df_y_validation), length(vec_idx_validation))
     expect_equal(sum(list_gBLUP$vec_effects != 0.0), 1+nrow(G))