Rshiny app roughly drafted + improving effects output by including th…

…e intercept and other fixed effects in the Bayesian models
jeffersonfparil · Jun 13, 2024 · 23a8b41 · 23a8b41
1 parent f252478
commit 23a8b41
Show file tree

Hide file tree

Showing 5 changed files with 1,154 additions and 13 deletions.
diff --git a/R/cross_validation.R b/R/cross_validation.R
@@ -770,6 +770,18 @@ fn_cross_validation_within_population = function(list_merged, n_folds=10, n_reps
         df_y_validation = NULL
         for (list_perf in list_list_perf) {
             # list_perf = list_list_perf[[1]]
+            if (methods::is(list_perf, "gpError")) {
+                error = chain(list_perf, methods::new("gpError",
+                    code=000,
+                    message=paste0(
+                        "Error in cross_validation::fn_cross_validation_within_population(...). ",
+                        "Error running cross-validation for population: ", population, " at ",
+                        "rep: ", list_cv_params$df_params$rep[i], ", ",
+                        "fold: ", list_cv_params$df_params$fold[i], ", and ",
+                        "model: ", list_cv_params$df_params$model[i], "."
+                    )))
+                return(error)
+            }
             if (is.null(df_metrics) & is.null(df_y_validation)) {
                 df_metrics = list_perf$df_metrics
                 df_y_validation = list_perf$df_y_validation

diff --git a/R/models.R b/R/models.R
@@ -722,14 +722,19 @@ fn_Bayes_A = function(list_merged, vec_idx_training, vec_idx_validation,
         ETA = list(MRK=list(X=list_merged$G, model="BayesA"),
                        list(X=list_merged$COVAR, model="FIXED"))
     } else {
-        ETA = list(MRK=list(X=list_merged$G, model="BayesA"))
+        ETA = list(MRK=list(X=list_merged$G, model="BayesA", saveEffects=TRUE))
     }
     ### Attempt at preventing overwrites to the running Gibbs samplers in parallel
     other_params$out_prefix = gsub(":", ".", gsub(" ", "-", paste(other_params$out_prefix, Sys.time(), stats::runif(1), sep="-")))
     ### Solve via Bayes A (a priori assume heritability at 50%, i.e. R2=0.5 below)
     sol = BGLR::BGLR(y=yNA, ETA=ETA, R2=0.5, nIter=other_params$nIter, burnIn=other_params$burnIn, saveAt=other_params$out_prefix, verbose=FALSE)
-    ### Extract effects
-    b_hat = sol$ETA$MRK$b
+    ### Extract effects including the intercept and fixed effects
+    if (!is.null(list_merged$COVAR)) {
+        b_hat = c(sol$mu, sol$ETA[[2]]$b, sol$ETA$MRK$b)
+    } else {
+        b_hat = c(sol$mu, sol$ETA$MRK$b)
+    }
+    names(b_hat)[1] = "intercept"
     n_non_zero = sum(abs(b_hat) >= .Machine$double.eps)
     if (verbose) {
         p = ncol(list_merged$G)
@@ -870,8 +875,13 @@ fn_Bayes_B = function(list_merged, vec_idx_training, vec_idx_validation,
     other_params$out_prefix = gsub(":", ".", gsub(" ", "-", paste(other_params$out_prefix, Sys.time(), stats::runif(1), sep="-")))
     ### Solve via Bayes B (a priori assume heritability at 50%, i.e. R2=0.5 below)
     sol = BGLR::BGLR(y=yNA, ETA=ETA, R2=0.5, nIter=other_params$nIter, burnIn=other_params$burnIn, saveAt=other_params$out_prefix, verbose=FALSE)
-    ### Extract effects
-    b_hat = sol$ETA$MRK$b
+    ### Extract effects including the intercept and fixed effects
+    if (!is.null(list_merged$COVAR)) {
+        b_hat = c(sol$mu, sol$ETA[[2]]$b, sol$ETA$MRK$b)
+    } else {
+        b_hat = c(sol$mu, sol$ETA$MRK$b)
+    }
+    names(b_hat)[1] = "intercept"
     n_non_zero = sum(abs(b_hat) >= .Machine$double.eps)
     if (verbose) {
         p = ncol(list_merged$G)
@@ -1012,8 +1022,13 @@ fn_Bayes_C = function(list_merged, vec_idx_training, vec_idx_validation,
     other_params$out_prefix = gsub(":", ".", gsub(" ", "-", paste(other_params$out_prefix, Sys.time(), stats::runif(1), sep="-")))
     ### Solve via Bayes C (a priori assume heritability at 50%, i.e. R2=0.5 below)
     sol = BGLR::BGLR(y=yNA, ETA=ETA, R2=0.5, nIter=other_params$nIter, burnIn=other_params$burnIn, saveAt=other_params$out_prefix, verbose=FALSE)
-    ### Extract effects
-    b_hat = sol$ETA$MRK$b
+    ### Extract effects including the intercept and fixed effects
+    if (!is.null(list_merged$COVAR)) {
+        b_hat = c(sol$mu, sol$ETA[[2]]$b, sol$ETA$MRK$b)
+    } else {
+        b_hat = c(sol$mu, sol$ETA$MRK$b)
+    }
+    names(b_hat)[1] = "intercept"
     n_non_zero = sum(abs(b_hat) >= .Machine$double.eps)
     if (verbose) {
         p = ncol(list_merged$G)