Merge pull request #92 from egouldo/67-model-params

Extract model params from supported models

.Rbuildignore

@@ -4,6 +4,7 @@
 ^.*\.Rproj$
 ^\.Rproj\.user$
 _targets.R
+_targets/
 
 ^_pkgdown\.yml$
 ^docs$

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: ManyEcoEvo
 Title: Meta-analyse data from 'Many-Analysts' style studies
-Version: 2.0.0.9001
+Version: 2.3.0
 Authors@R: c(person(given = "Elliot",
                     family = "Gould", 
                     email =  "[email protected]", 
@@ -32,22 +32,23 @@ Imports:
     tidyselect,
     pointblank,
     tibble,
-    cli, 
-    data.table, 
-    forcats, 
-    fs, 
-    glue, 
-    here, 
-    lme4, 
-    metafor, 
+    cli,
+    data.table,
+    forcats,
+    fs,
+    glue,
+    here,
+    lme4,
+    metafor,
     naniar,
     magrittr,
     tidyr,
     rlang,
     purrr,
-    tidyselect,
+    tidyselect
 Suggests:
     targets,
+    qs,
     broom.mixed,
     metaviz,
     ggeffects,
@@ -66,7 +67,7 @@ Remotes:
   daniel1noble/orchaRd,
   NightingaleHealth/ggforestplot
 Encoding: UTF-8
-LazyData: true
+LazyData: false
 Roxygen: list(markdown = TRUE)
 RoxygenNote: 7.3.2
 URL: https://github.com/egouldo/ManyEcoEvo,

NAMESPACE

@@ -42,6 +42,7 @@ export(fit_boxcox_ratings_ord)
 export(fit_metafor_mv)
 export(fit_metafor_mv_reduced)
 export(fit_metafor_uni)
+export(fit_multivar_MA)
 export(fit_sorensen_glm)
 export(fit_uni_mixed_effects)
 export(generate_collinearity_subset)
@@ -108,12 +109,14 @@ import(cli)
 import(dplyr)
 import(ggbeeswarm)
 import(ggplot2)
+import(lme4)
 import(metafor)
 import(recipes)
 import(rlang)
 import(see)
 importFrom(EnvStats,stat_n_text)
 importFrom(broom,tidy)
+importFrom(broom.mixed,tidy)
 importFrom(cli,cli_abort)
 importFrom(cli,cli_alert_info)
 importFrom(cli,cli_alert_warning)
@@ -136,28 +139,39 @@ importFrom(dplyr,right_join)
 importFrom(dplyr,select)
 importFrom(dplyr,summarise)
 importFrom(forcats,fct_relevel)
+importFrom(ggplot2,ggplot)
 importFrom(glue,glue)
+importFrom(lme4,lmer)
 importFrom(magrittr,"%>%")
+importFrom(metaviz,viz_funnel)
+importFrom(parameters,parameters)
+importFrom(performance,performance)
 importFrom(pointblank,col_vals_not_null)
 importFrom(pointblank,stop_if_not)
+importFrom(pointblank,test_col_vals_gte)
 importFrom(purrr,keep)
 importFrom(purrr,list_flatten)
 importFrom(purrr,list_rbind)
 importFrom(purrr,map)
 importFrom(purrr,map2)
 importFrom(purrr,map_chr)
 importFrom(purrr,map_dfr)
+importFrom(purrr,map_if)
 importFrom(purrr,pluck)
 importFrom(purrr,pmap)
 importFrom(purrr,possibly)
 importFrom(purrr,reduce)
 importFrom(purrr,reduce2)
 importFrom(purrr,set_names)
+importFrom(rlang,caller_env)
 importFrom(rlang,enquo)
+importFrom(rlang,expr)
 importFrom(rlang,is_na)
 importFrom(rlang,is_null)
 importFrom(rlang,na_chr)
+importFrom(rlang,new_formula)
 importFrom(sae,bxcx)
+importFrom(stringr,str_detect)
 importFrom(tibble,enframe)
 importFrom(tibble,tibble)
 importFrom(tidyr,hoist)

R/Z_VZ_preds.R

@@ -1,4 +1,4 @@
-#' Calculate Z and VZ of out-of-sample predictions
+#' Standardize out-of-sample predictions
 #'
 #' @param yi point-estimate prediction, on response scale
 #' @param yi_se standard error of \code{yi}
@@ -11,7 +11,7 @@
 Z_VZ_preds <- function(yi, yi_se, mu_p, sd_p ){
   #TODO should we pass in whole DF as arg instead of yi + yi_se??
   # We want to be able to keep the values linked to their corresponding
-  # scenario_ value!
+  # scenario_value!
   na_args <- purrr::discard(c(yi, yi_se, mu_p, sd_p), is.na) %>% 
     length()
 
@@ -26,8 +26,8 @@ Z_VZ_preds <- function(yi, yi_se, mu_p, sd_p ){
   }
 
 
-  Z <- (yi-mu_p)/sd_p
-  VZ <- yi_se/sd_p
+  Z <- (yi - mu_p) / sd_p
+  VZ <- yi_se / sd_p
 
   return(tibble(Z, VZ)) 
 }

R/calculate_descriptive_statistics.R

@@ -236,24 +236,24 @@ summarise_study <- function(ManyEcoEvo, ManyEcoEvo_results, id_subsets, subset_n
 
   # ------ Prepare Summary Data Subsets ------
 
-  subsets_tibble <- ManyEcoEvo::ManyEcoEvo %>% 
+  subsets_tibble <- ManyEcoEvo %>% 
     prepare_analyst_summary_data("all",
                                  id_subsets,
                                  subset_names)
 
-  subsets_tibble_sorensen <- ManyEcoEvo::ManyEcoEvo_results %>% 
+  subsets_tibble_sorensen <- ManyEcoEvo_results %>% 
     prepare_sorenson_summary_data("all",
                                   id_subsets,
                                   subset_names,
                                   filter_expressions = filter_vars)
 
-  subsets_tibble_variables <- ManyEcoEvo::ManyEcoEvo %>%
+  subsets_tibble_variables <- ManyEcoEvo %>%
     prepare_diversity_summary_data("all",
                                    id_subsets,
                                    subset_names)
 
   var_names <- 
-    ManyEcoEvo::ManyEcoEvo %>% 
+    ManyEcoEvo %>% 
     pull(diversity_data) %>% 
     map(~ .x %>% select(-id_col, -dataset) 
         %>% colnames) %>% 

R/filt_multivar_MA.R

@@ -0,0 +1,89 @@
+#' Fit a multivariate meta-regression model
+#' 
+#' @description Fit a multivariate meta-regression model that models the effect of peer-review ratings on the deviation from the meta-analytic mean (both continuous and categorical ratings), mean Sorensen's index, and/or whether the analysis uses a mixed effects model, or not.
+#' 
+#' @param data_tbl Data for model fitting
+#' @param ... Additional arguments passed to `lmer`
+#' @param env Environment in which to evaluate the formula, defaults to the calling environment
+#' @param N threshold for the number of analyses that must have been conducted using mixed effects models to include the binary predictor `mixed_model` in the meta-regression. Defaults to 5.
+#' 
+#' @return An object of class lmer.
+#' 
+#' @export
+#' @importFrom rlang new_formula
+#' @importFrom rlang caller_env
+#' @importFrom rlang expr
+#' @importFrom rlang expr
+#' @importFrom lme4 lmer
+#' @importFrom pointblank test_col_vals_gte
+#' @import dplyr
+#' @import cli
+#' @details
+#' Depending on whether enough analyses in `data_tbl` have been conducted with the `mixed_model` variable, the function will fit a model with or without the predictor `mixed_model`.
+#' 
+#' Expects the following columns in `data_tbl`:
+#' 
+#' - `RateAnalysis`: continuous peer-review ratings
+#' - `PublishableAsIs`: categorical peer-review ratings
+#' - `mean_diversity_index`: mean Sorensen's index
+#' - `box_cox_abs_deviation_score_estimate`: response variable, Box-Cox transformed deviation from the meta-analytic mean effect-size for each analysis
+#' - `mixed_model`: binary variable indicating whether the analysis used a mixed effects model or not
+#' - `ReviewerId`: reviewer identifier
+fit_multivar_MA <- function(data_tbl, N = 5, ..., env = rlang::caller_env()){
+
+  data_tbl %>% 
+    pointblank::expect_col_exists(columns = c(box_cox_abs_deviation_score_estimate, 
+                                              RateAnalysis, PublishableAsIs,
+                                              mean_diversity_index, 
+                                              ReviewerId, 
+                                              mixed_model))
+  # Define Models
+  f1 <- rlang::new_formula(rlang::expr(box_cox_abs_deviation_score_estimate), 
+                           rlang::expr(RateAnalysis + 
+                                         PublishableAsIs + 
+                                         mean_diversity_index + 
+                                         (1 | ReviewerId)), env = env)
+
+  f2 <- rlang::new_formula(rlang::expr(box_cox_abs_deviation_score_estimate), 
+                           rlang::expr(RateAnalysis + 
+                                         PublishableAsIs + 
+                                         mean_diversity_index + 
+                                         mixed_model +
+                                         (1 | ReviewerId)), env = env)
+
+  cli::cli_h2("Fitting multivariate meta-regression model")
+
+  pass_threshold <- 
+    data_tbl %>% 
+    count(mixed_model) %>% 
+    pointblank::test_col_vals_gte(n, N)
+
+  cur_group_bullets <- dplyr::cur_group() %>% 
+    transpose() %>% 
+    list_flatten() %>% 
+    enframe() %>%  
+    mutate(value = list_c(value)) %>% 
+    unite(group, everything(), 
+          sep = ": ") %>% 
+    pull(group)
+
+  if (pass_threshold == TRUE) {
+    cli::cli_alert_info(glue::glue("Presence of random effects in analyses ", 
+                                   cli::style_italic("included"), 
+                                   " as predictor in model for data subset:"))
+    cli::cli_bullets(c(setNames(cur_group_bullets, rep("*",length(cur_group_bullets)))))
+  } else {
+    cli::cli_alert_info(glue::glue("Presence of random effects in analyses ", 
+                                   cli::style_italic("excluded"), 
+                                   " as predictor in model for data subset:"))
+    cli::cli_bullets(c(setNames(cur_group_bullets, rep("*",length(cur_group_bullets)))))
+  }
+
+  #TODO MAKE SURE GIVES CORRECT EX
+  f <- if (pass_threshold) f2 else f1 # MAKE SURE RETURNS APPROPIRATELY
+
+  mod <- rlang::inject(lme4::lmer(!!f, data = data_tbl, ...))
+
+  return(mod)
+
+}

R/fit_MA_mv.R

@@ -19,13 +19,13 @@
 #' @examples
 #' ManyEcoEvo_results$effects_analysis[[1]] %>% 
 #'   fit_MA_mv(beta_estimate, beta_SE, "Zr")
-fit_MA_mv <- function(effects_analysis = data.frame(), Z_colname, VZr_colname, estimate_type = character(1L)){
+fit_MA_mv <- function(effects_analysis = data.frame(), Z_colname, VZ_colname, estimate_type = character(1L)){
   pointblank::stop_if_not(estimate_type %in% c("Zr", "yi", "y25", "y50", "y75"))
 
-  Zr <- effects_analysis %>%  dplyr::pull({{Z_colname}})
-  VZr <- effects_analysis %>%  dplyr::pull({{VZr_colname}})
-  mod <- ManyEcoEvo::fit_metafor_mv(estimate = Zr,
-                                    variance = VZr, 
+  Z <- effects_analysis %>%  dplyr::pull({{Z_colname}})
+  VZ <- effects_analysis %>%  dplyr::pull({{VZ_colname}})
+  mod <- ManyEcoEvo::fit_metafor_mv(estimate = Z,
+                                    variance = VZ, 
                                     estimate_type = estimate_type, 
                                     data = effects_analysis)
   return(mod)

R/fit_boxcox_ratings_cont.R

@@ -18,7 +18,7 @@
 #' #   fit_boxcox_ratings_cont(.,
 #' #                                   box_cox_abs_deviation_score_estimate,
 #' #                                   VZr )
-fit_boxcox_ratings_cont <- function(.data, outcome, outcome_var) {
+fit_boxcox_ratings_cont <- function(.data, outcome, outcome_var, ..., env = rlang::caller_env()) {
   cli::cli_h2(glue::glue("Fitting metaregression with continuous ratings predictor on box_cox_transformed outcomes"))
 
   # TODO @egouldo stopifnot data doesn't contain variables named eval(box_cox_outcome_var), eval(sampling_variance_var), review_data
@@ -38,15 +38,14 @@ fit_boxcox_ratings_cont <- function(.data, outcome, outcome_var) {
   f <- rlang::new_formula(rlang::ensym(outcome), 
                           expr(RateAnalysis + 
                                  (1 | study_id) #NOTE: ReviewerId removed due to singularity
-                          ))
-  mod <- lme4::lmer(f,
-                    data = data_tbl,
-                    weights = I(1/pull(data_tbl,{{outcome_var}})))
+                          ) , env = env)
 
+  mod <- rlang::inject(lme4::lmer(!!f, data = data_tbl, ...))
+
   return(mod)
 
 }
 
 poss_fit_boxcox_ratings_cont <- purrr::possibly(fit_boxcox_ratings_cont, 
                                                 otherwise = NA,
-                                                quiet = FALSE) #TODO export
+                                                quiet = FALSE) #TODO export

R/make_viz.R

@@ -1,12 +1,25 @@
 #' Make visualisations wrapper function
 #' @description Computes model summaries, tidy model summaries, model fit stats, funnel plots and forest plots for a dataframe of multiple fitted models
 #'
-#' @param data a nested dataframe with processed and standardised data stored in list-column `data`, grouped by variables `exclusion_set`, `dataset`, `estimate_type`
+#' @param data a nested dataframe with processed and standardised data stored in list-column `data`, grouped by variables `exclusion_set`, `dataset`, `estimate_type`, `publishable_subset`, `expertise_subset`, `collinearity_subset`. Each group contains a list-column `model` containing fitted models of class `rma.uni`, `rma.mv` or `merMod`.
 #'
-#' @return a nested dataframe grouped by variables `exclusion_set`, `dataset`, `estimate_type`, 
+#' @return a nested dataframe grouped by variables `exclusion_set`, `dataset`, `estimate_type`, `publishable_subset`, `expertise_subset`, `collinearity_subset` containing model summaries, tidy model summaries, model fit stats, funnel plots and forest plots
 #' @export
 #' @family targets-pipeline functions
 #' @family Multi-dataset Wrapper Functions
+#' @import dplyr
+#' @importFrom purrr map_if map2 pmap possibly
+#' @importFrom stringr str_detect
+#' @importFrom broom.mixed tidy
+#' @importFrom performance performance
+#' @importFrom metaviz viz_funnel
+#' @importFrom ggplot2 ggplot
+#' @importFrom parameters parameters
+#' @import metafor
+#' @import lme4
+#' @importFrom tidyr pivot_longer
+#' @importFrom tidyr unnest
+#' @importFrom rlang is_na
 make_viz <- function(data) {
   # targets wrapper function
   # define map helper fun
@@ -15,7 +28,7 @@ make_viz <- function(data) {
   }
   # remove unnecessary inputs, create summary tables and visualisations
   # repeat for yi and Zr
-  if(any(str_detect(unique(data$estimate_type),pattern = "Zr"))){
+  if (any(str_detect(unique(data$estimate_type),pattern = "Zr"))) {
     data_Zr <- data %>% 
       filter(estimate_type == "Zr") %>% 
       group_by(exclusion_set, dataset, estimate_type, publishable_subset, expertise_subset, collinearity_subset, data) %>% 
@@ -59,13 +72,14 @@ make_viz <- function(data) {
       mutate(publishable_subset = NA)
   }
 
-  if(exists("data_Zr") & exists("data_yi")){
+  if (exists("data_Zr") & exists("data_yi")) {
     all_data <- bind_rows(data_Zr, data_yi)
   } else if (exists("data_Zr")) {
     all_data <- data_Zr
   } else {
     all_data <- data_yi
   }
+
   viz_funnel_2 <- function(x){metaviz::viz_funnel(x, y_axis = "precision")}
 
   poss_viz_funnel <- possibly(viz_funnel_2, otherwise = NA)
@@ -114,8 +128,14 @@ make_viz <- function(data) {
                               .else = ~return(NA)
                             ),
                             NA
+      ),
+      model_params = map_if(
+        .x = model,
+        .p = possibly(\(x) class(x) %in% parameters::supported_models() %>% any(), otherwise = NA), 
+        .f = possibly(parameters::parameters, NA),
+        .else = ~ return(NA),
       )
     )
 
   return(viz_out)
-}
+}