diff --git a/supp_mat/SM3_ExplainingDeviation.qmd b/supp_mat/SM3_ExplainingDeviation.qmd
index d177a10..27f7990 100644
--- a/supp_mat/SM3_ExplainingDeviation.qmd
+++ b/supp_mat/SM3_ExplainingDeviation.qmd
@@ -1241,21 +1241,78 @@ patchwork::wrap_plots(yi_deviation_RE_plots, byrow = TRUE) +
 ```{r}
 #| label: multivariate-models-mod-summary
 
-# lmer(box_cox_abs_deviation_score_estimate ~ RateAnalysis + PublishableAsIs + mean_diversity_index + (1|study_id) + (1|ReviewerId), data = ManyEcoEvo_results$effects_analysis[[1]] %>% unnest(review_data)) 
-
-bt_multivar_mod <- lmer(box_cox_abs_deviation_score_estimate ~ RateAnalysis + PublishableAsIs + mean_diversity_index + mixed_model + (1|ReviewerId), 
-     data = ManyEcoEvo_results$effects_analysis[[1]] %>% 
-       unnest(review_data))
+fit_multivar_wrap <- function(data_tbl, ..., env = caller_env()){
+  f1 <- rlang::new_formula(expr(box_cox_abs_deviation_score_estimate), 
+                           expr(RateAnalysis + 
+                                  PublishableAsIs + 
+                                  mean_diversity_index + 
+                                  (1|ReviewerId)), env = env)
+  
+  f2 <- rlang::new_formula(expr(box_cox_abs_deviation_score_estimate), 
+                           expr(RateAnalysis + 
+                                  PublishableAsIs + 
+                                  mean_diversity_index + 
+                                  mixed_model +
+                                  (1|ReviewerId)), env = env)
+  
+  pass_threshold <- 
+    data_tbl %>% 
+    count(mixed_model) %>% 
+    pointblank::test_col_vals_gte(n, 5)
+  
+  if(pass_threshold == TRUE) {
+    cli::cli_alert_info("Model with random effects included")
+  } else (
+    cli::cli_alert_info("Model with random effects excluded")
+  )
+  #TODO MAKE SURE GIVES CORRECT EX
+  f <- if(pass_threshold) f2 else f1 # MAKE SURE RETURNS APPROPIRATELY
+  
+  mod <- inject(lme4::lmer(!!f, data = data_tbl, ...))
+  
+  return(mod)
+  
+}
 
-euc_multivar_mod <- 
-lmer(box_cox_abs_deviation_score_estimate ~ RateAnalysis + PublishableAsIs + mean_diversity_index + mixed_model + (1|ReviewerId), 
-     data = ManyEcoEvo_results$effects_analysis[[2]] %>% 
-       unnest(review_data))
+filter_vars <- exprs(exclusion_set == "complete",
+                    expertise_subset == "All",
+                    publishable_subset == "All",
+                    collinearity_subset == "All")
 
-bt_multivar_mod_R <- bt_multivar_mod %>% MuMIn::r.squaredGLMM()
-euc_multivar_mod_R <- euc_multivar_mod %>% MuMIn::r.squaredGLMM()
-bt_multivar_mod_sigma<- bt_multivar_mod %>% sigma()
-euc_multivar_mod_sigma<- euc_multivar_mod %>% sigma()
+multivar_mods <- 
+  ManyEcoEvo_results %>% 
+  dplyr::filter(!!!filter_vars) %>% 
+  group_by(dataset) %>% 
+  select(dataset, effects_analysis) %>% 
+  mutate(effects_analysis = map(effects_analysis, ~ .x %>% 
+                                  unnest(review_data))) %>% 
+  mutate(model = map(effects_analysis, fit_multivar_wrap)) %>%
+  mutate(model_perf = map(model, performance::model_performance)) %>% 
+  mutate(model_params = map(model, parameters::parameters)) %>% 
+  hoist(model_perf, "R2_conditional", "R2_marginal", "Sigma")
+
+bt_multivar_mod_R <- 
+  multivar_mods %>% ungroup %>%  
+  filter(dataset == "blue tit") %>% 
+  select(R2_marginal, R2_conditional) %>%  
+  transpose() %>% 
+  flatten_dbl()
+
+euc_multivar_mod_R <- 
+  multivar_mods %>% 
+  ungroup %>%  
+  filter(dataset == "eucalyptus") %>% 
+  select(R2_marginal, R2_conditional) %>%  
+  transpose() %>% 
+  flatten_dbl()
+
+bt_multivar_mod_sigma <- multivar_mods %>% 
+  filter(dataset == "blue tit") %>% 
+  pluck("Sigma")
+
+euc_multivar_mod_sigma <- multivar_mods %>% 
+  filter(dataset == "eucalyptus") %>% 
+  pluck("Sigma")
 ```
 
 ```{r}
@@ -1264,13 +1321,12 @@ euc_multivar_mod_sigma<- euc_multivar_mod %>% sigma()
 #| message: false
 #| column: page-right
 #| tbl-cap: "Parameter estimates from models explaining Box-Cox transformed deviation scores from the mean $Z_r$ as a function of continuous and categorical peer-review ratings in multivariate analyses. Standard Errors (SE), 95% Confidence Intervals (95%CI) are reported for all estimates, while t values, degrees of freedom and p-values are presented for fixed-effects."
-list(bt_multivar_mod, euc_multivar_mod) %>% 
-  set_names(c("blue tit", "eucalyptus")) %>% 
-  map_dfr(parameters::parameters, .id = "dataset") %>% 
-  group_by(dataset) %>% 
-  select(-CI) %>% 
-  mutate(dataset = str_replace(dataset, "eucalyptus", "*Eucalyptus*"),
-         Parameter = str_replace(Parameter, "mixed_model", "random_included")) %>% 
+
+
+multivar_mods %>% select(dataset, model_params) %>% unnest(model_params) %>% select(-CI) %>% 
+  mutate(
+    dataset = str_replace(dataset, "eucalyptus", "*Eucalyptus*"),             
+    Parameter = str_replace(Parameter, "mixed_model", "random_included"))  %>% 
   gt::gt() %>% 
   gt::fmt(columns = "p",
           fns = function(x) gtsummary::style_pvalue(x, prepend_p = TRUE)) %>% 
@@ -1290,7 +1346,7 @@ list(bt_multivar_mod, euc_multivar_mod) %>%
   },
   locations = cells_body(columns = Parameter)) %>% 
   gt::text_transform(fn = function(x) str_replace(x, "ReviewerId", "Reviewer ID")) %>% 
-    gt::text_transform(fn = function(x) map(x, gt::md), 
+  gt::text_transform(fn = function(x) map(x, gt::md), 
                      locations = gt::cells_row_groups()) %>% 
   gt::sub_missing(missing_text = "")