R/FamiliarDataComputationCalibrationData.R

In familiar: End-to-End Automated Machine Learning and Model Evaluation

#' @include FamiliarS4Generics.R
#' @include FamiliarS4Classes.R
#' @include PredictionTable.R
NULL


# familiarDataElementCalibrationData object ------------------------------------
setClass(
  "familiarDataElementCalibrationData",
  contains = "familiarDataElement"
)

# familiarDataElementCalibrationLinearFit object -------------------------------
setClass(
  "familiarDataElementCalibrationLinearFit",
  contains = "familiarDataElement"
)

# familiarDataElementCalibrationGoodnessOfFit object ---------------------------
setClass(
  "familiarDataElementCalibrationGoodnessOfFit",
  contains = "familiarDataElement"
)

# familiarDataElementCalibrationDensity object ---------------------------------
setClass(
  "familiarDataElementCalibrationDensity",
  contains = "familiarDataElement"
)



# extract_calibration_data (generic) -------------------------------------------

#'@title Internal function to extract calibration data.
#'
#'@description Computes calibration data from a `familiarEnsemble` object.
#'  Calibration tests are performed based on expected (predicted) and observed
#'  outcomes. For all outcomes, calibration-at-the-large and calibration slopes
#'  are determined. Furthermore, for all but survival outcomes, a repeated,
#'  randomised grouping Hosmer-Lemeshow test is performed. For survival
#'  outcomes, the Nam-D'Agostino and Greenwood-Nam-D'Agostino tests are
#'  performed.
#'
#'@inheritParams .extract_data
#'
#'@return A list with data.tables containing calibration test information for
#'  the ensemble model.
#'@md
#'@keywords internal
setGeneric(
  "extract_calibration_data",
  function(
    object,
    data,
    cl = NULL,
    ensemble_method = waiver(),
    evaluation_times = waiver(),
    detail_level = waiver(),
    estimation_type = waiver(),
    aggregate_results = waiver(),
    confidence_level = waiver(),
    bootstrap_ci_method = waiver(),
    is_pre_processed = FALSE,
    message_indent = 0L,
    verbose = FALSE,
    ...
  ) {
    standardGeneric("extract_calibration_data")
  }
)


# extract_calibration_data (ensemble) ------------------------------------------
setMethod(
  "extract_calibration_data",
  signature(object = "familiarEnsemble"),
  function(
    object,
    data,
    cl = NULL,
    ensemble_method = waiver(),
    evaluation_times = waiver(),
    detail_level = waiver(),
    estimation_type = waiver(),
    aggregate_results = waiver(),
    confidence_level = waiver(),
    bootstrap_ci_method = waiver(),
    is_pre_processed = FALSE,
    message_indent = 0L,
    verbose = FALSE,
    ...
  ) {
    
    # Message extraction start
    logger_message(
      paste0("Assessing model calibration."),
      indent = message_indent,
      verbose = verbose
    )
    
    require_package(
      x = "harmonicmeanp",
      purpose = "to compute p-values for model calibration tests",
      message_type = "warning"
    )
    
    # Load evaluation_times from the object settings attribute, if it is not
    # provided.
    if (is.waive(evaluation_times)) evaluation_times <- object@settings$eval_times
    if (is.waive(ensemble_method)) ensemble_method <- object@settings$ensemble_method
    if (is.waive(confidence_level)) confidence_level <- object@settings$confidence_level
    if (is.waive(bootstrap_ci_method)) bootstrap_ci_method <- object@settings$bootstrap_ci_method
    
    # Test if models are properly loaded
    if (!is_model_loaded(object = object)) ..error_ensemble_models_not_loaded()
    if (!model_is_trained(object = object)) return(NULL)
    
    # Check whether results should be aggregated.
    aggregate_results <- .parse_aggregate_results(
      x = aggregate_results,
      object = object,
      default = TRUE,
      data_element = "calibration_data"
    )
    
    proto_data_element <- .create_extract_calibration_data_object(
      object = object,
      ensemble_method = ensemble_method,
      evaluation_times = evaluation_times,
      detail_level = detail_level,
      estimation_type = estimation_type,
      confidence_level = confidence_level,
      bootstrap_ci_method = bootstrap_ci_method
    )
    
    # Generate elements to send to dispatch.
    calibration_data <- extract_dispatcher(
      FUN = .extract_calibration_data,
      has_internal_bootstrap = TRUE,
      cl = cl,
      object = object,
      data = data,
      proto_data_element = proto_data_element,
      is_pre_processed = is_pre_processed,
      ensemble_method = ensemble_method,
      evaluation_times = evaluation_times,
      aggregate_results = aggregate_results,
      message_indent = message_indent + 1L,
      verbose = verbose
    )
    
    return(calibration_data)
  }
)



# extract_calibration_data (prediction table) ----------------------------------
setMethod(
  "extract_calibration_data",
  signature(object = "familiarDataElementPredictionTable"),
  function(
    object,
    data,
    cl = NULL,
    ensemble_method = waiver(),
    evaluation_times = waiver(),
    detail_level = waiver(),
    estimation_type = waiver(),
    aggregate_results = waiver(),
    confidence_level = waiver(),
    bootstrap_ci_method = waiver(),
    is_pre_processed = FALSE,
    message_indent = 0L,
    verbose = FALSE,
    ...
  ) {

    # Message extraction start
    logger_message(
      paste0("Assessing model calibration."),
      indent = message_indent,
      verbose = verbose
    )
    
    require_package(
      x = "harmonicmeanp",
      purpose = "to compute p-values for model calibration tests",
      message_type = "warning"
    )
    
    if (is.waive(evaluation_times) && methods::.hasSlot(object, "time")) {
      evaluation_times <- object@time
    }
    if (is.waive(ensemble_method)) {
      ensemble_method <- "median"
      if (methods::.hasSlot(object, "ensemble_method")) ensemble_method <- object@ensemble_method
    } 
    
    # Default Values.
    if (is.waive(detail_level)) detail_level <- "ensemble"
    if (is.waive(estimation_type)) estimation_type <- "bootstrap_confidence_interval" 
    if (is.waive(confidence_level)) confidence_level <- 0.95
    if (is.waive(bootstrap_ci_method)) bootstrap_ci_method <- "bc"
    if (is.waive(aggregate_results)) aggregate_results <- TRUE
    
    # Check whether results should be aggregated.
    aggregate_results <- .parse_aggregate_results(
      x = aggregate_results,
      object = object,
      default = TRUE,
      data_element = "calibration_data"
    )
    
    proto_data_element <- .create_extract_calibration_data_object(
      object = object,
      ensemble_method = ensemble_method,
      evaluation_times = evaluation_times,
      detail_level = detail_level,
      estimation_type = estimation_type,
      confidence_level = confidence_level,
      bootstrap_ci_method = bootstrap_ci_method
    )
    
    # Generate elements to send to dispatch.
    calibration_data <- extract_dispatcher(
      FUN = .extract_calibration_data,
      has_internal_bootstrap = TRUE,
      cl = cl,
      object = object,
      data = data,
      proto_data_element = proto_data_element,
      is_pre_processed = is_pre_processed,
      ensemble_method = ensemble_method,
      evaluation_times = evaluation_times,
      aggregate_results = aggregate_results,
      message_indent = message_indent + 1L,
      verbose = verbose
    )
    
    return(calibration_data)
  }
)



.create_extract_calibration_data_object <- function(
    object,
    ensemble_method,
    evaluation_times,
    detail_level,
    estimation_type,
    confidence_level,
    bootstrap_ci_method
) {
    # Check evaluation_times argument
  if (object@outcome_type %in% c("survival")) {
    sapply(
      evaluation_times,
      .check_number_in_valid_range,
      var_name = "evaluation_times",
      range = c(0.0, Inf),
      closed = c(FALSE, TRUE)
    )
  }
  
  # Check ensemble_method argument
  .check_parameter_value_is_valid(
    x = ensemble_method,
    var_name = "ensemble_method",
    values = .get_available_ensemble_prediction_methods()
  )
  
  # Check confidence_level input argument
  .check_number_in_valid_range(
    x = confidence_level,
    var_name = "confidence_level",
    range = c(0.0, 1.0),
    closed = c(FALSE, FALSE)
  )
  
  # Check boostrap method
  .check_parameter_value_is_valid(
    x = bootstrap_ci_method,
    var_name = "bootstrap_ci_method",
    values = .get_available_bootstrap_confidence_interval_methods()
  )
  
  # Check the level detail.
  detail_level <- .parse_detail_level(
    x = detail_level,
    object = object,
    default = "hybrid",
    data_element = "calibration_data"
  )
  
  # Check the estimation type.
  estimation_type <- .parse_estimation_type(
    x = estimation_type,
    object = object,
    default = "bootstrap_confidence_interval",
    data_element = "calibration_data",
    detail_level = detail_level,
    has_internal_bootstrap = TRUE
  )
  
  # Generate a prototype data element.
  proto_data_element <- new(
    "familiarDataElementCalibrationData",
    detail_level = detail_level,
    estimation_type = estimation_type,
    confidence_level = confidence_level,
    bootstrap_ci_method = bootstrap_ci_method
  )
  
  return(proto_data_element)
}



.extract_calibration_data <- function(
    object,
    proto_data_element,
    evaluation_times = NULL,
    aggregate_results,
    cl,
    ...
) {
  # Ensure that the object is loaded
  object <- load_familiar_object(object)
  
  # Add model name.
  proto_data_element <- add_model_name(proto_data_element, object = object)
  
  # Add evaluation time as a identifier to the data element.
  if (length(evaluation_times) > 0L && object@outcome_type == "survival") {
    data_elements <- add_data_element_identifier(
      x = proto_data_element,
      evaluation_time = evaluation_times
    )
    
  } else {
    data_elements <- list(proto_data_element)
  }
  
  # Iterate over data elements.
  calibration_data <- mapply(
    ..extract_calibration_data,
    data_element = data_elements,
    MoreArgs = list(
      object = object,
      aggregate_results = aggregate_results,
      cl = cl,
      ...
    ),
    SIMPLIFY = FALSE
  )
  
  return(calibration_data)
}



# ..extract_calibration_data (generic) -----------------------------------------
setGeneric(
  "..extract_calibration_data",
  function(object, ...) standardGeneric("..extract_calibration_data")
)



# ..extract_calibration_data (character) ---------------------------------------
setMethod(
  "..extract_calibration_data",
  signature(object = "character"),
  function(
    object,
    ...
  ){
    # Ensure that the object is loaded
    object <- load_familiar_object(object)
    
    return(..extract_calibration_data(object = object, ...))
  }
)



# ..extract_calibration_data (model, ensemble) ---------------------------------
setMethod(
  "..extract_calibration_data",
  signature(object = "familiarModelUnion"),
  function(
    object,
    data,
    data_element,
    is_pre_processed,
    ensemble_method,
    ...
  ) {
    # Aggregate data.
    data <- aggregate_data(data)
    
    if (object@outcome_type %in% c("survival", "competing_risk")) {
      # Predict class probabilities or regression values.
      prediction_data <- .predict(
        object = object,
        data = data,
        type = "survival_probability",
        time = data_element@identifiers$evaluation_time,
        ensemble_method = ensemble_method,
        is_pre_processed = is_pre_processed
      )
      
    } else {
      # Predict class probabilities or regression values.
      prediction_data <- .predict(
        object = object,
        data = data,
        ensemble_method = ensemble_method,
        is_pre_processed = is_pre_processed
      )
    }
    
    return(..extract_calibration_data(
      object = prediction_data,
      data_element = data_element,
      ...
    ))
  }
)



# ..extract_calibration_data (prediction table) --------------------------------
setMethod(
  "..extract_calibration_data",
  signature(object = "familiarDataElementPredictionTable"),
  function(
    object,
    data_element,
    aggregate_results,
    cl = NULL,
    progress_bar = FALSE,
    verbose = FALSE,
    message_indent,
    ...
  ) {
    if (!(
      is(object, "predictionTableSurvivalProbability") || 
      is(object, "predictionTableClassification") ||
      is(object, "predictionTableRegression")
    )) {
      ..warning_no_data_extraction_from_prediction_table("calibration data")
      
      return(NULL)
    }
    
    # Message the user concerning the time at which metrics are computed. This
    # is only relevant for survival analysis.
    if (length(data_element@identifiers$evaluation_time) > 0L && progress_bar) {
      logger_message(
        paste0(
          "Assessing model calibration at time ",
          data_element@identifiers$evaluation_time, "."
        ),
        indent = message_indent,
        verbose = verbose
      )
    }
    
    
    # Check if any predictions are valid.
    if (!all_predictions_valid(object)) return(NULL)
    
    # Remove data with missing outcomes.
    object <- filter_missing_outcome(object)
    
    # Check that any prediction data remain.
    if (is_empty(object)) return(NULL)
    
    # Add positive class as an identifier.
    if (object@outcome_type %in% c("binomial")) {
      data_element <- add_data_element_identifier(
        x = data_element,
        positive_class = get_outcome_class_levels(object)[2L]
      )
      
    } else if (object@outcome_type %in% c("multinomial")) {
      data_element <- add_data_element_identifier(
        x = data_element,
        positive_class =  get_outcome_class_levels(object)
      )
    }

    # Add bootstrap data.
    bootstrap_data <- add_data_element_bootstrap(
      x = data_element,
      ...
    )
    
    # Add distribution data.
    if (!is.list(data_element)) data_element <- list(data_element)
    
    density_data <- lapply(
      data_element,
      .compute_calibration_data_density,
      object = object
    )
    
    # Iterate over elements.
    data_elements <- fam_mapply(
      cl = cl,
      assign = NULL,
      FUN = .compute_calibration_data,
      data_element = bootstrap_data$data_element,
      bootstrap = bootstrap_data$bootstrap,
      bootstrap_seed = bootstrap_data$seed,
      MoreArgs = list("object" = object),
      progress_bar = progress_bar && verbose,
      chopchop = TRUE
    )
    
    # Flatten list of data elements.
    data_elements <- unlist(data_elements)
    if (!is.list(data_elements)) data_elements <- list(data_elements)

    # Add in density data elements.
    data_elements <- c(data_elements, density_data)
    
    # Merge data elements
    data_elements <- merge_data_elements(data_elements)
    
    if (aggregate_results) data_elements <- .compute_data_element_estimates(x = data_elements)
    
    return(data_elements)
  }
)



# ..extract_calibration_data (NULL) --------------------------------------------
setMethod(
  "..extract_calibration_data",
  signature(object = "NULL"),
  function(object, ...) {
    return(NULL)
  }
)



.compute_calibration_data <- function(
    data_element,
    object,
    bootstrap,
    bootstrap_seed
) {
  
  # Suppress NOTES due to non-standard evaluation in data.table
  observed <- NULL
  
  # Start random number generator.
  rstream_object <- .start_random_number_stream(seed = ifelse(is.finite(bootstrap_seed), bootstrap_seed, 19L))
  
  # Bootstrap the data.
  if (bootstrap) {
    object <- get_bootstrap_sample(data = object, rstream_object = rstream_object)
  }
  
  # Compute calibration, linear and gof test data.
  calibration_data <- ..compute_calibration_data(
    object,
    data_element = data_element,
    rstream_object = rstream_object
  )
  
  # Extract data.
  calibration_at_large <- calibration_data$linear
  calibration_gof_test <- calibration_data$gof_test
  calibration_data <- calibration_data$calibration

  if (!is_empty(calibration_data) && data_element@estimation_type != "point") {
    calibration_data <- ..compute_calibration_data_interpolated(
      data = calibration_data
    )
    
    # Keep only finite predictions,
    calibration_data <- calibration_data[is.finite(observed)]
  }
  
  # Create additional data elements to contain linear fit data and GoF test
  # data.
  linear_fit_data_element <- methods::new(
    "familiarDataElementCalibrationLinearFit",
    data_element
  )  
  gof_data_element <- methods::new(
    "familiarDataElementCalibrationGoodnessOfFit",
    data_element
  )
  
  # Set for calibration data.
  data_element@data <- calibration_data
  data_element@value_column <- "observed"
  data_element@grouping_column <- "expected"
  
  # Set linear fit data, i.e. calibration-in-the-large and calibration slope.
  linear_fit_data_element@data <- calibration_at_large
  linear_fit_data_element@value_column <- "value"
  linear_fit_data_element@grouping_column <- "type"
  
  # Set goodness of fit data.
  gof_data_element@data <- calibration_gof_test
  gof_data_element@value_column <- "p_value"
  gof_data_element@grouping_column <- "type"
  
  return(list(
    data_element,
    linear_fit_data_element,
    gof_data_element
  ))
}



# ..compute_calibration_data (generic) -----------------------------------------
setGeneric(
  "..compute_calibration_data",
  function(object, ...) standardGeneric("..compute_calibration_data")
)




# ...compute_calibration_data (generic) ----------------------------------------
setGeneric(
  "...compute_calibration_data",
  function(object, ...) standardGeneric("...compute_calibration_data")
)



# ..compute_calibration_data (survival probability) ----------------------------
setMethod(
  "..compute_calibration_data",
  signature(object = "predictionTableSurvivalProbability"),
  function(object, data_element, rstream_object = NULL) {
    # Generate baseline calibration data for survival models from an input table
    # (data) containing survival probabilities for each sample and
    # corresponding time and event status.
    
    # Suppress NOTES due to non-standard evaluation in data.table
    exp_prob <- NULL
    
    # Set time.
    time <- data_element@identifiers$evaluation_time

    # Find sample identifiers.
    sample_identifiers <- get_id_columns(id_depth = "series")
    
    # Rename the survival column to standard name.
    data <- data.table::copy(.as_data_table(object))
    data.table::setnames(
      data,
      old = c("outcome_time", "outcome_event", "predicted_outcome"),
      new = c("time", "event", "exp_prob")
    )
    
    # Sort by survival probability.
    data <- data[order(exp_prob)]
    
    # Randomly split into groups. The number of groups is determined using
    # sturges rule.
    group_data <- .calibration_create_randomised_groups(
      x = data$exp_prob,
      y = data$event,
      sample_identifiers = data[, mget(sample_identifiers)],
      n_min_y_in_group = 4L,
      rstream_object = rstream_object
    )
    
    if (is_empty(group_data)) return(NULL)
    
    # Merge with prediction table.
    group_data <- merge(
      x = data,
      y = group_data,
      by = sample_identifiers,
      allow.cartesian = TRUE
    )
    
    # Iterate over groups and add details by comparing the kaplan-meier survival
    # curve within each group at time with the mean survival probability in
    # the group.
    calibration_data <- lapply(
      split(group_data, by = "group_id"),
      function(x, object, time) {
        ...compute_calibration_data(
          object = object,
          data = x,
          time = time
        )
      },
      object = object,
      time = time
    )
    
    # Combine to a single list.
    calibration_data <- data.table::rbindlist(
      calibration_data,
      use.names = TRUE
    )
    
    # Prevent processing of empty tables.
    if (is_empty(calibration_data)) return(NULL)
    
    # Set column order.
    data.table::setcolorder(
      x = calibration_data,
      neworder = c("expected", "observed", "km_var", "n_g")
    )
    
    # Calibration-in-the-large and calibration slope
    calibration_at_large <- .compute_calibration_linear_fit(
      calibration_data = calibration_data,
      outcome_type = object@outcome_type
    )
    
    # Nam-D'Agostino tests
    calibration_gof_test <- .compute_calibration_nam_dagostino(
      calibration_data = calibration_data
    )
    
    return(list(
      "calibration" = calibration_data,
      "linear" = calibration_at_large,
      "gof_test" = calibration_gof_test
    ))
  }
)



# ..compute_calibration_data (classification) ----------------------------------
setMethod(
  "..compute_calibration_data",
  signature(object = "predictionTableClassification"),
  function(object, data_element, rstream_object = NULL) {
    # For assessing the calibration of categorical outcomes, we require expected
    # and observed probabilities for 1 (binomial) or all classes (multinomial).
    # Expected probabilities are easy, as the model predicts them. Observed
    # probabilities are however based on class proportion within a group. This
    # means that we need to define groups, ordered by predicted probability. The
    # same groups need to be defined for the Hosmer-Lemeshow test, so we might
    # as well obtain all the data here.
    #
    # Groups are defined for each class level because the samples should be
    # ordered according to predicted class probability to create the groups.
    
    # Suppress NOTES due to non-standard evaluation in data.table
    observed_class <- exp_prob <- NULL
    
    positive_class <- data_element@identifiers$positive_class

    # Identify the real outcome columns
    outcome_column <- get_outcome_columns(x = "binomial")
    
    # Find sample identifiers.
    sample_identifiers <- get_id_columns(id_depth = "series")
    
    # Make local copy.
    data <- data.table::copy(
      .as_data_table(object)[, mget(c(sample_identifiers, outcome_column, positive_class))]
    )
    
    # Rename columns to standard names
    data.table::setnames(
      data,
      old = c(outcome_column, positive_class),
      new = c("observed_class", "exp_prob")
    )
    
    # Mask class so that positive outcomes are TRUE and the rest FALSE
    data[, "observed_class" := observed_class == positive_class]
    
    # Sort by probability
    data <- data[order(exp_prob)]
    
    # Randomly split into groups. The number of groups is determined using
    # sturges rule.
    group_data <- .calibration_create_randomised_groups(
      x = data$exp_prob,
      sample_identifiers = data[, mget(sample_identifiers)],
      rstream_object = rstream_object
    )
    
    if (is_empty(group_data)) return(NULL)
    
    # Merge with prediction table.
    group_data <- merge(
      x = data,
      y = group_data,
      by = sample_identifiers,
      allow.cartesian = TRUE
    )
    
    # Compute calibration data from each group.
    calibration_data <- lapply(
      split(group_data, by = "group_id"),
      function(x, object) {
        ...compute_calibration_data(
          object = object,
          data = x
        )
      },
      object = object
    )
    
    # Combine to a single list.
    calibration_data <- data.table::rbindlist(
      calibration_data,
      use.names = TRUE
    )
    
    # Check that any calibration data was generated.
    if (is_empty(calibration_data)) return(NULL)
    
    # Set column order
    data.table::setcolorder(
      x = calibration_data,
      neworder = c("expected", "observed", "n_g", "n_pos", "n_neg")
    )
    
    # Calibration-in-the-large and calibration slope
    calibration_at_large <- .compute_calibration_linear_fit(
      calibration_data = calibration_data,
      outcome_type = object@outcome_type
    )
    
    # Hosmer-Lemeshow tests
    calibration_gof_test <- .compute_calibration_hosmer_lemeshow(
      calibration_data = calibration_data
    )
    
    return(list(
      "calibration" = calibration_data,
      "linear" = calibration_at_large,
      "gof_test" = calibration_gof_test
    ))
  }
)


# ..compute_calibration_data (regression) --------------------------------------
setMethod(
  "..compute_calibration_data",
  signature(object = "predictionTableRegression"),
  function(object, data_element, rstream_object = NULL) {
    # Calibration for regression problems is pretty straightforward. However, for
    # goodness-of-fit tests, we need to constrain expected and observed value
    # ranges to [0, 1]. To do so, we use the range of outcome values from the
    # development data, as stored in the calibration_info slot of the input
    # object.
    
    # Suppress NOTES due to non-standard evaluation in data.table
    outcome <- predicted_outcome <- NULL
    n_groups <- 1L
    
    # Extract data.
    data <- data.table::copy(.as_data_table(object))
    
    # Find sample identifiers.
    sample_identifiers <- get_id_columns(id_depth = "series")
    
    # Determine the outcome range
    outcome_range <- object@observed_value_range
    if (anyNA(outcome_range)) outcome_range <- range(data$outcome, na.rm = TRUE, finite = TRUE)
    
    # Get the shift and scale parameters.
    norm_shift <- outcome_range[1L]
    norm_scale <- diff(outcome_range)
    
    # Set scale parameters equal to 0.0 to 1.0 to avoid division by 0.
    if (norm_scale == 0.0) norm_scale <- 1.0
    
    # Apply normalisation.
    data[, ":="(
      "expected" = (predicted_outcome - norm_shift) / norm_scale,
      "observed" = (outcome - norm_shift) / norm_scale
    )]
    
    # Randomly split into groups. The number of groups is determined using
    # sturges rule.
    group_data <- .calibration_create_randomised_groups(
      x = data$expected,
      sample_identifiers = data[, mget(sample_identifiers)],
      rstream_object = rstream_object
    )
    
    if (is_empty(group_data)) return(NULL)
    
    # Merge with prediction table.
    group_data <- merge(
      x = data,
      y = group_data,
      by = sample_identifiers,
      allow.cartesian = TRUE
    )
    
    # Compute calibration data from each group.
    calibration_data <- lapply(
      split(group_data, by = "group_id"),
      function(x, object) {
        ...compute_calibration_data(
          object = object,
          data = x
        )
      },
      object = object
    )
    
    # Combine to a single list.
    calibration_data <- data.table::rbindlist(
      calibration_data,
      use.names = TRUE
    )
    
    # Check if the resulting table contains data.
    if (is_empty(calibration_data)) return(NULL)
    
    # Make columns match the expected order
    data.table::setcolorder(
      x = calibration_data,
      neworder = c("expected", "observed", "n_g")
    )
    
    # Calibration-in-the-large and calibration slope
    calibration_at_large <- .compute_calibration_linear_fit(
      calibration_data = calibration_data,
      outcome_type = object@outcome_type
    )
    
    # Hosmer-Lemeshow tests
    calibration_gof_test <- .compute_calibration_hosmer_lemeshow(
      calibration_data = calibration_data
    )
    
    return(list(
      "calibration" = calibration_data,
      "linear" = calibration_at_large,
      "gof_test" = calibration_gof_test
    ))
  }
)



# ...compute_calibration_data (survival probability) ---------------------------
setMethod(
  "...compute_calibration_data",
  signature(object = "predictionTableSurvivalProbability"),
  function(
    object,
    data,
    time
  ) {
    # Set default values.
    obs_prob <- exp_prob <- km_var <- NA_real_
    n_g <- NA_integer_
    
    if (nrow(data) >= 2L) {
      
      # Fit a Kaplan-Meier curve for the current group
      km_fit <- survival::survfit(Surv(time, event) ~ 1, data = data)
      
      if (length(km_fit$time) >= 2L) {
        
        # Get observed probability
        obs_prob <- stats::approx(
          x = km_fit$time,
          y = km_fit$surv,
          xout = time,
          yleft = 1.0,
          method = "constant",
          rule = 2L
        )$y
        
        # Get expected probability
        exp_prob <- mean(data$exp_prob)
        
        # Get group size
        n_g <- nrow(data)
        
        # Get greenwood variance estimate.
        km_var <- stats::approx(
          x = km_fit$time,
          y = km_fit$std.err,
          xout = time,
          method = "constant",
          rule = 2L
        )$y^2.0
        
      } 
    }
    
    # Create table.
    calibration_table <- list(
      "expected" = exp_prob,
      "observed" = obs_prob,
      "n_g" = n_g,
      "km_var" = km_var
    )
    
    return(calibration_table)
  }
)



# ...compute_calibration_data (classification) ---------------------------------
setMethod(
  "...compute_calibration_data",
  signature(object = "predictionTableClassification"),
  function(
    object,
    data
  ) {
    # Check that the groups list contains at least one entry.
    if (is_empty(data)) return(NULL)
  
    # Create table
    calibration_table <- list(
      # Mean expected probability in a group.
      "expected" = mean(data$exp_prob),
      # Observed proportion of positive class in a group.
      "observed" = mean(data$observed_class),
      # Number of samples in the group
      "n_g" = nrow(data),
      # Number of samples with the positive class in each group.
      "n_pos" = sum(data$observed_class),
      # Number of samples with the negative class in each group.
      "n_neg" = sum(!data$observed_class)
    )
    
    return(calibration_table)
  }
)
    


# ...compute_calibration_data (regression) ---------------------------------
setMethod(
  "...compute_calibration_data",
  signature(object = "predictionTableRegression"),
  function(
    object,
    data
  ) {
    # Check that the groups list contains at least one entry.
    if (is_empty(data)) return(NULL)
    
    # Create table
    calibration_table <- list(
      # Mean expected value in the group.
      "expected" = mean(data$expected),
      # Observed value in the group.
      "observed" = mean(data$observed),
      # Number of samples in the group
      "n_g" = nrow(data)
    )
    
    return(calibration_table)
  }
)



#' Create randomised groups Creates randomised groups, e.g. for tests that
#' depend on splitting (continuous) data into groups, such as the
#' Hosmer-Lemeshow test
#'
#' The default fast mode is based on random sampling, whereas the slow mode is
#' based on probabilistic joining of adjacent groups. As the name suggests, fast
#' mode operates considerably more efficient.
#'
#' @param x Vector with data used for sorting. Groups are formed based on
#'   adjacent values.
#' @param y Vector with markers, e.g. the events. Should be 0 or 1 (for an
#'   event).
#' @param sample_identifiers data.table with sample_identifiers. If provide, a
#'   list of grouped sample_identifiers will be returned, and integers
#'   otherwise.
#' @param n_max_groups Maximum number of groups that need to be formed.
#' @param n_min_groups Minimum number of groups that need to be formed.
#' @param n_min_y_in_group Minimum number of y=1 in each group for a valid
#'   group.
#'
#' @details Creates randomised groups, e.g. for tests that depend on splitting
#'   (continuous) data into groups, such as the Hosmer-Lemeshow test
#'
#' @return data.table with sample identifiers and random group ids.
#' @md
#' @keywords internal
.calibration_create_randomised_groups <- function(
    x, 
    y = NULL, 
    sample_identifiers, 
    n_max_groups = NULL, 
    n_min_groups = NULL, 
    n_min_y_in_group = NULL, 
    unique_samples_only = TRUE,
    rstream_object = NULL
) {
  # Suppress NOTES due to non-standard evaluation in data.table
  group_id <- cum_y <- weight <- cum_prob_lower <- cum_prob_upper <- NULL
  exclude <- y_in_group <- NULL
  
  # Populate the generic table.
  data <- data.table::copy(sample_identifiers[, mget(get_id_columns(id_depth = "series"))])
  data[, "x" := x]
  if (is.null(y)) {
    data[, "y" := 0L]
  } else {
    data[, "y" := y]
  }
  if (unique_samples_only) data <- unique(data)
  
  # - Get number of x
  n_x <- nrow(data)
  
  # - Get number of y
  n_y <- sum(data$y)
  
  # - Get the maximum of groups
  if (is.null(n_max_groups)) {
    n_max_groups <- ceiling(2.5 * n_x^(1.0 / 3.0))
  }
  
  # - Update maximum number of groups
  if (n_y > 0.0 && !is.null(n_min_y_in_group)) {
    n_max_groups <- min(c(n_max_groups, floor(n_y / n_min_y_in_group)))
  }
  
  # - Update mininum number of groups based on n_max_groups
  if (is.null(n_min_groups)) {
    n_min_groups <- min(c(n_max_groups, ceiling(1.0 * n_x^(1.0 / 3.0))))
  }
  
  # Some checks
  if (n_max_groups < n_min_groups || n_max_groups > n_x || n_max_groups < 2L) {
    return(NULL)
  }
  data <- data[order(x)]
  
  # Initial loop counter
  loop_iter <- 0L
  
  while (TRUE) {
    # Draw a random number of groups between n_min_groups and n_max_groups
    if (n_min_groups == n_max_groups) {
      n_group_draw <- n_max_groups
    } else {
      n_group_draw <- fam_sample(
        x = seq.int(from = n_min_groups, to = n_max_groups, by = 1L),
        size = 1L, 
        replace = FALSE,
        rstream_object = rstream_object
      )
    }
    
    # Draw a randomised groups assignment
    random_group_id <- sort(
      fam_sample(
        x = seq_len(n_group_draw),
        size = n_x, 
        replace = TRUE,
        rstream_object = rstream_object
      )
    )
    
    # Assign group id
    data[, "group_id" := random_group_id]
    
    # Check if all groups have the minimum required y (e.g. events/group size)
    if (n_y > 0L && !is.null(n_min_y_in_group)) {
      dt_y <- data[, list(y_in_group = sum(y)), by = group_id]
      if (all(dt_y$y_in_group >= n_min_y_in_group)) break
      
    } else {
      # If there is no minimum number of events required, only run one
      # iteration
      break
    }
    
    # Update loop_iter in case the sample was not good enough
    loop_iter <- loop_iter + 1L
    
    # Break from outer loop after 10 unsuccessful random samplings and create
    # a static split instead
    if (loop_iter >= 10L) {
      # Cumulative sum over y
      data[, "cum_y" := cumsum(y)]
      
      # Assign static group ids
      data[, "group_id" := findInterval(
        x = cum_y,
        vec = c(
          -Inf,
          stats::quantile(x = cum_y, (1L:(n_min_groups - 1L)) / n_min_groups),
          Inf
        )
      )]
      
      # Break from loop
      break
    }
  }
  
  # Get sample identifiers for each group
  return(data[, mget(c(get_id_columns(id_depth = "series"), "group_id"))])
}



.compute_calibration_data_density <- function(
    data_element,
    object
) {
  
  # Suppress NOTES due to non-standard evaluation in data.table
  predicted_outcome <- NULL
  
  # Copy data element.
  data_element <- methods::new(
    "familiarDataElementCalibrationDensity",
    data_element
  )
  
  # Copy data
  data <- data.table::copy(.as_data_table(object))
  
  if (object@outcome_type %in% c("survival", "competing_risk")) {
    # Rename the survival column to standard name.
    data.table::setnames(
      data,
      old = "predicted_outcome",
      new = "expected"
    )
    
  } else if (object@outcome_type %in% c("binomial", "multinomial")) {
    # Rename class probability columns to standard names.
    data.table::setnames(
      data,
      old = data_element@identifiers$positive_class,
      new = "expected"
    )
    
  } else if (object@outcome_type %in% c("continuous")) {
    # Determine the outcome range
    outcome_range <- range(object@observed_value_range)
    
    # Get the shift and scale parameters.
    norm_shift <- outcome_range[1L]
    norm_scale <- diff(outcome_range)
    
    # Set scale parameters equal to 0.0 to 1.0 to avoid division by 0.
    if (norm_scale == 0.0) norm_scale <- 1.0
    
    # Apply normalisation.
    data[, ":="("expected" = (predicted_outcome - norm_shift) / norm_scale)]
    
  } else {
    ..error_outcome_type_not_implemented(object@outcome_type)
  }
  
  # Get expected values.
  expected <- data$expected
  
  # Set interpolation range for expected values. This allows for aggregating
  # observed values based on the expected value.
  expected_interpolated <- seq(
    from = 0.000,
    to = 1.000,
    by = 0.005
  )
  
  # Determine the bin for expected.
  n <- length(expected_interpolated)
  
  # Discretise bins.
  group_bin <- as.integer(floor((n - 1L) * expected)) + 1L
  
  # Check if the group_bin still falls within the valid range.
  group_bin <- ifelse(group_bin > n, n, group_bin)
  group_bin <- ifelse(group_bin < 1L, 1L, group_bin)
  
  # Compute the frequency of each entry.
  data <- data.table::data.table("expected" = expected_interpolated[group_bin])
  data <- data[, list("frequency" = .N / length(group_bin)), by = "expected"]
  
  # Set linear fit data, i.e. calibration-in-the-large and calibration slope.
  data_element@data <- data
  data_element@value_column <- "frequency"
  data_element@grouping_column <- "expected"
  
  return(data_element)
}



.compute_calibration_linear_fit <- function(
    calibration_data, 
    outcome_type
) {
  # Tests calibration-at-the-large and calibration slope
  
  # Suppress NOTES due to non-standard evaluation in data.table.
  expected <- observed <- p_value <- NULL
  
  # Check that calibration data are not empty.
  if (is_empty(calibration_data)) return(NULL)
  
  # Filter out non-finite expected and observed values.
  calibration_data <- calibration_data[is.finite(expected) & is.finite(observed)]
  
  # Check if the input is empty or only contains one entry
  if (nrow(calibration_data) < 2L) return(NULL)
  
  # Fit per repeated measurement.
  fit_data <- ..compute_calibration_fit(
    calibration_data,
    outcome_type = outcome_type
  )
  
  # Remove NaNs
  fit_data <- fit_data[is.finite(p_value)]
  if (is_empty(fit_data)) return(NULL)
  
  return(fit_data)
}



..compute_calibration_fit <- function(
    calibration_data,
    outcome_type
) {
  
  # Suppress NOTES due to non-standard evaluation in data.table.
  type <- value <- ci_low <- ci_up <- NULL
  
  # Perform a linear fit. Note that the slope is offset by 1*expected to
  # directly compare with the expected slope of 1.
  fit <- stats::lm(
    observed ~ expected + 1 + offset(expected),
    data = calibration_data
  )
  
  # Get coefficients for intercept and slope
  fit_coef <- stats::coef(fit)
  
  # Determine if there is no slope.
  no_slope <- is.na(fit_coef["expected"]) && data.table::uniqueN(calibration_data$expected) == 1L
  if (no_slope) fit_coef["expected"] <- 0.0
  
  # Get confidence intervals
  fit_conf_int <- suppressWarnings(stats::confint(fit))
  if (no_slope) fit_conf_int["expected", ] <- c(-Inf, Inf)
  
  # Get summary
  fit_summary <- suppressWarnings(stats::summary.lm(fit))$coefficients
  if (no_slope) {
    fit_summary <- rbind(
      fit_summary,
      matrix(data = c(0.0, Inf, 0.0, 1.0), ncol = 4L)
    )
    rownames(fit_summary) <- c("(Intercept)", "expected")
  }
  
  # Correct for having multiple overlapping groups in binomial and multinomial
  # tests.
  if (outcome_type %in% c("binomial", "multinomial", "survival")) {
    
    # Determine the number of groups.
    n_groups <- 1L
    
    # Recompute the standard deviation
    fit_summary[, 2L] <- fit_summary[, 2L] * sqrt(n_groups)
    
    # Recompute the t-score
    fit_summary[, 3L] <- fit_summary[, 1L] / fit_summary[, 2L]
    fit_summary[, 3L][fit_summary[, 2L] == 0.0] <- Inf
    
    # Recompute the p-value
    fit_summary[, 4L] <- 2.0 * stats::pt(abs(fit_summary[, 3L]), fit$df.residual, lower.tail = FALSE)
    
    # Adapt the confidence intervals
    fit_conf_int[, 1L] <- fit_summary[, 1L] + fit_summary[, 2L] * stats::qnorm(0.025)
    fit_conf_int[, 2L] <- fit_summary[, 1L] + fit_summary[, 2L] * stats::qnorm(0.975)
  }
  
  # Construct data table and update the slope so that the expected slope
  # (slope+1) is shown instead.
  calibration_at_large <- data.table::data.table(
    "type" = c("offset", "slope"),
    "value" = fit_coef + c(0.0, 1.0),
    "ci_low" = fit_conf_int[, 1L] + c(0.0, 1.0),
    "ci_up" = fit_conf_int[, 2L] + c(0.0, 1.0),
    "p_value" = fit_summary[, 4L]
  )
  
  return(calibration_at_large)
}



.compute_calibration_hosmer_lemeshow <- function(calibration_data) {
  
  # Suppress NOTES due to non-standard evaluation in data.table
  observed <- expected <- n_g <- hm_group <- statistic <- n_groups <- NULL
  
  # Check for empty calibration tables
  if (is_empty(calibration_data)) return(NULL)
  
  # Make local copy of the calibration data.
  calibration_data <- data.table::copy(calibration_data)
  
  # Remove NA predictions.
  calibration_data <- calibration_data[is.finite(expected)]
  if (is_empty(calibration_data)) return(NULL)
  
  # Compute test statistic for each group
  calibration_data[
    ,
    "hm_group" := ..compute_calibration_test_statistic(observed, expected, n_g),
    by = seq_len(nrow(calibration_data))
  ]
  
  # Remove NA values in hm_group
  calibration_data <- calibration_data[is.finite(hm_group)]
  if (is_empty(calibration_data)) return(NULL)
  
  # Compute test statistic for each time point
  gof_table <- calibration_data[
    ,
    list("statistic" = sum(hm_group, na.rm = TRUE), n_groups = .N)
  ]
  
  # Remove all entries with n_groups < 3
  gof_table <- gof_table[n_groups >= 3L]
  
  if (is_empty(gof_table)) return(NULL)
  
  # Perform chi-square test
  gof_table <- gof_table[
    ,
    list("p_value" = stats::pchisq(
      q = statistic,
      df = n_groups - 2L,
      lower.tail = FALSE
    ))
  ]
  
  # Add type.
  gof_table[, "type" := "hosmer_lemeshow"]

  # Reorder columns
  data.table::setcolorder(gof_table, c("type", "p_value"))
  
  return(gof_table)
}



..compute_calibration_test_statistic <- function(
    observed,
    expected,
    n_g,
    method = "approximate"
) {
  
  if (!method %in% c("exact", "approximate")) {
    ..error_reached_unreachable_code(
      "..compute_hm_test_statistic: method should be exact or approximate"
    )
  }
  
  # Check for infinite or NA values.
  if (!is.finite(expected) || !is.finite(observed)) return(NA_real_)
  
  if (method == "exact") {
    # Use exact computation for comparison. Note that this has asymptotic
    # behaviour for expected -> 0 and expected -> 1.
    value <- (observed - expected)^2.0 * n_g / (expected * (1.0 - expected))
    
    return(value)
    
  } else if (method == "approximate") {
    # Use exact computation for comparison.
    exact_value <- (observed - expected)^2.0 * n_g / (expected * (1.0 - expected))
    if (!is.finite(exact_value)) exact_value <- Inf
    
    if (expected <= 0.05) {
      
      # Use Taylor series for 1/(x * (1-x)) to second order at 0.05.
      value <- n_g * (expected - observed)^2.0 * 400.0 *
        (1.0 / 19.0 - 360.0 / 361.0 * (expected - 1.0 / 20.0))
      
      if (exact_value < value) value <- exact_value
      if (value < 0.0) value <- 0.0
      
    } else if (expected >= 0.95) {
      
      # Use Taylor series for 1/(x * (1-x)) to second order at 0.95.
      value <- n_g * (expected - observed)^2.0 * 400.0 *
        (1.0 / 19.0 + 360.0 / 361.0 * (expected - 19.0 / 20.0))
      
      if (exact_value < value) value <- exact_value
      if (value < 0.0) value <- 0.0
      
    } else {
      # Use exact computation in between.
      value <- exact_value
    }
    
    return(value)
  }
}



.compute_calibration_nam_dagostino <- function(calibration_data) {
  # Nam-D'Agostino and Greenwood-Nam-D'Agostino tests. See Demler et al. 2015.
  
  # Suppress NOTES due to non-standard evaluation in data.table
  observed <- expected <- n_g <- km_var <- nd_group <- gnd_group <- n_groups <- NULL
  nam_dagostino <- greenwood_nam_dagostino <- NULL
  
    # Check for empty calibration tables
  if (is_empty(calibration_data)) return(NULL)
  
  # Make local copy of calibration_data to avoid updating the main copy by
  # reference.
  calibration_data <- data.table::copy(calibration_data)
  
  # Remove NA predictions.
  calibration_data <- calibration_data[is.finite(expected)]
  if (is_empty(calibration_data)) return(NULL)
  
  # Compute test statistic for each group.
  calibration_data[
    ,
    ":="(
      "nd_group" = ..compute_calibration_test_statistic(observed, expected, n_g),
      "gnd_group" = (observed - expected)^2.0 / km_var
    ),
    by = seq_len(nrow(calibration_data))
  ]
  
  # Remove NA values in hm_group
  calibration_data <- calibration_data[is.finite(nd_group) & is.finite(gnd_group)]
  if (is_empty(calibration_data)) return(NULL)
  
  # Compute test statistic for each time point
  gof_table <- calibration_data[
    ,
    list(
      "nam_dagostino" = sum(nd_group, na.rm = TRUE),
      "greenwood_nam_dagostino" = sum(gnd_group, na.rm = TRUE),
      "n_groups" = .N
    )
  ]
  
  # Compute test score for both test tests.
  gof_table <- gof_table[
    n_groups > 1L,
    list(
      "nam_dagostino" = stats::pchisq(
        q = nam_dagostino,
        df = n_groups - 1L,
        lower.tail = FALSE
      ),
      "greenwood_nam_dagostino" = stats::pchisq(
        q = greenwood_nam_dagostino,
        df = n_groups - 1L,
        lower.tail = FALSE
      )
    )
  ]
  
  # Check for an empty goodness-of-fit table.
  if (is_empty(gof_table)) return(NULL)
  
  # Melt so that each test is on a separate row
  gof_table <- data.table::melt(
    gof_table,
    variable.name = "type",
    value.name = "p_value",
    measure.vars = colnames(gof_table),
    variable.factor = FALSE
  )
  
  # Reorder columns
  data.table::setcolorder(
    x = gof_table,
    neworder = c("type", "p_value")
  )
  
  return(gof_table)
}



..compute_calibration_data_interpolated <- function(
    data,
    method = "loess"
) {
  
  # Suppress NOTES due to non-standard evaluation in data.table
  observed <- expected <- NULL
  
  # Set interpolation range for expected values. This allows for aggregating
  # observed values based on the expected value.
  expected_interpolated <- seq(
    from = 0.000,
    to = 1.000,
    by = 0.005
  )
  
  # Select only unique entries.
  data <- unique(data)
  
  # Select only instances where expected and observed are finite values.
  data <- data[is.finite(expected) & is.finite(observed)]
  
  # Determine the number of unique expected values.
  n_instances <- data.table::uniqueN(data, by = "expected")
  
  # Switch interpolation method for small number of predictions.
  if (method == "loess") {
    if (n_instances >= 4L) {
      degree <- 2L
      
    } else if (n_instances == 3L) {
      degree <- 1L
      
    } else if (n_instances == 2L) {
      method <- "linear"
      
    } else {
      return(list(
        "expected" = NA_real_,
        "observed" = NA_real_
      ))
    }
    
  } else if (method == "linear") {
    
    if (n_instances < 2L) {
      return(list(
        "expected" = NA_real_,
        "observed" = NA_real_
      ))
    }
    
  } else {
    ..error_reached_unreachable_code(paste0(
      "..compute_calibration_data_interpolated: unexpected interpolation ",
      "method encountered: ", method
    ))
  }
  
  if (method == "linear") {
    # Interpolate observed values. Return NA values outside the known expected
    # range.
    observed_interpolated <- suppressWarnings(stats::approx(
      x = data$expected,
      y = data$observed,
      xout = expected_interpolated,
      method = "linear",
      rule = 1L
    )$y)
    
  } else if (method == "loess") {
    # Set up the loess model.
    loess_predictor <- suppressWarnings(stats::loess(
      observed ~ expected,
      data = data,
      span = 1.0,
      degree = degree,
      normalize = FALSE,
      control = stats::loess.control(surface = "direct")
    ))
    
    # Predicted interpolated values.
    observed_interpolated <- suppressWarnings(predict(
      loess_predictor,
      newdata = expected_interpolated
    ))
  }
  
  # Set expected values and the corresponding interpolated observed values.
  return(data.table::data.table(
    "expected" = expected_interpolated,
    "observed" = observed_interpolated
  ))
}



# ..compute_data_element_estimates (familiarDataElementCalibrationData) --------
setMethod(
  "..compute_data_element_estimates",
  signature(x = "familiarDataElementCalibrationData"),
  function(x, x_list = NULL, ...) {
    
    # It might be that x was only used to direct to this method.
    if (!is.null(x_list)) x <- x_list
    if (!is.list(x)) x <- list(x)
    
    # Identify the estimation types of the current data elements.
    estimation_type <- sapply(x, function(x) (x@estimation_type))
    
    # Determine the bootstrap_ci_method and the aggregation function
    if (any(estimation_type %in% c("bci", "bootstrap_confidence_interval"))) {
      
      # Point estimates need to be for all expected probability values,
      # which means that we have to compute it from the data.
      x <- .add_point_estimate_from_elements(
        x = x[estimation_type %in% c("bci", "bootstrap_confidence_interval")]
      )
    }
    
    # Check that x is not empty.
    if (is_empty(x)) return(NULL)
    
    return(callNextMethod(x_list = x))
  }
)



# ..compute_data_element_estimates (familiarDataElementCalibrationLinearFit) ----
setMethod(
  "..compute_data_element_estimates",
  signature(x = "familiarDataElementCalibrationLinearFit"),
  function(x, x_list = NULL, ...) {
    
    # Suppress NOTES due to non-standard evaluation in data.table
    p_value <- NULL
    
    # It might be that x was only used to direct to this method.
    if (!is.null(x_list)) x <- x_list
    if (!is.list(x)) x <- list(x)
    
    # Process data normally.
    y <- callNextMethod()
    
    # Identify the estimation types of the current data elements.
    estimation_type <- sapply(x, function(x) (x@estimation_type))
    
    # Determine the bootstrap_ci_method and the aggregation function
    if (any(estimation_type %in% c("bci", "bootstrap_confidence_interval"))) {
      
      # Obtain data from bootstrapped data.
      data <- x[estimation_type %in% c("bci", "bootstrap_confidence_interval")][[1L]]@data
      
    } else {
      data <- x[[1L]]@data
    }
    
    # Get the grouping column.
    grouping_column <- x[[1L]]@grouping_column
    
    # Compute p-value by grouping column.
    p_value <- data[
      ,
      list("p_value" = harmonic_p_value(p_value)),
      by = c(grouping_column)
    ]
    
    # Merge data into 
    y@data <- merge(
      x = y@data,
      y = p_value,
      by = grouping_column
    )
    
    # Update value column
    y@value_column <- setdiff(
      names(y@data),
      y@grouping_column
    )
    
    return(y)
  }
)


# ..compute_data_element_estimates (familiarDataElementCalibrationGoodnessOfFit) -----
setMethod(
  "..compute_data_element_estimates",
  signature(x = "familiarDataElementCalibrationGoodnessOfFit"),
  function(x, x_list = NULL, ...) {
    
    # Suppress NOTES due to non-standard evaluation in data.table
    p_value <- NULL
    
    # It might be that x was only used to direct to this method.
    if (!is.null(x_list)) x <- x_list
    if (!is.list(x)) x <- list(x)
    
    # Identify the estimation types of the current data elements.
    estimation_type <- sapply(x, function(x) (x@estimation_type))
    
    # Determine the bootstrap_ci_method and the aggregation function
    if (any(estimation_type %in% c("bci", "bootstrap_confidence_interval"))) {
      
      # Obtain data from bootstrapped data.
      x <- x[estimation_type %in% c("bci", "bootstrap_confidence_interval")][[1L]]
      
    } else {
      x <- x[[1L]]
    }
    
    # Get the grouping column.
    grouping_column <- x@grouping_column
    
    # Since samples overlap, combination methods that assume independence cannot
    # be used (e.g. Fisher's method). Hence we calculate a harmonic mean p value
    # according to Wilson (2019), 10.1073/pnas.1814092116.
    x@data <- x@data[
      ,
      list("p_value" = harmonic_p_value(p_value)),
      by = c(grouping_column)
    ]
    
    # Update value column
    x@value_column <- setdiff(names(x@data), x@grouping_column)
    
    return(x)
  }
)


# ..compute_data_element_estimates (familiarDataElementCalibrationDensity) -----
setMethod(
  "..compute_data_element_estimates",
  signature(x = "familiarDataElementCalibrationDensity"),
  function(x, x_list = NULL, ...) {
    
    # Suppress NOTES due to non-standard evaluation in data.table
    frequency <- expected <- NULL
    
    # It might be that x was only used to direct to this method.
    if (!is.null(x_list)) x <- x_list
    if (!is.list(x)) x <- list(x)
    
    # Identify the estimation types of the current data elements.
    estimation_type <- sapply(x, function(x) (x@estimation_type))
    
    # Determine the bootstrap_ci_method and the aggregation function
    if (any(estimation_type %in% c("bci", "bootstrap_confidence_interval"))) {
      
      # Obtain data from bootstrapped data.
      x <- x[estimation_type %in% c("bci", "bootstrap_confidence_interval")][[1L]]
      
    } else {
      x <- x[[1L]]
    }
    
    # Get the grouping column.
    grouping_column <- x@grouping_column
    
    # Compute frequency by grouping column.
    x@data <- x@data[
      ,
      list("frequency" = sum(frequency)),
      by = c(grouping_column)
    ]
    
    # Normalise frequency by all grouping columns minus expected.
    grouping_column <- setdiff(grouping_column, "expected")
    
    if (length(grouping_column) > 0L) {
      x@data <- x@data[
        ,
        "frequency" := frequency / sum(frequency),
        by = c(grouping_column)
      ][order(expected)]
      
    } else {
      x@data <- x@data[, "frequency" := frequency / sum(frequency)][order(expected)]
    }
    
    # Update value column
    x@value_column <- setdiff(names(x@data), x@grouping_column)
    
    return(x)
  }
)



# export_calibration_data (generic) --------------------------------------------

#'@title Extract and export calibration and goodness-of-fit tests.
#'
#'@description Extract and export calibration and goodness-of-fit tests for data
#'  in a familiarCollection.
#'
#'@inheritParams export_all
#'@inheritParams export_univariate_analysis_data
#'
#'@inheritDotParams extract_calibration_data
#'@inheritDotParams as_familiar_collection
#'
#'@details Data is usually collected from a `familiarCollection` object.
#'  However, you can also provide one or more `familiarData` objects, that will
#'  be internally converted to a `familiarCollection` object. It is also
#'  possible to provide a `familiarEnsemble` or one or more `familiarModel`
#'  objects together with the data from which data is computed prior to export.
#'  Paths to the previous files can also be provided.
#'
#'  All parameters aside from `object` and `dir_path` are only used if `object`
#'  is not a `familiarCollection` object, or a path to one.
#'
#'  Calibration tests are performed based on expected (predicted) and observed
#'  outcomes. For all outcomes, calibration-at-the-large and calibration slopes
#'  are determined. Furthermore, for all but survival outcomes, a repeated,
#'  randomised grouping Hosmer-Lemeshow test is performed. For survival
#'  outcomes, the Nam-D'Agostino and Greenwood-Nam-D'Agostino tests are
#'  performed.
#'
#'@return A list of data.tables (if `dir_path` is not provided), or nothing, as
#'  all data is exported to `csv` files.
#'@exportMethod export_calibration_data
#'@md
#'@rdname export_calibration_data-methods
setGeneric(
  "export_calibration_data",
  function(
    object,
    dir_path = NULL,
    aggregate_results = TRUE,
    export_collection = FALSE,
    ...
  ) {
    standardGeneric("export_calibration_data")
  }
)



# export_calibration_data (collection) -----------------------------------------

#'@rdname export_calibration_data-methods
setMethod(
  "export_calibration_data",
  signature(object = "familiarCollection"),
  function(
    object,
    dir_path = NULL,
    aggregate_results = TRUE,
    export_collection = FALSE,
    ...
  ) {
    
    # Make sure the collection object is updated.
    object <- update_object(object = object)
    
    # Obtain calibration data.
    calibration_data <- .export(
      x = object,
      data_slot = "calibration_data",
      dir_path = dir_path,
      aggregate_results = aggregate_results,
      type = "calibration",
      subtype = "data",
      object_class = "familiarDataElementCalibrationData"
    )
    
    # Obtain linear fit data
    linear_fit_data <- .export(
      x = object,
      data_slot = "calibration_data",
      dir_path = dir_path,
      aggregate_results = aggregate_results,
      type = "calibration",
      subtype = "at_large",
      object_class = "familiarDataElementCalibrationLinearFit"
    )
    
    # Obtain goodness of fit data
    gof_data <- .export(
      x = object,
      data_slot = "calibration_data",
      dir_path = dir_path,
      aggregate_results = aggregate_results,
      type = "calibration",
      subtype = "gof_test",
      object_class = "familiarDataElementCalibrationGoodnessOfFit"
    )
    
    # Obtain density data
    density_data <- .export(
      x = object,
      data_slot = "calibration_data",
      dir_path = dir_path,
      aggregate_results = aggregate_results,
      type = "calibration",
      subtype = "density",
      object_class = "familiarDataElementCalibrationDensity"
    )
    
    # Set data list.
    data_list <- list(
      "data" = calibration_data,
      "density" = density_data,
      "linear_test" = linear_fit_data,
      "gof_test" = gof_data
    )
    
    if (!is.null(dir_path)) data_list <- NULL
    if (export_collection) data_list <- c(data_list, list("collection" = object))
    
    return(data_list)
  }
)



# export_calibration_data (general) --------------------------------------------

#'@rdname export_calibration_data-methods
setMethod(
  "export_calibration_data",
  signature(object = "ANY"),
  function(
    object,
    dir_path = NULL,
    aggregate_results = TRUE,
    export_collection = FALSE,
    ...
  ) {
    
    # Attempt conversion to familiarCollection object.
    object <- do.call(
      as_familiar_collection,
      args = c(
        list(
          "object" = object,
          "data_element" = "calibration_data",
          "aggregate_results" = aggregate_results
        ),
        list(...)
      )
    )
    
    return(do.call(
      export_calibration_data,
      args = c(
        list(
          "object" = object,
          "dir_path" = dir_path,
          "aggregate_results" = aggregate_results,
          "export_collection" = export_collection
        ),
        list(...)
      )
    ))
  }
)