R/advanced_validation.R

In imedR: Guarded Evidence Synthesis for the Indian Microplastics Evidence Database

#' Guard analytical size-window compatibility
#'
#' Uses linked SIZE_METHODS records to identify unresolved or heterogeneous
#' analytical/reporting lower bounds. It never substitutes mesh/filter pore
#' size for an analytical lower size limit.
#' @param x An imed_db.
#' @param sample_id Optional sample IDs.
#' @param action error, warn or return.
#' @return A `tibble` with columns `compatible` (logical) and `issue` (character;
#'   `"None"` when no problem is found) if `action = "return"`. For `action =
#'   "error"` or `"warn"` the function returns invisibly a single logical value
#'   (`TRUE` if the size windows are compatible) and stops or warns when they are
#'   not.
#' @examples
#' db <- imed_data(tables = "SIZE_METHODS", validate = FALSE)
#' imed_check_size_window(db, action = "return")
#' @export
imed_check_size_window <- function(x,sample_id=NULL,action=c("error","warn","return")) {
  action=match.arg(action)
  if (!inherits(x,"imed_db")) .imed_stop("Requires an imed_db.")
  z=x$SIZE_METHODS
  if (!is.null(sample_id)) z=z[z$Sample_ID %in% sample_id,,drop=FALSE]
  issues=character()
  al=.imed_num(z$Analytical_Lower_Limit_um)
  oc_col=.imed_col(z,c("Operational_Reporting_Cutoff_um","Operational_Cutoff_um"))
  oc=.imed_num(z[[oc_col]])
  known=unique(al[is.finite(al)])
  if (length(known)>1) issues=c(issues,paste("Multiple analytical lower limits:",paste(known,collapse=", "),"um"))
  if (any(!is.finite(al) & !is.finite(oc)))
    issues=c(issues,paste(sum(!is.finite(al)&!is.finite(oc)),"size-method record(s) have neither analytical lower limit nor operational cutoff."))
  ok=length(issues)==0
  out=tibble::tibble(compatible=ok,issue=if(ok)"None" else issues)
  if(action=="return") return(out)
  if(!ok && action=="error") .imed_stop("Size-window synthesis blocked:\n- ",paste(issues,collapse="\n- "))
  if(!ok && action=="warn") .imed_warn(paste(issues,collapse=" "))
  invisible(ok)
}

#' Sensitivity summary by evidence restriction
#'
#' Compares all eligible records with a stricter subset when evidence fields
#' are available. This is descriptive sensitivity analysis, not causal proof.
#' @param x An `imed_db`, or a measurement-level data frame such as the output of
#'   `imed_analysis_set()`. The comparability guard is applied, so the
#'   records must belong to one compatible comparability class.
#' @param value Name of the numeric measurement field to summarise (default `"Mean"`).
#' @return A `tibble` with columns `stratum`, `n` (number of finite values), `mean` and
#'   `median`. The first row summarises all compatible records; a second row for
#'   strict-synthesis-eligible records is added when the field
#'   `Strict_Synthesis_Eligibility` is present.
#' @examples
#' toy <- data.frame(
#'   Study_ID = paste0("S", 1:12),
#'   Comparability_Class = "CLASS_A",
#'   Reported_Unit = "items/L",
#'   Metric_Type = "Abundance",
#'   Strict_Synthesis_Eligibility = "YES",
#'   Group = rep(c("g1", "g2"), 6),
#'   Mean = c(2.1, 3.4, 1.8, 4.0, 2.9, 3.1, 2.4, 3.8, 1.9, 3.3, 2.7, 4.2),
#'   SD = rep(c(0.6, 0.9), 6),
#'   N = rep(c(10, 12, 8), 4)
#' )
#' imed_sensitivity(toy)
#' @export
imed_sensitivity <- function(x,value="Mean") {
  d=.imed_table(x,"MEASUREMENTS")
  imed_check_comparability(d,action="error")
  if(!value %in% names(d)) .imed_stop("Value field absent: ",value)
  y=.imed_num(d[[value]])
  base=tibble::tibble(stratum="All compatible",n=sum(is.finite(y)),
    mean=mean(y,na.rm=TRUE),median=stats::median(y,na.rm=TRUE))
  if("Strict_Synthesis_Eligibility" %in% names(d)) {
    keep=grepl("^YES",as.character(d$Strict_Synthesis_Eligibility),ignore.case=TRUE)
    yy=y[keep]
    strict=tibble::tibble(stratum="Strict synthesis eligible",n=sum(is.finite(yy)),
      mean=mean(yy,na.rm=TRUE),median=stats::median(yy,na.rm=TRUE))
    return(dplyr::bind_rows(base,strict))
  }
  base
}

#' Leave-one-study-out influence analysis for a guarded mean meta-analysis
#' @param x A guarded analysis set that passes `imed_meta_eligibility()`, as for
#'   `imed_meta_mean()`.
#' @param method Between-study variance estimator passed to `metafor` (default `"REML"`).
#' @return A `tibble` with one row per omitted study. The first column, `omitted`,
#'   identifies the study that was left out; the remaining columns are the
#'   pooled estimate and heterogeneity statistics returned by
#'   `metafor::leave1out()` after that study is removed.
#' @examples
#' if (requireNamespace("metafor", quietly = TRUE)) {
#'   toy <- data.frame(
#'     Study_ID = paste0("S", 1:12),
#'     Comparability_Class = "CLASS_A",
#'     Reported_Unit = "items/L",
#'     Metric_Type = "Abundance",
#'     Strict_Synthesis_Eligibility = "YES",
#'     Group = rep(c("g1", "g2"), 6),
#'     Mean = c(2.1, 3.4, 1.8, 4.0, 2.9, 3.1, 2.4, 3.8, 1.9, 3.3, 2.7, 4.2),
#'     SD = rep(c(0.6, 0.9), 6),
#'     N = rep(c(10, 12, 8), 4)
#'   )
#'   imed_leave_one_out(toy)
#' }
#' @export
imed_leave_one_out <- function(x,method="REML") {
  if(!requireNamespace("metafor",quietly=TRUE)) .imed_stop("Install 'metafor'.")
  fit=imed_meta_mean(x,method=method)
  z=metafor::leave1out(fit)
  tibble::as_tibble(as.data.frame(z),rownames="omitted")
}

#' Forest plot for an imedR meta-analysis
#' @param fit A fitted `metafor` `rma` object, for example from `imed_meta_mean()`.
#' @param ... Further arguments passed to `metafor::forest()`.
#' @return The fitted model `fit`, returned invisibly. The function is called for its
#'   side effect of drawing a forest plot on the active graphics device.
#' @examples
#' if (requireNamespace("metafor", quietly = TRUE)) {
#'   toy <- data.frame(
#'     Study_ID = paste0("S", 1:12),
#'     Comparability_Class = "CLASS_A",
#'     Reported_Unit = "items/L",
#'     Metric_Type = "Abundance",
#'     Strict_Synthesis_Eligibility = "YES",
#'     Group = rep(c("g1", "g2"), 6),
#'     Mean = c(2.1, 3.4, 1.8, 4.0, 2.9, 3.1, 2.4, 3.8, 1.9, 3.3, 2.7, 4.2),
#'     SD = rep(c(0.6, 0.9), 6),
#'     N = rep(c(10, 12, 8), 4)
#'   )
#'   fit <- imed_meta_mean(toy)
#'   imed_forest(fit)
#' }
#' @export
imed_forest <- function(fit,...) {
  if(!inherits(fit,"rma")) .imed_stop("fit must be a metafor rma object.")
  metafor::forest(fit,...)
  invisible(fit)
}

#' Funnel plot for an imedR meta-analysis
#'
#' A funnel plot is diagnostic only and should not be interpreted as a direct
#' test of publication bias, especially with small or heterogeneous evidence.
#' @param fit A fitted `metafor` `rma` object, for example from `imed_meta_mean()`.
#' @param ... Further arguments passed to `metafor::funnel()`.
#' @return The fitted model `fit`, returned invisibly. The function is called for its
#'   side effect of drawing a funnel plot on the active graphics device.
#' @examples
#' if (requireNamespace("metafor", quietly = TRUE)) {
#'   toy <- data.frame(
#'     Study_ID = paste0("S", 1:12),
#'     Comparability_Class = "CLASS_A",
#'     Reported_Unit = "items/L",
#'     Metric_Type = "Abundance",
#'     Strict_Synthesis_Eligibility = "YES",
#'     Group = rep(c("g1", "g2"), 6),
#'     Mean = c(2.1, 3.4, 1.8, 4.0, 2.9, 3.1, 2.4, 3.8, 1.9, 3.3, 2.7, 4.2),
#'     SD = rep(c(0.6, 0.9), 6),
#'     N = rep(c(10, 12, 8), 4)
#'   )
#'   fit <- imed_meta_mean(toy)
#'   imed_funnel(fit)
#' }
#' @export
imed_funnel <- function(fit,...) {
  if(!inherits(fit,"rma")) .imed_stop("fit must be a metafor rma object.")
  metafor::funnel(fit,...)
  invisible(fit)
}

#' Meta-regression with minimum-information guard
#'
#' Requires at least 10 independent effects by default and at least 5 effects
#' per model coefficient. These are guardrails, not guarantees of adequate power.
#' @param x A guarded analysis set that passes `imed_meta_eligibility()`.
#' @param moderators Character vector of field names in the measurement table to use as
#'   moderators.
#' @param min_k Minimum number of independent effects required (default 10). At least
#'   5 effects per model coefficient are also required.
#' @param method Between-study variance estimator passed to `metafor` (default `"REML"`).
#' @return A fitted random-effects meta-regression model of class `c("rma.uni",
#'   "rma")` from `metafor::rma.uni()`, with the study identifiers as labels. The
#'   coefficients describe the change in the mean outcome associated with each
#'   moderator.
#' @examples
#' if (requireNamespace("metafor", quietly = TRUE)) {
#'   toy <- data.frame(
#'     Study_ID = paste0("S", 1:12),
#'     Comparability_Class = "CLASS_A",
#'     Reported_Unit = "items/L",
#'     Metric_Type = "Abundance",
#'     Strict_Synthesis_Eligibility = "YES",
#'     Group = rep(c("g1", "g2"), 6),
#'     Mean = c(2.1, 3.4, 1.8, 4.0, 2.9, 3.1, 2.4, 3.8, 1.9, 3.3, 2.7, 4.2),
#'     SD = rep(c(0.6, 0.9), 6),
#'     N = rep(c(10, 12, 8), 4)
#'   )
#'   imed_meta_regression(toy, moderators = "Group")
#' }
#' @export
imed_meta_regression <- function(x,moderators,min_k=10,method="REML") {
  if(!requireNamespace("metafor",quietly=TRUE)) .imed_stop("Install 'metafor'.")
  d=.imed_table(x,"MEASUREMENTS")
  e=imed_meta_eligibility(d)
  if(!all(e$eligible)) .imed_stop("Meta-regression blocked by base eligibility.")
  moderators=as.character(moderators)
  miss=setdiff(moderators,names(d))
  if(length(miss)) .imed_stop("Moderator field(s) absent: ",paste(miss,collapse=", "))
  k=nrow(d); p=length(moderators)+1
  if(k<min_k || k<5*p) .imed_stop("Insufficient independent effects for guarded meta-regression: k=",k,", coefficients=",p,".")
  yi=.imed_num(d$Mean); vi=(.imed_num(d$SD)^2)/.imed_num(d$N)
  form=stats::as.formula(paste("~",paste(moderators,collapse="+")))
  metafor::rma.uni(yi=yi,vi=vi,mods=form,data=d,method=method,slab=d$Study_ID)
}

#' Validation report for manuscript/release audit
#'
#' Returns database integrity, HOLD burden, missingness and comparability counts.
#' @param x An `imed_db` object, as returned by `imed_data()`.
#' @return A named list with five elements: `integrity` (the `tibble` returned by
#'   `imed_validate()`), `cardinalities` (table, row and column counts from
#'   `imed_tables()`), `hold_summary` (a one-row `tibble` with the number of
#'   measurements, the number flagged as HOLD and the number eligible),
#'   `comparability` (counts of records by `Comparability_Class` and
#'   `Strict_Synthesis_Eligibility`) and `measurement_missingness` (the
#'   `tibble` returned by `imed_missingness()` for the measurements table).
#' @examples
#' db <- imed_data(validate = FALSE)
#' rep <- imed_validation_report(db)
#' names(rep)
#' rep$hold_summary
#' @export
imed_validation_report <- function(x) {
  if(!inherits(x,"imed_db")) .imed_stop("Requires an imed_db.")
  v=imed_validate(x)
  h=imed_flag_holds(x)
  cmp=x$MEASUREMENT_COMPARABILITY
  list(
    integrity=v,
    cardinalities=imed_tables(x),
    hold_summary=tibble::tibble(
      measurements=nrow(h),
      holds=sum(h$.imed_hold,na.rm=TRUE),
      eligible=nrow(h)-sum(h$.imed_hold,na.rm=TRUE)
    ),
    comparability=dplyr::count(cmp,Comparability_Class,Strict_Synthesis_Eligibility,sort=TRUE),
    measurement_missingness=imed_missingness(x,"MEASUREMENTS")
  )
}