feat!: update meta-analysis function arguments cf. changes in bdf5375 #…

…118

- `fit_MA_mv()` switch from `Z_colname` to `outcome_colname` and `VZ_colname` to `outcome_SE_colname` #118
- `calculate_deviation_score()` #118:  expose colname selection column to user , return meta-analytic mean in column, add `outcome_colname` to specify relevant outcome column name, add argument checking #116
- `apply_slice_conditionally()` create column `exclusion_set` if does not exist and switch to `se_log` as `VZ_colname` in example for eucalyptus` #188
- `box_cox_transform()` expose `outcome_SE_colname` arg, add relevant arg checks, update cli output, parse outcome_SE_columns` to `folded_params()` function call #188
- replace `std.error_log` arg with `se_log` for any logged `yi` data #188
- and update wrapper function calls in `meta_analyse_datasets()`, provide NULL default for new arg `outcome_SE` #188
- `fit_box_cox_ratings_cat()`: add arg checks #116, replace arg name `.data` with `data`, switch to `inject` approach as per `fit_box_cox_ratings_cont()` #97
- fit_box_cox_ratings_cont()` switch to `data` arg from `.data`, add arg checks #116
- add arg checks to `standardise_response()`, `fit_sorensen_glm()`, `fit_multivar_MA()` , `fit_uni_mixed_effects()` #116

docs!: update corresponding docs

- delete obsolete `@examples` code #102
- add #TODO comments #97
- add missing `@import` calls
- switch to `@importFrom` from `@import` for `purrr::` `rlang::` and `tidyr::` et al
- don't import `rlang::` in package.R

bug!: minor bug fixes

- unquote bare colname in pointblank check for `exclude_extreme_VZ()`
- `fit_metafor_mv_reduced()` replace obs_id with study_id as random effect

R/ManyEcoEvo-package.R

@@ -2,7 +2,6 @@
 "_PACKAGE"
 
 ## usethis namespace: start
-#' @import rlang
 #' @importFrom glue glue
 #' @importFrom lifecycle deprecated
 ## usethis namespace: end

R/anonymise_teams.R

@@ -6,7 +6,7 @@
 #' @return A `df` with anonymised values of `id_col` based on the `New_Identifier` colum of `lookup`
 #' @export
 #' @importFrom pointblank col_vals_not_null
-#' @import tidyr
+#' @importFrom tidyr separate unite
 #' @import dplyr
 anonymise_teams <- function(df, lookup) { # TODO actually... this is anonymise_id_col()
   df %>%

R/apply_VZ_exclusions.R

@@ -27,7 +27,7 @@
 #'     dataset_log_transform = "eucalyptus") %>%
 #'   generate_yi_subsets() %>% 
 #'   apply_VZ_exclusions(VZ_colname = 
-#'                         list("eucalyptus" = "std.error_log", 
+#'                         list("eucalyptus" = "se_log", 
 #'                              "blue tit" = "VZ"), 
 #'                       VZ_cutoff = 3)
 apply_VZ_exclusions <- function(df = data.frame(), VZ_colname, VZ_cutoff) {

R/assign_transformation_type.R

@@ -8,7 +8,7 @@
 #' Based on the response transformation and link function, the function assigns the back-transformation type to be applied to the analysts' point-estimates. The function and assigns the identity transformation if the effects were reported on the link-scale and the estimates are already back-transformed the original response variable prior to modelling. When either of these cases is not true for a given analysis, the function returns the value of the `link_fun` or `response_transformation` argument. When an analysis has been reported on the link-scale and the analyst transformed the response variable prior to modelling, the function assigns the `"double-transformation"`  value for that analysis. When the `response_transformation` and `link_fun` arguments are missing, the function assigns the `"identity"` value to the analysis, assuming that `NA` values are equivalent to the identity transformation. 
 #' @export
 #' @import dplyr
-#' @import rlang
+#' @importFrom rlang is_na
 #' @family Back-transformation
 #' @seealso [prepare_response_variables_yi(), standardise_response()]. To be called prior to [clean_response_transformation()].
 assign_transformation_type <- function(response_transformation = character(1L),

R/box_cox_transform.R

@@ -1,7 +1,7 @@
 #' Box-cox transform absolute deviation from the meta-analytic mean scores
 #'
 #' @param data Dataset for model fitting, must contain columns `"abs_deviation_score_estimate"` and standard error
-#' @param dataset character string of either "blue tit" or "eucalyptus"
+#' @param dataset character string of length 1 either "blue tit" or "eucalyptus"
 #' @family Box-Cox transformation
 #' @family Analysis-level functions
 #' @return data with additional columns of box-cox transformed deviation scores and variance
@@ -16,62 +16,73 @@
 #' @importFrom stringr str_starts
 #' @importFrom tidyr hoist
 #' @seealso [variance_boxcox()], [folded_params()]
-box_cox_transform <- function(data, dataset) {
+box_cox_transform <- function(data, dataset, outcome_SE_colname) {
+  # ----- Argument Checking -----
   if ( any(rlang::is_na(data), rlang::is_null(data))) {
-    cli::cli_alert_warning(text = 
-                             glue::glue(
-                               "Cannot box-cox transform data for",
-                               paste(names(dplyr::cur_group()),
-                                     dplyr::cur_group(),
-                                     sep = " = ",
-                                     collapse = ", "
-                               )
-                             ))
+    cli::cli_alert_warning(
+      text = 
+        glue::glue(
+          "Cannot box-cox transform data for",
+          paste(names(dplyr::cur_group()),
+                dplyr::cur_group(),
+                sep = " = ",
+                collapse = ", "
+          )
+        ))
     result <- NA
   } else {
+
+    # ---- Compute Box-Cox Transformed absolute deviation scores ----
     cli::cli_h2(
-      c("Box-cox transforming absolute deviation scores for {.arg dataset} = {.val {dataset}}.")
+      c("Box-cox transforming absolute deviation scores for ",
+        "{.arg dataset} = {.val {dataset}}."
+      )
     )
 
-    box_cox_recipe <- recipes::recipe(~.,
-                                      data = select(
-                                        data,
-                                        starts_with("abs_deviation_score_")
-                                      )
-    ) %>%
+    box_cox_recipe <- 
+      recipes::recipe( ~ .,
+                       data = select(
+                         data,
+                         starts_with("abs_deviation_score_")
+                       )
+      ) %>%
       timetk::step_box_cox(everything(), limits = c(0, 10)) %>%
       recipes::prep(training = data, retain = TRUE) # estimate lambda + box cox transform vars
 
     if (box_cox_recipe %>%
-        recipes::tidy(number = 1) %>% nrow() > 0) { # TODO pull execution of if/else and check result in if() so not executing twice (next line)
+        recipes::tidy(number = 1) %>% 
+        nrow() > 0) { # TODO pull execution of if/else and check result in if() so not executing twice (next line)
       lambda <- box_cox_recipe %>%
         recipes::tidy(number = 1) %>%
         pull(., lambda) %>%
         `names<-`(., pull(box_cox_recipe %>%
-                            recipes::tidy(number = 1), terms))
+                            recipes::tidy(number = 1), 
+                          terms))
 
       if (!is.null(dataset)) {
         cli::cli_alert_info(
           c(
             "Optimised Lambda used in Box-Cox Transformation of ",
             "{dataset} dataset variables ",
-            "is {round(lambda, 4)} for `{names(lambda)}`."
+            "is {round(lambda, 2)} for `{names(lambda)}`."
           )
         )
       }
 
-      # Z_colname <- data %>% colnames %>% keep(., str_starts(., "Z"))
-      VZ_colname <- data %>%
-        colnames() %>%
-        keep(., str_starts(., "VZ"))
-
       result <- recipes::juice(box_cox_recipe) %>%
         rename_with(.fn = ~ paste0("box_cox_", .x)) %>%
         bind_cols(data, .) %>%
-        mutate(fold_params = map2(.x = abs_deviation_score_estimate, .y = !!as.name(VZ_colname), .f = folded_params)) %>%
-        hoist(fold_params, folded_mu_val = 1, folded_v_val = 2) %>%
+        mutate(fold_params = 
+                 map2(.x = abs_deviation_score_estimate, 
+                      .y = !!as.name(outcome_SE_colname), 
+                      .f = folded_params)) %>%
+        hoist(fold_params, 
+              folded_mu_val = 1, 
+              folded_v_val = 2) %>%
         mutate(
-          box_cox_var = variance_box_cox(folded_mu_val, folded_v_val, lambda[[1]]),
+          box_cox_var = variance_box_cox(folded_mu_val, 
+                                         folded_v_val, 
+                                         lambda[[1]]),
           lambda = lambda[[1]]
         )
     } else {
@@ -99,15 +110,17 @@ variance_box_cox <- function(folded_mu, folded_v, lambda) {
 #' Calculate the folded parameters for the Box-Cox transformation
 #' @param abs_dev_score The absolute deviation score
 #' @param VZr The variance of the standardised effect size
-#' @return A list containing the mean and variance of the folded parameters
+#' @return A named list containing the mean `fold_mu` and variance `fold_v` of the folded parameters
 #' @export
 #' @family Box-Cox transformation
 #' @family Analysis-level functions
 folded_params <- function(abs_dev_score, VZr) {
+
   mu <- abs_dev_score
   sigma <- sqrt(VZr)
   fold_mu <- sigma * sqrt(2 / pi) * exp((-mu^2) / (2 * sigma^2)) + mu * (1 - 2 * pnorm(-mu / sigma)) # folded abs_dev_score
   fold_se <- sqrt(mu^2 + sigma^2 - fold_mu^2)
   fold_v <- fold_se^2 # folded VZr
   return(list(fold_mu = fold_mu, fold_v = fold_v))
+
 }

R/calculate_deviation_score.R

@@ -3,38 +3,52 @@
 #'
 #' @param data Dataset containing the column `"Z"` \(standardised out of sample prediction piont estimates\) or standardised effect-sizes, `"Zr"`.
 #' @param meta_analytic_model Fitted meta-analytic model of class `rma`
-#'
+#' @param outcome_colname Column name in `data` containing outcomes used in meta-analysis
 #' @return `data` but with an additional column of deviation score estimates `"abs_deviation_score_estimate"`
 #' @export
 #'
-#' @examples
-#' # Example Usage: TODO - this is incorrect bc box-cox transformation is included in meta_analyse_datasets_wrapper
-#' # library(tidyverse); library(targets)
-#' # tar_load(meta_analysis_outputs)
-#' # calculate_deviation_score(meta_analysis_outputs$data[[1]],
-#' #                           meta_analysis_outputs$MA_mod[[1]])
 #' @import dplyr
 #' @importFrom cli cli_h2
-calculate_deviation_score <- function(data, meta_analytic_model) {
-  cli::cli_h2(c("Calculating absolute deviation scores from standardised effect sizes"))
-
-  stopifnot("rma" %in% class(meta_analytic_model))
-  # TODO build in checks:
-  # stop if not data is not a dataframe
-  # stop if not data does not contain the column Zr/Z
-
-
-  meta_analytic_mean <- meta_analytic_model$beta[[1]] # TODO should this be the fitted val from the model, or the raw score?
-  # meta_analytic_se <- meta_analytic_model$se[[1]]
-
+calculate_deviation_score <- function(data, meta_analytic_model, outcome_colname) {
+
+  # ----- Argument Checks -----
+  stopifnot(
+    "rma" %in% class(meta_analytic_model),
+    is.data.frame(data)
+  )
+
+  pointblank::test_col_exists(data, outcome_colname)
+
+  meta_analytic_mean <- meta_analytic_model$beta[[1]]
+
+  if (any(rlang::is_null(meta_analytic_mean), 
+          rlang::is_na(meta_analytic_mean))) {
+    cli::cli_abort(
+      text = 
+        c("{.arg meta_analytic_mean} could not be calculated ",
+          "({.val {meta_analytic_mean}}): ",
+          "check the meta-analytic model in {.arg meta_analytic_model}."
+        )
+    )
+  }
+
+  # ----- Calculate deviation scores -----
+  cli::cli_h2(c("Calculating absolute deviation scores from meta-analytic mean"))
+
+  cli::cli_alert_info(
+    text = 
+      c("Using the meta-analytic mean outcome as the reference point:",
+        "{.arg meta_analytic_mean} = {.val {round(meta_analytic_mean, 4)}}")
+  )
+
   out <- data %>%
     ungroup() %>%
     mutate(across(
-      .cols = starts_with("Z"),
-      .fns = ~ abs(.x - {
-        meta_analytic_mean
-      }),
+      .cols = {{outcome_colname}},
+      .fns = ~ abs(.x - meta_analytic_mean),
       .names = "abs_deviation_score_estimate"
-    )) # abs_deviation_score_se = abs({VZr - meta_analytic_se})?
+    ), 
+    meta_analytic_mean = meta_analytic_mean) 
+
   return(out)
 }

R/compute_MA_inputs.R

@@ -24,6 +24,7 @@
 #' @importFrom purrr map2 map pmap
 #' @importFrom rlang is_null
 #' @importFrom pointblank col_exists test_col_exists
+#' @importFrom cli cli_alert_warning cli_h1
 compute_MA_inputs <- function(ManyEcoEvo, estimate_type = NULL) {
   # TODO should be renamed something to do with diversity indices... that's the
   # only thing happening here!!

R/exclude_extreme_VZ.R

@@ -7,12 +7,13 @@
 #' @export
 #' @import dplyr
 #' @importFrom pointblank col_exists
-#' @importFrom cli cli_alert
+#' @importFrom cli cli_alert cli_h2
 #' @importFrom tidyselect all_of
+#' @importFrom purrr pluck
 exclude_extreme_VZ <- function(df = data.frame(), VZ_colname, VZ_cutoff = numeric(1L)) {
 
   # ----- Argument Checking -----
-  pointblank::col_exists(df, VZ_colname)
+  pointblank::col_exists(df, !!{{VZ_colname}})
   stopifnot(
     length(VZ_cutoff) == 1, 
     is.numeric(VZ_cutoff),

R/filt_multivar_MA.R

@@ -3,7 +3,7 @@
 #' @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 ... Additional arguments passed to [lme4::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.
 #'
@@ -32,69 +32,89 @@
 #' - `ReviewerId`: reviewer identifier
 #' @family Model fitting and meta-analysis
 fit_multivar_MA <- function(data_tbl, N = 5, ..., env = rlang::caller_env()) {
+
+  # ----- Argument Checks -----
+
   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),
+    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)),
+                  PublishableAsIs +
+                  mean_diversity_index +
+                  (1 | ReviewerId)),
     env = env
   )
-
-  f2 <- rlang::new_formula(rlang::expr(box_cox_abs_deviation_score_estimate),
+
+  f2 <- rlang::new_formula(
+    rlang::expr(box_cox_abs_deviation_score_estimate),
     rlang::expr(RateAnalysis +
-      PublishableAsIs +
-      mean_diversity_index +
-      mixed_model +
-      (1 | ReviewerId)),
+                  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 = ": "
+          sep = ": "
     ) %>%
     pull(group)
-
+
+  # ---- Conditionally Fit Models -----
+
   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)))))
+
+    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)))))
+
+    cli::cli_bullets(c(
+      setNames(cur_group_bullets, 
+               rep("*", length(cur_group_bullets)))
+    ))
   }
-
-  # TODO MAKE SURE GIVES CORRECT EX
+  
+  # TODO MAKE SURE GIVES CORRECT EXPECTED OUTPUT
   f <- if (pass_threshold) f2 else f1 # MAKE SURE RETURNS APPROPIRATELY
-
+  
   mod <- rlang::inject(lme4::lmer(!!f, data = data_tbl, ...))
-
+  
   return(mod)
 }