R/mean-by-group.R

In lifecycle-project/ds-helper: Helper Functions for Use with 'DataSHIELD'

#' Describes a numeric variable by strata of another numeric grouping variable.
#'
#' This has similar functionality to `tapply` or the dplyr chain `group_by` %>%
#' `summarise`. It offers additional flexilibity over `ds.tapply` in that it
#' allows you to specify upper and lower values for each strata. By contrast,
#' ds.tapply will produce a summary based on every unique value of the grouping
#' variable, which may not always be what is required.
#'
#' @importFrom DSI datashield.connections_find datashield.assign
#' @importFrom dsBaseClient ds.meanSdGp
#' @importFrom purrr map
#' @importFrom tibble as_tibble
#' @importFrom dplyr bind_rows %>% slice 
#' @importFrom tidyr separate
#'
#' @template conns
#' @template df
#' @param outcome String specifying outcome variable within `df`.
#' @param group_var String specifying grouping variable within 'df'.
#' @param intervals Optionally, numeric vector defining how to stratify
#' `group_var`. Values should specify alternately lower and upper values for
#' each strata. If NULL, `outcome` is summarised by every unique value of
#' `group_var`.
#' @template checks
#'
#' @return Tibble containing mean values for each strata of `group_var`.
#'
#' @family descriptive functions
#'
#' @export
dh.meanByGroup <- function(df = NULL, outcome = NULL, group_var = NULL,
                           intervals = NULL, conns = NULL, checks = FALSE) {
  value <- op <- tmp <- varname <- new_df_name <- age <- group <- cohort <-
    . <- enough_obs <- variable <- level <- std.dev <- nvalid <- ntotal <-
    x <- nvalid_2 <- nvalid_1 <- mean_2 <- std.dev_2 <- nmissing <- NULL

  if (is.null(df)) {
    stop("`df` must not be NULL.", call. = FALSE)
  }

  if (is.null(outcome)) {
    stop("`outcome` must not be NULL.", call. = FALSE)
  }

  if (is.null(group_var)) {
    stop("`group_var` must not be NULL.", call. = FALSE)
  }

  if (is.null(conns)) {
    conns <- datashield.connections_find()
  }

  if (checks == TRUE) {
    .isDefined(df = df, vars = vars, conns = conns)
  }

  ## There is an easy way and a hard way. If we bin based on integer units of
  ## the binning variable it is quite quick

  if (is.null(intervals)) {


    ## ---- First we round up the age variable -----------------------------------------------
    DSI::datashield.assign(
      conns, "age_tmp", as.symbol(paste0(df, "$", group_var, "+0.5"))
    )

    calltext <- call("asIntegerDS", "age_tmp")
    DSI::datashield.assign(conns, "age_round", calltext)

    ## ---- Now we get the mean values for the outcome by age --------------------------------
    calltext <- paste0(
      "meanSdGpDS(", paste0(df, "$", outcome), ",",
      "age_round", ")"
    )
    
    mean_tmp <- DSI::datashield.aggregate(conns, as.symbol(calltext))
    ## ---- Deal with disclosure issues ----------------------------------------
    warnings <- .checkDisclosureMeanGp(mean_tmp)
    
    if(nrow(warnings$issues) > 0){
      
      warning(
"The following cohorts have at least one group with between 1 and 2 observations 
so it is not possible to return grouped values. Try collapsing some categories or 
re-run the function using the `intervals` argument. \n\n", 
      
paste0(warnings$issues$cohort, collapse = ", ")

)
      
    }
    
    ## ---- Now put into neat long format -------------------------------------------------------
    out <- warnings$no_issues_data %>%
      map(function(x) {
        x$Mean_gp %>%
          as_tibble(rownames = "age") %>%
          mutate(age = as.numeric(str_remove(age, "age_round_")))
      }) %>%
      bind_rows(.id = "cohort") %>%
      dplyr::rename(mean = value)

    ## ---- Remove temporary objects ------------------------------------------------------------
    dh.tidyEnv(
      obj = c("age_tmp", "age_round"),
      type = "remove", 
      conns = conns
    )
  } else if (!is.null(intervals)) {

    ## ---- This is the harder one -------------------------------------------------------------------

    ## First we need to create a table defining our age bands.
    pairs <- split(intervals, ceiling(seq_along(intervals) / 2))

    subnames <- unlist(
      pairs %>% map(~ paste0("grp", "_", paste0(., collapse = "_"))),
      use.names = FALSE
    )

    cats <- tibble(
      varname = rep(subnames, each = 2),
      value = intervals,
      op = rep(c(">=", "<="), times = (length(intervals) / 2)),
      tmp = ifelse(op == ">=", "gte", "lte"),
      new_df_name = paste0("grp", tmp, value)
    )


    ## Now we create vectors of 1s and 0s indicating whether the above criteria
    ## are met for each subject.
    cats %>%
      pmap(function(value, op, new_df_name, ...) {
        ds.Boole(
          V1 = paste0(df, "$", group_var),
          V2 = value,
          Boolean.operator = op,
          newobj = new_df_name,
          datasources = conns
        )
      })

    suppressMessages(
      assign_conditions <- cats %>%
        group_by(varname) %>%
        summarise(condition = paste(new_df_name, collapse = "*"))
    )

    assign_conditions %>%
      pmap(function(condition, varname) {
        DSI::datashield.assign(conns, varname, as.symbol(condition))
      })
    
    ## Check disclosure issues
    disclosure <- assign_conditions$varname %>%
      map(~.checkDisclosure(., conns = conns)) %>%
      set_names(assign_conditions$varname) %>%
      bind_rows(.id = "variable")
    

    no_obs <- disclosure %>% dplyr::filter(enough_obs == FALSE)
    
    if(nrow(no_obs) > 0){
      
      warning(paste(capture.output({
        cat("These cohorts have no observations within the following
              intervals and will be excluded\n\n")
        no_obs %>% dplyr::select(variable, cohort) %>%
          print(n = Inf)
      }), collapse = "\n"))

    }
    
    valid_obs <- disclosure %>% dplyr::filter(enough_obs == TRUE)
    
    ## We then use these vectors to summarise mean observed height at the age
    ## periods we are interested in.
    obs_by_agecat_comb <- valid_obs %>%
      pmap(function(variable, cohort, ...){
        
        ds.meanSdGp(
          x = paste0(df, "$", outcome),
          y = variable,
          type = "split",
          datasources = conns[cohort])
      })
    

    ## Now we take these values and put them into a neater table
    out <- obs_by_agecat_comb %>%
      map(function(x){
        
        tibble(
          group = dimnames(x$Mean_gp_study)[[1]],
          mean = as.numeric(x$Mean_gp_study),
          std.dev = as.numeric(x$StDev_gp_study), 
          nvalid = as.numeric(x$Nvalid_gp_study),
          cohort = colnames(x$Mean_gp_study)) 
      }) %>%
      bind_rows() %>%
      separate(group, into = c("group", "level"), sep="_(?=[^_]+$)") %>%
      pivot_wider(
        names_from = "level", 
        values_from = c("nvalid", "mean", "std.dev")) %>%
      dplyr::rename(
        "nvalid" = nvalid_2, 
        "nmissing" = nvalid_1, 
        "mean" = mean_2,
        "std.dev" = std.dev_2) %>%
      mutate(group = str_remove(group, "grp_")) %>%
      dplyr::select(cohort, group, mean, std.dev, nvalid, nmissing)
    
        
    ## ---- Remove temporary objects -------------------------------------------
    dh.tidyEnv(
      obj = c(cats$new_df_name, assign_conditions$varname),
      type = "remove", 
      conns = conns
    )
  }

  return(out)
}


#' Internal function to check disclosure issues in this function
#' 
#' @param mean_gp output from ds.meanSdGp
#' 
#' @importFrom purrr map
#' @importFrom tibble tibble
#' @importFrom dplyr %>% bind_rows filter
#' @noRd
.checkDisclosureMeanGp <- function(mean_gp){
  
  warnings <- mean_gp %>% map(function(x){
    tibble(warning = !is.null(x$Warning))}) %>%
    bind_rows(.id = "cohort")
  
  issues <- warnings %>% dplyr::filter(warning == TRUE)  
  no_issues <- mean_gp[!names(mean_gp) %in% issues$cohort] 
  
  return(list(
    summary = warnings, 
    issues = issues, 
    no_issues_data = no_issues))
}