diff --git a/R/ManyEcoEvo-package.R b/R/ManyEcoEvo-package.R index 97c07e1..78a8bfd 100644 --- a/R/ManyEcoEvo-package.R +++ b/R/ManyEcoEvo-package.R @@ -2,7 +2,6 @@ "_PACKAGE" ## usethis namespace: start -#' @import rlang #' @importFrom glue glue #' @importFrom lifecycle deprecated ## usethis namespace: end diff --git a/R/anonymise_teams.R b/R/anonymise_teams.R index 96302f6..0936cf1 100644 --- a/R/anonymise_teams.R +++ b/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 %>% diff --git a/R/apply_VZ_exclusions.R b/R/apply_VZ_exclusions.R index a17322a..4dc627d 100644 --- a/R/apply_VZ_exclusions.R +++ b/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) { diff --git a/R/assign_transformation_type.R b/R/assign_transformation_type.R index 4ed5c1f..5db3799 100644 --- a/R/assign_transformation_type.R +++ b/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), diff --git a/R/box_cox_transform.R b/R/box_cox_transform.R index 6cc846f..a52f4cb 100644 --- a/R/box_cox_transform.R +++ b/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)) + } diff --git a/R/calculate_deviation_score.R b/R/calculate_deviation_score.R index c997e2b..bf4969d 100644 --- a/R/calculate_deviation_score.R +++ b/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) } diff --git a/R/compute_MA_inputs.R b/R/compute_MA_inputs.R index a8b2397..5805a2f 100644 --- a/R/compute_MA_inputs.R +++ b/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!! diff --git a/R/exclude_extreme_VZ.R b/R/exclude_extreme_VZ.R index c413f5d..04ffb07 100644 --- a/R/exclude_extreme_VZ.R +++ b/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), diff --git a/R/filt_multivar_MA.R b/R/filt_multivar_MA.R index fe4d290..89c1480 100644 --- a/R/filt_multivar_MA.R +++ b/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) } diff --git a/R/fit_MA_mv.R b/R/fit_MA_mv.R index 67d952d..d335168 100644 --- a/R/fit_MA_mv.R +++ b/R/fit_MA_mv.R @@ -1,9 +1,10 @@ #' Fit Meta-regression with random-effects -#' @description Fit a meta-regression model with random effects using the [metafor::rma.mv()] function from the `metafor` package to a dataframe containing the estimates and variances for the meta-analysis. +#' +#' @description Fit a meta-regression model with random effects using the [metafor::rma.mv()] function from the [metafor::] package to a data.frame containing the estimates and variances for the meta-analysis. #' #' @param effects_analysis A dataframe containing the estimates and variances for the meta-analysis. -#' @param Z_colname The name of the column containing the estimates. -#' @param VZ_colname The name of the column containing the variances. +#' @param outcome_colname The name of the column containing the estimates. +#' @param outcome_SE_colname The name of the column containing the variances. #' @param estimate_type The type of estimate to be used in the model. One of `c("Zr", "y50", "y25", "y75", or "yi")`. #' #' @return A fitted model of class `c("rma.mv","rma")`. @@ -21,14 +22,19 @@ #' fit_MA_mv(beta_estimate, beta_SE, "Zr") #' @seealso [fit_metafor_mv()], [metafor::rma.mv()] #' @family Model fitting and meta-analysis -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")) +fit_MA_mv <- function(effects_analysis = data.frame(), + outcome_colname, + outcome_SE_colname, + estimate_type = character(1L)) { + + pointblank::stop_if_not(estimate_type %in% + c("Zr", "yi", "y25", "y50", "y75")) - Z <- effects_analysis %>% dplyr::pull({{ Z_colname }}) - VZ <- effects_analysis %>% dplyr::pull({{ VZ_colname }}) + outcome <- effects_analysis %>% dplyr::pull({{ outcome_colname }}) + outcome_variance <- effects_analysis %>% dplyr::pull({{ outcome_SE_colname }}) mod <- ManyEcoEvo::fit_metafor_mv( - estimate = Z, - variance = VZ, + estimate = outcome, + variance = outcome_variance, estimate_type = estimate_type, data = effects_analysis ) diff --git a/R/fit_boxcox_ratings_cat.R b/R/fit_boxcox_ratings_cat.R index 7dd73b7..d31e6db 100644 --- a/R/fit_boxcox_ratings_cat.R +++ b/R/fit_boxcox_ratings_cat.R @@ -1,7 +1,7 @@ #' Fit model of boxcox deviation scores as function of continuous ratings #' @description Fit an lmer model of the box-cox transformed deviation from the meta-analytic mean scores as a function of continuous peer-review ratings #' -#' @param .data Data for model fitting +#' @param data Data for model fitting #' @param outcome outcome variable, unquoted. #' @param outcome_var Variance of the `outcome` variable #' @param interceptless A logical relating to whether the model should be interceptless or not. Defaults to `FALSE`. @@ -13,63 +13,87 @@ #' @importFrom rlang ensym new_formula expr #' @import dplyr #' @importFrom forcats fct_relevel -fit_boxcox_ratings_cat <- function(.data, outcome, outcome_var, interceptless = FALSE) { - cli::cli_h2(c("Fitting lmer with categorical ratings predictor on box_cox_transformed outcomes")) - # Example Usage: - # library(tidyverse);library(targets);library(metafor) - # tar_load(meta_analysis_outputs) - # meta_analysis_outputs$data[[1]] %>% - # fit_boxcox_ratings_cat(., - # outcome = box_cox_abs_deviation_score_estimate, - # outcome_var = VZr, interceptless = FALSE) - - # TODO @egouldo stopifnot data doesn't contain variables named eval(box_cox_outcome_var), eval(sampling_variance_var), review_data - # TODO @egouldo unnest and then check stopifnot: RateAnalysis, ReviewerId, study_id. +#' @importFrom tidyr unnest +#' @importFrom cli cli_h2 +#' @importFrom pointblank expect_col_exists +#' @examples +#' # Example Usage: +#' # library(tidyverse);library(targets);library(metafor) +#' # tar_load(meta_analysis_outputs) +#' # meta_analysis_outputs$data[[1]] %>% +#' # fit_boxcox_ratings_cat(., +#' # outcome = box_cox_abs_deviation_score_estimate, +#' # outcome_var = VZr, interceptless = FALSE) +fit_boxcox_ratings_cat <- function(data, outcome, outcome_var, interceptless = FALSE, ..., env = rlang::caller_env()) { + + # ----- Argument Checks ----- + stopifnot( + is.data.frame(data), + is.logical(interceptless) + ) + + pointblank::expect_col_exists( + data, + columns = c(starts_with("box_cox_abs_"), + {{outcome}}, + {{outcome_var}}, + study_id, + review_data) + ) + + # ----- Fit model ----- + cli::cli_h2(c("Fitting {.fn lmer} with categorical ratings predictor {.arg PublishableAsIs} on Box-Cox transformed outcome: outcome: {.val {outcome}}")) + data_tbl <- - .data %>% + data %>% unnest(cols = c(review_data)) %>% + pointblank::col_exists( + columns = c("PublishableAsIs", "ReviewerId") + ) %>% select( study_id, ReviewerId, PublishableAsIs, - starts_with("box_cox_"), # @TODO - delete this row? + starts_with("box_cox_"), #TODO - delete this row? {{ outcome }}, {{ outcome_var }} ) %>% ungroup() %>% mutate( - PublishableAsIs = forcats::fct_relevel(PublishableAsIs, c( - "deeply flawed and unpublishable", - "publishable with major revision", - "publishable with minor revision", - "publishable as is" - )), + PublishableAsIs = + forcats::fct_relevel( + PublishableAsIs, + c("deeply flawed and unpublishable", + "publishable with major revision", + "publishable with minor revision", + "publishable as is") + ), obs_id = 1:n() ) - + if (interceptless == FALSE) { f <- rlang::new_formula( rlang::ensym(outcome), - expr( + rlang::expr( PublishableAsIs + - (1 | ReviewerId) # + (1 | study_id ) RE ommitted due to convergence issues - ) + (1 | ReviewerId) + ), + env = env ) - - mod <- lme4::lmer(formula = f, data = data_tbl) + } else { - ( # interceptless: for plotting - - mod <- lme4::lmer( #TODO implement same inject approach as fit_boxcox_ratings_cont - rlang::new_formula( - rlang::ensym(outcome), - expr(-1 + PublishableAsIs + (1 | ReviewerId)) - ), - data = data_tbl - ) + # interceptless: for plotting + f <- rlang::new_formula( + rlang::ensym(outcome), + rlang::expr( + -1 + PublishableAsIs + (1 | ReviewerId) + ), + env = env ) } - + + mod <- rlang::inject(lme4::lmer(!!f, data = data_tbl, ...)) + return(mod) } @@ -79,6 +103,6 @@ fit_boxcox_ratings_cat <- function(.data, outcome, outcome_var, interceptless = #' @family Model fitting and meta-analysis #' @importFrom purrr possibly poss_fit_boxcox_ratings_cat <- purrr::possibly(fit_boxcox_ratings_cat, - otherwise = NA, - quiet = FALSE + otherwise = NA, + quiet = FALSE ) diff --git a/R/fit_boxcox_ratings_cont.R b/R/fit_boxcox_ratings_cont.R index 4494ddb..9529b44 100644 --- a/R/fit_boxcox_ratings_cont.R +++ b/R/fit_boxcox_ratings_cont.R @@ -1,36 +1,37 @@ #' Fit model of boxcox deviation scores as function of continuous ratings #' @description Fit an lmer model of the box-cox transformed deviation from the meta-analytic mean scores as a function of continuous peer-review ratings #' -#' @param .data Data for model fitting +#' @param data Data for model fitting #' @param outcome outcome variable, unquoted. #' @param outcome_var Variance of the `outcome` variable #' #' @return An object of class `lme4::lmerMod-class` #' @export #' @family Model fitting and meta-analysis -#' @examples -#' # Example Usage: -#' # Note, outcome = the name using NSE of the response variable, otucome_var = the -#' # variance associated with that variable. -#' # library(tidyverse);library(targets);library(metafor);library(tidymodels);library(multilevelmod) -#' # tar_load(meta_analysis_outputs) -#' # meta_analysis_outputs$data[[1]] %>% -#' # fit_boxcox_ratings_cont(., -#' # box_cox_abs_deviation_score_estimate, -#' # VZr ) #' @import dplyr #' @importFrom lme4 lmer #' @importFrom rlang ensym new_formula inject expr caller_env #' @importFrom tidyr unnest #' @importFrom cli cli_h2 -fit_boxcox_ratings_cont <- function(.data, outcome, outcome_var, ..., env = rlang::caller_env()) { - cli::cli_h2(c("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 - # TODO @egouldo unnest and then check stopifnot: RateAnalysis, ReviewerId, study_id. - +fit_boxcox_ratings_cont <- function(data, outcome, outcome_var, ..., env = rlang::caller_env()) { + + # ----- Argument Checks ----- + stopifnot( + is.data.frame(data) + ) + pointblank::expect_col_exists( + data, + columns = c(starts_with("box_cox_abs_"), + {{outcome}}, + {{outcome_var}}, + study_id, + review_data) + ) + + # ----- Fit model ----- + cli::cli_h2(c("Fitting {.fn lmer} with continuous ratings predictor {.arg RateAnalysis} on Box-Cox transformed outcome: {.val {outcome}}")) data_tbl <- # TODO, consider extracting unnesting outside of this fn. - .data %>% + data %>% select( study_id, review_data, @@ -40,18 +41,20 @@ fit_boxcox_ratings_cont <- function(.data, outcome, outcome_var, ..., env = rlan ) %>% unnest(cols = c(review_data)) %>% ungroup() %>% - mutate(., obs_id = 1:nrow(.)) #TODO Is this correct, or should it be done BEFORE the unnesting?? - - f <- rlang::new_formula(rlang::ensym(outcome), - expr( + pointblank::col_exists(columns = c("RateAnalysis", "ReviewerId")) %>% + mutate(., obs_id = 1:nrow(.)) + + f <- rlang::new_formula( + rlang::ensym(outcome), + rlang::expr( RateAnalysis + (1 | study_id) # NOTE: ReviewerId removed due to singularity ), env = env ) - + mod <- rlang::inject(lme4::lmer(!!f, data = data_tbl, ...)) - + return(mod) } @@ -61,6 +64,6 @@ fit_boxcox_ratings_cont <- function(.data, outcome, outcome_var, ..., env = rlan #' @family Model fitting and meta-analysis #' @importFrom purrr possibly poss_fit_boxcox_ratings_cont <- purrr::possibly(fit_boxcox_ratings_cont, - otherwise = NA, - quiet = FALSE + otherwise = NA, + quiet = FALSE ) diff --git a/R/fit_metafor_mv.R b/R/fit_metafor_mv.R index 576626c..a95dce6 100644 --- a/R/fit_metafor_mv.R +++ b/R/fit_metafor_mv.R @@ -12,29 +12,19 @@ #' @importFrom metafor rma.mv #' @export #' -#' @examples -#' # TODO -- is this the best way of setting up this fun?? (i.e. to take numeric vectors)? -#' # Example Usage: -#' # library(tidyverse);library(targets);library(metafor) # NOT RUN, TODO: remove after create pkg -#' # source("R/functions.R") #NOT RUN, TODO: remove after create pkg -#' tar_read(round_2_survey_meta_analysis) %>% -#' filter(dataset == "eucalyptus") %>% -#' filter( -#' !is.na(Zr), -#' !is.na(VZr), -#' !is.infinite(Zr), -#' !is.infinite(VZr) -#' ) %>% -#' fit_metafor_mv(estimate = .$Zr, variance = .$VZr, estimate_type = "Zr", data = .) #' @family Model fitting and meta-analysis fit_metafor_mv <- function(estimate, variance, estimate_type = character(1L), data) { + cli::cli_h2(c("Fitting metaregression")) - match.arg(estimate_type, choices = c("Zr", "y50", "y25", "y75", "yi"), several.ok = FALSE) - + + match.arg(estimate_type, + choices = c("Zr", "y50", "y25", "y75", "yi"), + several.ok = FALSE) + if (estimate_type != "Zr") { # you need to put SE^2 for y variance <- variance^2 } - + mod <- metafor::rma.mv( yi = estimate, # of type "Zr" or "ymed", "y25", "y75" V = variance, @@ -54,46 +44,32 @@ fit_metafor_mv <- function(estimate, variance, estimate_type = character(1L), da #' @param estimate Numeric vector #' @param variance Numeric vector #' @param estimate_type Character vector of either "Zr", "y50", "y25", "y75". -#' @param data Dataframe containing estimates and variances with case id column \code{study_id}. +#' @param data Dataframe containing estimates and variances with case id column `study_id`. #' -#' @return Object of class \code{rma.mv} -#' @import metafor +#' @return Object of class `rma.mv` +#' @importFrom metafor rma.mv #' @import dplyr -#' @importFrom glue glue #' @importFrom cli cli_h2 #' @export -#' -#' @examples -#' # TODO -- is this the best way of setting up this fun?? (i.e. to take numeric vectors)? -#' # Example Usage: -#' # library(tidyverse);library(targets);library(metafor) # NOT RUN, TODO: remove after create pkg -#' # source("R/functions.R") #NOT RUN, TODO: remove after create pkg -#' tar_read(round_2_survey_meta_analysis) %>% -#' filter(dataset == "eucalyptus") %>% -#' filter( -#' !is.na(Zr), -#' !is.na(VZr), -#' !is.infinite(Zr), -#' !is.infinite(VZr) -#' ) %>% -#' fit_metafor_mv(estimate = .$Zr, variance = .$VZr, estimate_type = "Zr", data = .) #' @family Model fitting and meta-analysis fit_metafor_mv_reduced <- function(estimate, variance, estimate_type = character(1L), data) { + cli::cli_h2(c("Fitting multivariate metaregression")) + match.arg(estimate_type, choices = c("Zr", "y50", "y25", "y75"), several.ok = FALSE) - + data <- data %>% ungroup() %>% mutate(obs_id = 1:n()) # metafor idiosyncrasy - + if (estimate_type != "Zr") { # you need to put SE^2 for y variance <- variance^2 } - + metafor::rma.mv( yi = estimate, # of type "Zr" or "ymed", "y25", "y75" V = variance, - random = list(~ 1 | TeamIdentifier / obs_id), + random = list(~ 1 | TeamIdentifier / study_id), data = data, sparse = TRUE, sigma2 = c(0, NA) @@ -102,10 +78,11 @@ fit_metafor_mv_reduced <- function(estimate, variance, estimate_type = character #' Possibly [fit_metafor_mv()] #' @description Wrapper for [fit_metafor_mv()] that returns `NA` if an error is thrown +#' @importFrom purrr possibly #' @keywords internal #' @family Model fitting and meta-analysis #' @seealso [fit_metafor_mv()] poss_fit_metafor_mv <- purrr::possibly(fit_metafor_mv, - otherwise = NA, - quiet = FALSE + otherwise = NA, + quiet = FALSE ) diff --git a/R/fit_sorensen_glm.R b/R/fit_sorensen_glm.R index aa75429..aaa7757 100644 --- a/R/fit_sorensen_glm.R +++ b/R/fit_sorensen_glm.R @@ -1,39 +1,44 @@ #' Fit univariate glm of deviation scores on sorensen diversity index #' -#' @param data A dataframe containing box-cox transformed absolute deviation from the meta-analytic mean scores +#' @param data A `data.frame` containing Box-Cox transformed absolute deviation from the meta-analytic mean scores and the mean Sorensen's scores for each analysis #' -#' @return An fitted model object of class `glm` and `parsnip` +#' @return A fitted model object of class `glm` and `parsnip` #' @export #' @family Model fitting and meta-analysis -#' @examples -#' # library(tidyverse);library(targets);library(metafor);library(tidymodels) -#' # tar_load(meta_analysis_outputs) -#' # fit_sorensen_glm(meta_analysis_results$data[[1]]) +#' @importFrom parsnip linear_reg fit +#' @importFrom recipes recipe update_role +#' @importFrom workflows workflow add_model add_recipe extract_fit_parsnip +#' @import dplyr +#' @importFrom cli cli_h2 cli_alert_info +#' @importFrom pointblank expect_col_exists +#' @importFrom purrr simplify fit_sorensen_glm <- function(data) { cli::cli_h2(c("Fitting glm for box-cox transformed outcome with sorensen diversity index as predictor")) cli::cli_alert_info(dplyr::cur_group() %>% purrr::simplify()) # TODO only run when applied within a tibble on a list-col... want fn available on its own - - - # TODO: stopifnot(does not contain box_cox_abs_ or mean_diversity_index) - + + pointblank::expect_col_exists(data, + columns = c(starts_with("box_cox_abs_"), + "mean_diversity_index")) + data <- data %>% - dplyr::select(dplyr::starts_with("box_cox_abs_"), mean_diversity_index) - + dplyr::select(dplyr::starts_with("box_cox_abs_"), + mean_diversity_index) + glm_recipe <- recipes::recipe(~., - data = data + data = data ) %>% recipes::update_role(starts_with("box_cox_abs_"), new_role = "outcome") - + glm_mod <- parsnip::linear_reg(engine = "glm") - + fitted_mod <- workflows::workflow() %>% workflows::add_model(glm_mod) %>% workflows::add_recipe(glm_recipe) %>% parsnip::fit(data = data) %>% workflows::extract_fit_parsnip() - + return(fitted_mod) } @@ -42,6 +47,6 @@ fit_sorensen_glm <- function(data) { #' @keywords internal #' @importFrom purrr possibly poss_fit_sorensen_glm <- purrr::possibly(fit_sorensen_glm, - otherwise = NA, - quiet = FALSE -) # TODO export + otherwise = NA, + quiet = FALSE +) diff --git a/R/fit_uni_mixed_effects.R b/R/fit_uni_mixed_effects.R index dca16c4..6c0a576 100644 --- a/R/fit_uni_mixed_effects.R +++ b/R/fit_uni_mixed_effects.R @@ -1,4 +1,4 @@ -#' Fit univariate glm of deviation scores on random effects inclusion +#' Fit model of Box-Cox transformed deviation scores as a function random-effects inclusion in analyses #' @description Fits a univariate glm of box-cox transformed absolute deviation from the meta-analytic mean scores as a function of whether the analysis was a mixed effects model \(i.e. included random effects\) or not. #' #' @param data Dataframe containing box-cox transformed absolute deviation scores and binary column called `mixed_model` describing whether or not the analysis used a mixed-effects model. @@ -20,51 +20,56 @@ #' @importFrom workflows workflow add_model add_recipe extract_fit_parsnip #' @seealso [parsnip::details_linear_reg_glm] for details on the [parsnip::linear_reg] engine. fit_uni_mixed_effects <- function(data) { + cli::cli_h2(c("Fitting glm for box-cox transformed outcome with inclusion of random effects (binary variable) as predictor")) - if (pointblank::test_col_exists( - data, - c( - "mixed_model", - data %>% - select(starts_with("box_cox_abs_")) %>% #TODO will this not fail if missing?? Seems illogical to derive the col name we are testing to look for from the object?? - colnames() + if (pointblank::test_col_exists(data, + columns = c("mixed_model", + starts_with("box_cox_abs_")))) { + + cli::cli_alert_warning( + c("Columns {.var mixed_model} and ", + data %>% + dplyr::select(starts_with("box_cox_abs_")) %>% + colnames(), + " missing. Returning {.val {NA}}.") ) - )) { - if (length(unique(data$mixed_model)) == 1) { - cli::cli_warn(message = "More than 1 unique value of {.var mixed_model} is needed to fit model with {.var mixed_model} as predictor variable.") - - return(NA) - } else { - data <- data %>% - dplyr::select(dplyr::starts_with("box_cox_abs_"), mixed_model) - - glm_recipe <- - recipes::recipe(~., - data = data - ) %>% - recipes::update_role(starts_with("box_cox_abs_"), new_role = "outcome") %>% - recipes::step_mutate(mixed_model = as.factor(mixed_model)) %>% - recipes::step_naomit() - - glm_mod <- parsnip::linear_reg(engine = "glm") - - fitted_mod <- - workflows::workflow() %>% - workflows::add_model(glm_mod) %>% - workflows::add_recipe(glm_recipe) %>% - parsnip::fit(data = data) %>% - workflows::extract_fit_parsnip() - - return(fitted_mod) - } - } else { - cli::cli_alert_warning(c("Columns {.var mixed_model} and ", data %>% - dplyr::select(starts_with("box_cox_abs_")) %>% - colnames(), " missing. Returning {.val {NA}}.")) return(NA) - } + + } else if ( length(unique(data$mixed_model)) == 1) { + + cli::cli_warn(message = "More than 1 unique value of {.var mixed_model} ", + "is needed to fit model with {.var mixed_model} ", + "as predictor variable. Returning {.val {NA}}") + + return(NA) + + } else { + + data <- data %>% + dplyr::select(dplyr::starts_with("box_cox_abs_"), + mixed_model) + + glm_recipe <- + recipes::recipe(~., + data = data + ) %>% + recipes::update_role(starts_with("box_cox_abs_"), new_role = "outcome") %>% + recipes::step_mutate(mixed_model = as.factor(mixed_model)) %>% + recipes::step_naomit() + + glm_mod <- parsnip::linear_reg(engine = "glm") + + fitted_mod <- + workflows::workflow() %>% + workflows::add_model(glm_mod) %>% + workflows::add_recipe(glm_recipe) %>% + parsnip::fit(data = data) %>% + workflows::extract_fit_parsnip() + } + + return(fitted_mod) } #' Possibly [fit_uni_mixed_effects()] diff --git a/R/generate_outlier_subsets.R b/R/generate_outlier_subsets.R index 2a06025..3bf73c5 100644 --- a/R/generate_outlier_subsets.R +++ b/R/generate_outlier_subsets.R @@ -38,7 +38,7 @@ #' dataset_log_transform = "eucalyptus") %>% #' generate_yi_subsets() %>% #TODO: must be run after prepare_response_variables?? #' apply_VZ_exclusions( -#' VZ_colname = list("eucalyptus" = "std.error_log", +#' VZ_colname = list("eucalyptus" = "se_log", #' "blue tit" = "VZ"), #' VZ_cutoff = 3) %>% #' generate_exclusion_subsets() %>% @@ -62,9 +62,10 @@ generate_outlier_subsets <- function(data, outcome_variable = NULL, n_min = NULL "estimate_type", "dataset") - if (!is.null(enexpr(ignore_subsets))) { - ignore_subsets_columns <- - rlang::call_args(enquo(ignore_subsets)) %>% + #TODO consider switching to exprs instead of list as input + # see meta_analyse_datasets + if (!is.null(enexpr(ignore_subsets))) { + ignore_subsets_columns <- rlang::call_args(enquo(ignore_subsets)) %>% map(rlang::f_lhs) %>% map(rlang::as_string) %>% list_c() %>% @@ -211,10 +212,14 @@ slice_conditionally <- function(data, n_min, n_max, outcome_variable) { #' @param filter_vars A list of quosures to be used in [dplyr::filter()] to subset `y` #' @param n_min integer, the number of bottom outliers to remove #' @param n_max integer, the number of top outliers to remove -#' @details the tibble `x` must contain the columns `data`, `outcome_colname`, `n_min`, and `n_max` +#' @details the tibble `x` must contain the columns `data`, `outcome_colname`, `n_min`, and `n_max`. Will create the columns `exclusion_set` if not present in the dataset, assigning `"complete"` to `x` and `"complete-rm_outliers"` to the subsetted data bound to `x`. #' @keywords internal apply_slice_conditionally <- function(x, filter_vars){ - out <- bind_rows(x, { + out <- bind_rows(x %>% + mutate( exclusion_set = + if ("exclusion_set" %in% colnames(.)) { + exclusion_set } + else {"complete"}), { x %>% filter(!!!filter_vars) %>% mutate(data = @@ -224,8 +229,10 @@ apply_slice_conditionally <- function(x, filter_vars){ n_max = ..4, outcome_variable = ..2 ))) %>% - mutate( - exclusion_set = paste0(exclusion_set, "-rm_outliers"), + mutate( exclusion_set = + if ("exclusion_set" %in% colnames(.)) { + paste0(exclusion_set, "-rm_outliers") } + else {"complete-rm_outliers"}, data = map2( .x = data, diff --git a/R/get_forest_plot_data.R b/R/get_forest_plot_data.R index a98d359..39b1b17 100644 --- a/R/get_forest_plot_data.R +++ b/R/get_forest_plot_data.R @@ -6,13 +6,14 @@ #' @family Plotting functions #' @examples #' get_forest_plot_data(model) -#' @import broom.mixed +#' @importFrom broom tidy #' @import dplyr -#' @import forcats -#' @import stringr +#' @importFrom forcats fct_reorder +#' @importFrom stringr str_detect +#' @importFrom purrr keep_at get_forest_plot_data <- function(model){ model %>% - broom.mixed::tidy(conf.int = TRUE, include_studies = TRUE) %>% + broom::tidy(conf.int = TRUE, include_studies = TRUE) %>% dplyr::mutate( point_shape = ifelse(stringr::str_detect(term, "overall"), @@ -39,8 +40,8 @@ get_forest_plot_data <- function(model){ #' @return A ggplot object #' @export #' @import ggplot2 -#' @import ggforestplot -#' @import NatParksPalettes +#' @importFrom ggforestplot theme_forest +#' @importFrom NatParksPalettes scale_color_natparks_d #' @family Plotting functions #' @examples #' data(ManyEcoEvo_results) diff --git a/R/log_transform.R b/R/log_transform.R index 84ccad9..c61e52f 100644 --- a/R/log_transform.R +++ b/R/log_transform.R @@ -3,12 +3,12 @@ #' @param estimate Point estimates on the original response scale #' @param std.error Standard error of the point estimates on the original response scale #' @param sim Number of simulations to run when generating distribution to sample from -#' @return A dataframe containing the `mean_log`, `std.error_log`, `lower` and `upper` 95% CI on the log scale +#' @return A dataframe containing the `mean_log`, `se_log`, `lower` and `upper` 95% CI on the log scale #' @export #' @examples log_transform(1.77, 1.01, 1000) #' @import dplyr #' @importFrom cli cli_alert_danger cli_alert_success -#' @import purrr +#' @importFrom purrr map_lgl flatten_dbl #' @seealso Equivalent in workflow data hierarchy to [Z_VZ_preds()]. log_transform <- function(estimate = numeric(1L), std.error = numeric(1L), @@ -23,7 +23,7 @@ log_transform <- function(estimate = numeric(1L), quantiles <- quantile(log_simulated, c(0.025, 0.975), na.rm = TRUE) out <- data.frame(mean_log = m_est, - std.error_log = std.error_est, + se_log = std.error_est, lower = quantiles[[1]], upper = quantiles[[2]]) diff --git a/R/log_transform_yi.R b/R/log_transform_yi.R index 7293c07..392587d 100644 --- a/R/log_transform_yi.R +++ b/R/log_transform_yi.R @@ -7,16 +7,20 @@ #' @param sim Number of simulations to run when generating distribution to sample from #' @export #' @import dplyr -#' @import purrr -#' @import cli -#' @import rlang +#' @importFrom purrr map2 +#' @importFrom tidyr unnest +#' @importFrom cli cli_warn +#' @importFrom rlang is_na log_transform_yi <- function(back_transformed_data, sim = 10000L, ...){ if (any(rlang::is_na(sim), rlang::is_na(back_transformed_data))) { + cli::cli_warn("Argument {.arg sim} or {.arg back_transformed_data} is {.val {NA}}. Returning {.val {NA}} for log-transformed predictions.") + return(NA) + } names_lookup <- c(estimate = "estimate", #blue tit diff --git a/R/make_param_table.R b/R/make_param_table.R index 96198c1..eddf269 100644 --- a/R/make_param_table.R +++ b/R/make_param_table.R @@ -11,7 +11,7 @@ #' make_param_table(ManyEcoEvo::blue_tit_data) #' make_param_table(ManyEcoEvo::euc_data) #' @import dplyr -#' @import tidyr +#' @importFrom tidyr pivot_longer #' @details The parameter table is used in the computation of Z-values during the standardisation of out-of-sample predictions with `pred_to_Zr()`. `make_param_table()` returns a tidy, long-form tibble with the variable names of `analysis_data` stored in column`variable`, the corresponding `parameter` ("mean" or "sd"), and the `value` of that `parameter`. #' #' @details # Note diff --git a/R/meta_analyse_datasets.R b/R/meta_analyse_datasets.R index 9f298ab..90da72a 100644 --- a/R/meta_analyse_datasets.R +++ b/R/meta_analyse_datasets.R @@ -44,7 +44,7 @@ #' #' Function assumes that if argument `outcome_variable` is supplied, then `outcome_SE` is also supplied, and conversely, if `outcome_SE` is not supplied, then neither is `outcome_variable` (*TODO* not yet checked in function). #' -meta_analyse_datasets <- function(data, outcome_variable = NULL, outcome_SE, filter_vars = NULL) { +meta_analyse_datasets <- function(data, outcome_variable = NULL, outcome_SE = NULL, filter_vars = NULL) { # ----- Argument Checks ----- stopifnot( @@ -129,18 +129,18 @@ meta_analyse_datasets <- function(data, outcome_variable = NULL, outcome_SE, fil matched_formulae_outcome_SE, .f = ~ map_match_formulae( {data %>% - ungroup %>% - select(names(outcome_SE))}, - .x, - .y, - col_name = "outcome_SE_colname")) %>% + ungroup %>% + select(names(outcome_SE))}, + .x, + .y, + col_name = "outcome_SE_colname")) %>% bind_rows() %>% drop_na(outcome_SE_colname) %>% distinct() }, by = unique(names(outcome_SE)) ) - + } # ----- Fit Meta-Models & Create Plots ----- @@ -171,19 +171,19 @@ meta_analyse_datasets <- function(data, outcome_variable = NULL, outcome_SE, fil effects_analysis = ifelse(is.na(MA_mod), NA, - purrr::map2( - .x = effects_analysis, - .y = MA_mod, - .f = ~ calculate_deviation_score(.x, .y) + purrr::pmap( + .l = list(effects_analysis, MA_mod, outcome_colname), + .f = calculate_deviation_score ) ), effects_analysis = ifelse(rlang::is_na(effects_analysis), NA, - purrr::map2( - .x = effects_analysis, - .y = dataset, - .f = ~ box_cox_transform(.x, .y) + purrr::pmap( + .l = list(effects_analysis, + dataset, + outcome_SE_colname), + .f = box_cox_transform ) ), sorensen_glm = @@ -259,8 +259,9 @@ meta_analyse_datasets <- function(data, outcome_variable = NULL, outcome_SE, fil } else { out <- out %>% - mutate(effects_analysis = map(effects_analysis, ~ .x %>% - unnest(review_data))) %>% + mutate(effects_analysis = + map(effects_analysis, ~ .x %>% + unnest(review_data))) %>% mutate(MA_mod_mv = map(effects_analysis, fit_multivar_MA)) %>% select(-ends_with("_colname")) } diff --git a/R/plot_model_means_box_cox_cat.R b/R/plot_model_means_box_cox_cat.R index 94bd4aa..3421277 100644 --- a/R/plot_model_means_box_cox_cat.R +++ b/R/plot_model_means_box_cox_cat.R @@ -11,7 +11,7 @@ #' @export #' @import ggplot2 #' @import dplyr -#' @import see +#' @importFrom see geom_jitter2 scale_fill_material_d theme_modern #' @importFrom EnvStats stat_n_text #' @importFrom forcats fct_relevel #' @importFrom sae bxcx diff --git a/R/prepare_response_variables.R b/R/prepare_response_variables.R index 4b9f3db..8be6a93 100644 --- a/R/prepare_response_variables.R +++ b/R/prepare_response_variables.R @@ -13,11 +13,12 @@ #' @family targets-pipeline functions. #' @family Multi-dataset Wrapper Functions #' @export -#' @import cli +#' @importFrom cli cli_abort cli_alert_info cli_alert_warning cli_h1 cli_h2 #' @import dplyr -#' @import purrr -#' @import rlang -#' @importFrom pointblank expect_col_exists expect_col_is_numeric expect_col_is_character expect_col_vals_in_set +#' @importFrom purrr set_names map2 pmap +#' @importFrom rlang !! := enquo +#' @importFrom pointblank expect_col_exists expect_col_is_numeric expect_col_is_character expect_col_vals_in_set test_col_exists +#' @importFrom tidyr tibble drop_na #' @examples #' data(ManyEcoEvo) #' ManyEcoEvo %>% prepare_response_variables(estimate_type = "Zr") @@ -33,48 +34,70 @@ prepare_response_variables <- function(ManyEcoEvo, several.ok = FALSE) stopifnot(is.data.frame(ManyEcoEvo)) + pointblank::expect_col_exists(object = ManyEcoEvo, columns = c(dataset, data)) if (!pointblank::test_col_exists(ManyEcoEvo, "estimate_type")) { + ManyEcoEvo <- dplyr::mutate(ManyEcoEvo, estimate_type = estimate_type) + } if (!is.null(dataset_standardise)) { + stopifnot( is.character(dataset_standardise), length(dataset_standardise) >= 1, - length(dataset_standardise) <= length(unique(ManyEcoEvo$dataset))) + length(dataset_standardise) <= length(unique(ManyEcoEvo$dataset)) + ) + match.arg(dataset_standardise, choices = ManyEcoEvo$dataset, several.ok = TRUE) } if (!is.null(dataset_log_transform)) { + stopifnot( is.character(dataset_log_transform), length(dataset_log_transform) >= 1, - length(dataset_log_transform) <= length(unique(ManyEcoEvo$dataset))) + length(dataset_log_transform) <= length(unique(ManyEcoEvo$dataset)) + ) + match.arg(dataset_log_transform, choices = ManyEcoEvo$dataset, several.ok = TRUE) } if (!is.null(param_table)) { + stopifnot(is.data.frame(param_table)) - pointblank::expect_col_exists(object = param_table, - columns = c("variable", - "parameter", - "value")) - pointblank::expect_col_is_numeric(object = param_table, - columns = "value") - pointblank::expect_col_is_character(object = param_table, - columns = c("variable", - "parameter")) - pointblank::expect_col_vals_in_set(object = param_table, - columns = "parameter", - set = c("mean", "sd")) + + pointblank::expect_col_exists( + object = param_table, + columns = c("variable", + "parameter", + "value") + ) + + pointblank::expect_col_is_numeric( + object = param_table, + columns = "value" + ) + + pointblank::expect_col_is_character( + object = param_table, + columns = c("variable", + "parameter") + ) + + pointblank::expect_col_vals_in_set( + object = param_table, + columns = "parameter", + set = c("mean", "sd") + ) } # ------ Prepare Response Variables ------ @@ -82,8 +105,11 @@ prepare_response_variables <- function(ManyEcoEvo, out <- ManyEcoEvo if (estimate_type != "Zr") { + if (all(is.null(param_table), !is.null(dataset_standardise))) { + cli::cli_abort("{.arg param_table} must be supplied for {.val {estimate_type}} data when {.arg dataset_standardise} is not {.val NULL}.") + } # ------ Back transform if estimate_type is yi only ------ @@ -96,9 +122,7 @@ prepare_response_variables <- function(ManyEcoEvo, .y = dataset, .f = ~ back_transform_response_vars_yi( dat = .x, - estimate_type = !!{ - estimate_type - }, + estimate_type = !!{estimate_type}, dataset = .y ) ), @@ -126,7 +150,9 @@ prepare_response_variables <- function(ManyEcoEvo, ) } else { #yi + standardise and/or log-transform + cli::cli_alert_info("Standardising and/or log-transforming response variables for {.val {estimate_type}} estimates.") + transform_datasets <- bind_rows( tibble( @@ -140,9 +166,12 @@ prepare_response_variables <- function(ManyEcoEvo, set_names("standardise_response") )) %>% drop_na() # in case of NULLs + } } else { # Zr + if (!is.null(param_table)) { + cli::cli_abort("{.arg param_table} must be {.val NULL} for {.val {estimate_type}} data") } @@ -170,8 +199,10 @@ prepare_response_variables <- function(ManyEcoEvo, select(-fns) return(out) + } +#' Wrapper function to standardise response variables #' @rdname process_analyst_data #' @seealso This internal helper function is called by [prepare_response_variables()] #' @importFrom rlang caller_env exec diff --git a/R/standardise_response.R b/R/standardise_response.R index d2305b3..a9e3720 100644 --- a/R/standardise_response.R +++ b/R/standardise_response.R @@ -4,7 +4,7 @@ #' @param ... Ignored #' @importFrom cli cli_h1 cli_h2 #' @import dplyr -#' @import tidyr +#' @importFrom tidyr unnest #' @importFrom purrr pmap map2 map_int map #' @importFrom rlang is_na #' @importFrom pointblank col_exists vars @@ -54,6 +54,9 @@ standardise_response <- function(data, param_table = NULL, dataset = character(1L), ...) { + stopifnot(is.data.frame(data), + is.character(dataset)) + # TODO insert checks that appropriate columns exist # TODO apply to data and check that all cases accounted for! match.arg(estimate_type, choices = c("Zr", "yi", "y25", "y50", "y75"), several.ok = FALSE) @@ -63,7 +66,9 @@ standardise_response <- function(data, "{.val {estimate_type}}")) if (estimate_type == "Zr") { + # ------ Convert Effect Sizes to Zr ------- + cli::cli_h2(paste0("Computing standardised effect sizes ", "{.code Zr}", " and variance ", "{.code VZr}")) data <- data %>% @@ -78,7 +83,9 @@ standardise_response <- function(data, )) %>% tidyr::unnest(cols = c(Zr_VZr)) } else { + # ------ Convert predictions to Z ------- + cli::cli_h2(paste0("Standardising out-of-sample predictions")) data <- data %>% @@ -145,8 +152,7 @@ standardise_response <- function(data, #' and for eucalyptus data, `data$back_transformed_data$estimate` is renamed `Z`. #' `se.fit` is renamed `VZ`. #' @import dplyr -#' @import purrr -#' @importFrom tidyr any_of +#' @importFrom purrr map process_response <- function(data, ...){ Z_names_lookup <- c(Z = "estimate", #blue tit @@ -180,7 +186,7 @@ process_response <- function(data, ...){ #' @export #' @import dplyr #' @importFrom cli cli_h1 cli_h2 -#' @importFrom pointblank col_exists vars +#' @importFrom pointblank col_exists #' @importFrom purrr map #' @importFrom rlang is_na #' @describeIn process_analyst_data Standardise response data for meta-analysis diff --git a/R/summarise_analysis_types.R b/R/summarise_analysis_types.R index f5f1ecc..5755b9e 100644 --- a/R/summarise_analysis_types.R +++ b/R/summarise_analysis_types.R @@ -13,16 +13,17 @@ #' @return A summarised tibble with the variables `subset`, `dataset`, `num_teams`, `total_analyses`, `sum_linear`, `sum_mixed`, `sum_Bayesian`. #' @export #' @import dplyr -#' @import broom -#' @import tidyr +#' @importFrom broom tidy +#' @importFrom tidyr unnest #' @import metafor -#' @import purrr +#' @importFrom purrr map set_names map_dfr #' @family Multi-dataset Wrapper Functions #' @author Hannah S. Fraser #' @author Elliot Gould #' @examples #' summarise_analysis_types(ManyEcoEvo_results, ManyEcoEvo_yi_results, ManyEcoEvo) summarise_analysis_types <- function(ManyEcoEvo_results, ManyEcoEvo_yi_results, ManyEcoEvo) { + effect_ids <- ManyEcoEvo_results %>% filter( exclusion_set == "complete", @@ -30,21 +31,22 @@ summarise_analysis_types <- function(ManyEcoEvo_results, ManyEcoEvo_yi_results, ) %>% select(MA_mod, effects_analysis) %>% group_by(estimate_type, dataset) %>% - mutate(tidy_mod = map( - MA_mod, - ~ broom::tidy(.x, - conf.int = TRUE, - include_studies = TRUE - ) %>% - rename(study_id = term) - ), .keep = "none") %>% + mutate(tidy_mod = + map( + MA_mod, + ~ broom::tidy(.x, + conf.int = TRUE, + include_studies = TRUE + ) %>% + rename(study_id = term) + ), .keep = "none") %>% unnest(tidy_mod) %>% filter(type == "study") %>% ungroup() %>% select(study_id) %>% rename(id_col = study_id) %>% # TODO duplicates for "Bell-2-2-1" and "Bonalbo-1-1-1 WHY? distinct() - + prediction_ids <- ManyEcoEvo_yi_results %>% # TODO Euc mod_data_logged not here! filter( exclusion_set == "complete", @@ -52,18 +54,19 @@ summarise_analysis_types <- function(ManyEcoEvo_results, ManyEcoEvo_yi_results, ) %>% select(MA_mod, effects_analysis, -exclusion_set) %>% group_by(estimate_type, dataset) %>% - mutate(tidy_mod = map( - MA_mod, - ~ broom::tidy(.x, conf.int = TRUE, include_studies = TRUE) %>% - rename(study_id = term) - ), .keep = "none") %>% + mutate(tidy_mod = + map( + MA_mod, + ~ broom::tidy(.x, conf.int = TRUE, include_studies = TRUE) %>% + rename(study_id = term) + ), .keep = "none") %>% unnest(tidy_mod) %>% filter(type == "study") %>% ungroup() %>% select(study_id) %>% rename(id_col = study_id) %>% distinct() - + Master <- ManyEcoEvo %>% ungroup() %>% @@ -84,13 +87,13 @@ summarise_analysis_types <- function(ManyEcoEvo_results, ManyEcoEvo_yi_results, lm = ifelse(linear_model == "linear", 1, 0), glm = ifelse(linear_model == "generalised", 1, 0) ) # TODO move this into master processing so don't have to repeat else where!! - + effects <- Master %>% right_join(effect_ids, by = c("id_col")) # repeat for each - + predictions <- Master %>% right_join(prediction_ids, by = c("id_col")) - + summarised_data <- full_join( map_dfr( .x = list(effects, predictions) %>% @@ -105,12 +108,11 @@ summarise_analysis_types <- function(ManyEcoEvo_results, ManyEcoEvo_yi_results, .id = "subset" ) ) - + return(summarised_data) # TODO next: set up so can run on just one object ManyEcoEvo_results, and account for subsets too! } - #' Summarise number of analyst teams and total analyses per dataset #' #' @param data A dataframe containing the variables `TeamIdentifier` and `dataset` diff --git a/R/summarise_conclusions.R b/R/summarise_conclusions.R index 7b47a19..a795c92 100644 --- a/R/summarise_conclusions.R +++ b/R/summarise_conclusions.R @@ -9,9 +9,9 @@ #' @return A dataframe with count values for each unique `Conclusion` in columns for each `subset` ("effects", "predictions", "all"), for each `dataset`. #' @export #' @import dplyr -#' @import purrr -#' @import broom -#' @import tidyr +#' @importFrom purrr map_dfr map +#' @importFrom broom tidy +#' @importFrom tidyr unnest pivot_wider #' @import metafor #' @family Multi-dataset Wrapper Functions #' @author Hannah S. Fraser diff --git a/R/summarise_model_composition.R b/R/summarise_model_composition.R index cdd0b8c..1efd364 100644 --- a/R/summarise_model_composition.R +++ b/R/summarise_model_composition.R @@ -5,7 +5,6 @@ #' @details #' Generates summary data with [summarise_model_composition_data()]. #' -#' #' @param ManyEcoEvo_results A tibble of `ManyEcoEvo_results` #' @param ManyEcoEvo_yi_results A tibble of `ManyEcoEvo_yi_results` #' @param ManyEcoEvo A tibble of `ManyEcoEvo` @@ -13,8 +12,9 @@ #' @return A dataframe #' @export #' @import dplyr -#' @import tidyr -#' @import purrr +#' @importFrom broom tidy +#' @importFrom tidyr unnest pivot_longer +#' @importFrom purrr map map_dfr set_names #' @author Hannah S. Fraser #' @author Elliot Gould #' @family Multi-dataset Wrapper Functions @@ -28,21 +28,22 @@ summarise_model_composition <- function(ManyEcoEvo_results, ManyEcoEvo_yi_result ) %>% select(MA_mod, effects_analysis) %>% group_by(estimate_type, dataset) %>% - mutate(tidy_mod = map( - MA_mod, - ~ broom::tidy(.x, - conf.int = TRUE, - include_studies = TRUE - ) %>% - rename(study_id = term) - ), .keep = "none") %>% + mutate(tidy_mod = + map( + MA_mod, + ~ broom::tidy(.x, + conf.int = TRUE, + include_studies = TRUE + ) %>% + rename(study_id = term) + ), .keep = "none") %>% unnest(tidy_mod) %>% filter(type == "study") %>% ungroup() %>% select(study_id) %>% rename(id_col = study_id) %>% # TODO duplicates for "Bell-2-2-1" and "Bonalbo-1-1-1 WHY? distinct() - + prediction_ids <- ManyEcoEvo_yi_results %>% # TODO Euc mod_data_logged not here! filter( exclusion_set == "complete", @@ -50,19 +51,20 @@ summarise_model_composition <- function(ManyEcoEvo_results, ManyEcoEvo_yi_result ) %>% select(MA_mod, effects_analysis, -exclusion_set) %>% group_by(estimate_type, dataset) %>% - mutate(tidy_mod = map( - MA_mod, - ~ broom::tidy(.x, conf.int = TRUE, include_studies = TRUE) %>% - rename(study_id = term) - ), .keep = "none") %>% + mutate(tidy_mod = + map( + MA_mod, + ~ broom::tidy(.x, conf.int = TRUE, include_studies = TRUE) %>% + rename(study_id = term) + ), .keep = "none") %>% unnest(tidy_mod) %>% filter(type == "study") %>% ungroup() %>% select(study_id) %>% rename(id_col = study_id) %>% distinct() - - + + Master <- ManyEcoEvo %>% ungroup() %>% @@ -83,13 +85,13 @@ summarise_model_composition <- function(ManyEcoEvo_results, ManyEcoEvo_yi_result lm = ifelse(linear_model == "linear", 1, 0), glm = ifelse(linear_model == "generalised", 1, 0) ) # TODO move this into master processing so don't have to repeat else where!! - + effects <- Master %>% # TODO consider generalising, so for each 'estimate_type' group: repeat. right_join(effect_ids, by = c("id_col")) # repeat for each - + predictions <- Master %>% right_join(prediction_ids, by = c("id_col")) - + summarised_data <- map_dfr( .x = list(effects, predictions, Master) %>% @@ -103,7 +105,7 @@ summarise_model_composition <- function(ManyEcoEvo_results, ManyEcoEvo_yi_result names_pattern = "(.*)_(.*)", values_to = "value" ) - + return(summarised_data) } diff --git a/R/summarise_review_data.R b/R/summarise_review_data.R index 18859bf..0eef712 100644 --- a/R/summarise_review_data.R +++ b/R/summarise_review_data.R @@ -9,9 +9,9 @@ #' @return A named list of length two of `summarised_data_analyses` containing a tibble of summary statistics for each outcome `subset` (`effects` or `predictions`) generated by [summarise_analyses_by_reviewer()], and `summarised_data_reviews` containing a tibble of sumary statistics for each outcome `subset`, for each `dataset` generated by [summarise_reviews_per_analysis()]. #' @export #' @import dplyr -#' @import tidyr -#' @import purrr -#' @import broom +#' @importFrom tidyr unnest +#' @importFrom purrr set_names map_dfr +#' @importFrom broom tidy #' @import metafor #' @author Hannah S. Fraser #' @author Elliot Gould @@ -30,7 +30,7 @@ summarise_reviews <- function(ManyEcoEvo, ManyEcoEvo_results, ManyEcoEvo_yi_resu unnest(data) %>% select(ends_with("_id"), id_col, dataset, review_data) %>% unnest(review_data) - + effect_ids <- ManyEcoEvo_results %>% filter( exclusion_set == "complete", @@ -38,21 +38,22 @@ summarise_reviews <- function(ManyEcoEvo, ManyEcoEvo_results, ManyEcoEvo_yi_resu ) %>% select(MA_mod, effects_analysis) %>% group_by(estimate_type, dataset) %>% - mutate(tidy_mod = map( - MA_mod, - ~ broom::tidy(.x, - conf.int = TRUE, - include_studies = TRUE - ) %>% - rename(study_id = term) - ), .keep = "none") %>% + mutate(tidy_mod = + map( + MA_mod, + ~ broom::tidy(.x, + conf.int = TRUE, + include_studies = TRUE + ) %>% + rename(study_id = term) + ), .keep = "none") %>% unnest(tidy_mod) %>% filter(type == "study") %>% ungroup() %>% select(study_id) %>% rename(id_col = study_id) %>% # TODO duplicates for "Bell-2-2-1" and "Bonalbo-1-1-1 WHY? distinct() - + prediction_ids <- ManyEcoEvo_yi_results %>% # TODO Euc mod_data_logged not here! filter( exclusion_set == "complete", @@ -60,39 +61,40 @@ summarise_reviews <- function(ManyEcoEvo, ManyEcoEvo_results, ManyEcoEvo_yi_resu ) %>% select(MA_mod, effects_analysis, -exclusion_set) %>% group_by(estimate_type, dataset) %>% - mutate(tidy_mod = map( - MA_mod, - ~ broom::tidy(.x, conf.int = TRUE, include_studies = TRUE) %>% - rename(study_id = term) - ), .keep = "none") %>% + mutate(tidy_mod = + map( + MA_mod, + ~ broom::tidy(.x, conf.int = TRUE, include_studies = TRUE) %>% + rename(study_id = term) + ), .keep = "none") %>% unnest(tidy_mod) %>% filter(type == "study") %>% ungroup() %>% select(study_id) %>% rename(id_col = study_id) %>% distinct() - + effects <- review_data %>% # TODO consider generalising, so for each 'estimate_type' group: repeat. right_join(effect_ids, by = c("id_col")) %>% # repeat for each filter(!is.na(dataset)) # TODO id source of NAs in dataset! - + predictions <- review_data %>% right_join(prediction_ids, by = c("id_col")) - + summarised_data_analyses <- map_dfr( .x = list(effects, predictions) %>% purrr::set_names("effects", "predictions"), .f = summarise_analyses_by_reviewer, .id = "subset" ) - + summarised_data_reviews <- map_dfr( .x = list(effects, predictions) %>% purrr::set_names("effects", "predictions"), .f = summarise_reviews_per_analysis, .id = "subset" ) - + summarised_data <- list( summarised_data_analyses, summarised_data_reviews @@ -101,7 +103,7 @@ summarise_reviews <- function(ManyEcoEvo, ManyEcoEvo_results, ManyEcoEvo_yi_resu "summarised_data_analyses", "summarised_data_reviews" ) - + return(summarised_data) } diff --git a/R/summarise_sorensen_index.R b/R/summarise_sorensen_index.R index 733d66c..f8c962a 100644 --- a/R/summarise_sorensen_index.R +++ b/R/summarise_sorensen_index.R @@ -10,8 +10,9 @@ #' @return A tibble of aggregate summary statistics (`mean`, `sd`, `min`, `max`) for mean Sorensen's index estimates across each `subset` and `dataset`. #' @export #' @import dplyr -#' @import tidyr -#' @import purrr +#' @importFrom tidyr unnest +#' @importFrom purrr map map_dfr set_names +#' @importFrom broom tidy #' @author Hannah S. Fraser #' @author Elliot Gould #' @family Multi-dataset Wrapper Functions @@ -28,14 +29,14 @@ summarise_sorensen_index <- function(ManyEcoEvo_results, ManyEcoEvo_yi_results) ungroup() %>% select(dataset, diversity_indices) %>% unnest(diversity_indices) - + sorensen_index_yi <- ManyEcoEvo_yi_results %>% filter(exclusion_set == "complete") %>% ungroup() %>% select(dataset, diversity_indices) %>% unnest(diversity_indices) - + effect_ids <- ManyEcoEvo_results %>% # TODO ensure for other related functions are properly filtering ManyEcoEvo_results filter( exclusion_set == "complete", @@ -44,43 +45,45 @@ summarise_sorensen_index <- function(ManyEcoEvo_results, ManyEcoEvo_yi_results) ) %>% select(MA_mod, effects_analysis) %>% group_by(estimate_type, dataset) %>% - mutate(tidy_mod = map( - MA_mod, - ~ broom::tidy(.x, - conf.int = TRUE, - include_studies = TRUE - ) %>% - rename(study_id = term) - ), .keep = "none") %>% + mutate(tidy_mod = + map( + MA_mod, + ~ broom::tidy(.x, + conf.int = TRUE, + include_studies = TRUE + ) %>% + rename(study_id = term) + ), .keep = "none") %>% unnest(tidy_mod) %>% filter(type == "study") %>% ungroup() %>% select(study_id) %>% rename(id_col = study_id) %>% # TODO duplicates for "Bell-2-2-1" and "Bonalbo-1-1-1 WHY? distinct() - + prediction_ids <- ManyEcoEvo_yi_results %>% # TODO Euc mod_data_logged not here! filter(exclusion_set == "complete") %>% select(MA_mod, effects_analysis, -exclusion_set) %>% group_by(estimate_type, dataset) %>% - mutate(tidy_mod = map( - MA_mod, - ~ broom::tidy(.x, conf.int = TRUE, include_studies = TRUE) %>% - rename(study_id = term) - ), .keep = "none") %>% + mutate(tidy_mod = + map( + MA_mod, + ~ broom::tidy(.x, conf.int = TRUE, include_studies = TRUE) %>% + rename(study_id = term) + ), .keep = "none") %>% unnest(tidy_mod) %>% filter(type == "study") %>% ungroup() %>% select(study_id) %>% rename(id_col = study_id) %>% distinct() - + effects <- sorensen_index_zr %>% # TODO consider generalising, so for each 'estimate_type' group: repeat. right_join(effect_ids, by = c("id_col")) # repeat for each - + predictions <- sorensen_index_yi %>% right_join(prediction_ids, by = c("id_col")) - + summarised_data <- map_dfr( .x = list(predictions, effects) %>% @@ -88,7 +91,7 @@ summarise_sorensen_index <- function(ManyEcoEvo_results, ManyEcoEvo_yi_results) .f = summarise_sorensen_index_data, .id = "subset" ) - + return(summarised_data) } @@ -99,7 +102,7 @@ summarise_sorensen_index <- function(ManyEcoEvo_results, ManyEcoEvo_yi_results) #' @param data A dataframe containing `mean_diversity_index` for the Sorensen's index estimates for each analysis `id_col`, for each `dataset`. #' @return A dataframe with the `mean`, `sd`, `min`, `max` mean Sorensen's index values for each `dataset`. #' @export -#' +#' @import dplyr #' @examples #' ManyEcoEvo_results %>% #' filter( diff --git a/R/summarise_variable_counts.R b/R/summarise_variable_counts.R index 1bb6b26..3a86040 100644 --- a/R/summarise_variable_counts.R +++ b/R/summarise_variable_counts.R @@ -13,9 +13,9 @@ #' @return A dataframe of count values `n` or summary statistic values `n_mean`,`n_sd`,`n_min`,`n_max` of counts depending on the the value supplied to the `output` argument. #' @export #' @import dplyr -#' @import tidyr -#' @import purrr -#' @import broom +#' @importFrom broom tidy +#' @importFrom tidyr unnest +#' @importFrom purrr map map_dfr set_names #' @author Hannah S. Fraser #' @author Elliot Gould #' @family Multi-dataset Wrapper Functions @@ -25,11 +25,11 @@ #' summarise_variable_counts(ManyEcoEvo, ManyEcoEvo_results, ManyEcoEvo_yi_results, "aggregate") summarise_variable_counts <- function(ManyEcoEvo, ManyEcoEvo_results, ManyEcoEvo_yi_results, output = "count") { stopifnot(output == "count" | output == "aggregate") - + diversity <- ManyEcoEvo %>% select(diversity_data) %>% unnest(diversity_data) - + effect_ids <- ManyEcoEvo_results %>% filter( exclusion_set == "complete", @@ -38,14 +38,15 @@ summarise_variable_counts <- function(ManyEcoEvo, ManyEcoEvo_results, ManyEcoEvo select(MA_mod, effects_analysis) %>% group_by(estimate_type, dataset) %>% mutate( - tidy_mod = map( - MA_mod, - ~ broom::tidy(.x, - conf.int = TRUE, - include_studies = TRUE - ) %>% - rename(study_id = term) - ), + tidy_mod = + map( + MA_mod, + ~ broom::tidy(.x, + conf.int = TRUE, + include_studies = TRUE + ) %>% + rename(study_id = term) + ), .keep = "none" ) %>% unnest(tidy_mod) %>% @@ -54,7 +55,7 @@ summarise_variable_counts <- function(ManyEcoEvo, ManyEcoEvo_results, ManyEcoEvo select(study_id) %>% rename(id_col = study_id) %>% # TODO duplicates for "Bell-2-2-1" and "Bonalbo-1-1-1 WHY? distinct() - + prediction_ids <- ManyEcoEvo_yi_results %>% # TODO Euc mod_data_logged not here! filter( exclusion_set == "complete", @@ -62,24 +63,25 @@ summarise_variable_counts <- function(ManyEcoEvo, ManyEcoEvo_results, ManyEcoEvo ) %>% select(MA_mod, effects_analysis, -exclusion_set) %>% group_by(estimate_type, dataset) %>% - mutate(tidy_mod = map( - MA_mod, - ~ broom::tidy(.x, conf.int = TRUE, include_studies = TRUE) %>% - rename(study_id = term) - ), .keep = "none") %>% + mutate(tidy_mod = + map( + MA_mod, + ~ broom::tidy(.x, conf.int = TRUE, include_studies = TRUE) %>% + rename(study_id = term) + ), .keep = "none") %>% unnest(tidy_mod) %>% filter(type == "study") %>% ungroup() %>% select(study_id) %>% rename(id_col = study_id) %>% distinct() - + effects <- diversity %>% # TODO consider generalising, so for each 'estimate_type' group: repeat. right_join(effect_ids, by = c("id_col")) # repeat for each - + predictions <- diversity %>% right_join(prediction_ids, by = c("id_col")) - + # repeat for all, predictions, effects summarised_data <- map_dfr( @@ -87,22 +89,27 @@ summarise_variable_counts <- function(ManyEcoEvo, ManyEcoEvo_results, ManyEcoEvo purrr::set_names("effects", "predictions", "all"), .f = count_analyses_variables_used, .id = "subset" ) - + if (output == "count") { + return(summarised_data) + } else { + summarised_data <- summarised_data %>% group_by(dataset, subset) %>% summarise(across(n, - .fns = list( - mean = ~ mean(.x, na.rm = TRUE), - sd = ~ sd(.x, na.rm = TRUE), - min = ~ min(.x, na.rm = TRUE), # TODO is value of 0 correct? - max = ~ max(.x, na.rm = TRUE) - ) + .fns = list( + mean = ~ mean(.x, na.rm = TRUE), + sd = ~ sd(.x, na.rm = TRUE), + min = ~ min(.x, na.rm = TRUE), # TODO is value of 0 correct? + max = ~ max(.x, na.rm = TRUE) + ) )) } + return(summarised_data) + } #' Count number of analyses each variable is used @@ -115,7 +122,7 @@ summarise_variable_counts <- function(ManyEcoEvo, ManyEcoEvo_results, ManyEcoEvo #' @return A dataframe of counts `n` each `variable` is used across all analyses within a given `dataset` #' @export #' @import dplyr -#' @import tidyr +#' @importFrom tidyr pivot_longer #' @author Hannah S. Fraser #' @author Elliot Gould #' @examples diff --git a/R/utils.R b/R/utils.R index ad64226..dd5933a 100644 --- a/R/utils.R +++ b/R/utils.R @@ -134,8 +134,8 @@ rm_inf_na <- function(effects_analysis, Z_colname, VZ_colname) { #' @return A named list of tibbles. Each tibble contains the rows of .tbl for the associated group and all the columns, including the grouping variables. Note that this returns a list_of which is slightly stricter than a simple list but is useful for representing lists where every element has the same type. #' @export #' @import dplyr -#' @import rlang -#' @import purrr +#' @importFrom rlang ensym +#' @importFrom purrr set_names map_chr pluck map #' @examples #' data(euc_data) #' data(blue_tit_data)