R/varform.R

In R4VN: Health Data Analysis and Publication-Ready Reporting

# ==========================================================================
# Variable transformation diagnostics (R4VN equivalent of Stata gladder)
# ==========================================================================

.r4vn_varform_transformations <- function(x, shift = c("none", "auto")) {
  shift <- match.arg(shift)
  z <- x
  add <- 0
  finite <- z[is.finite(z)]
  if (shift == "auto" && length(finite) && min(finite) <= 0) {
    add <- abs(min(finite)) + max(.Machine$double.eps^0.5, diff(range(finite)) * 1e-6)
  }
  zp <- z + add
  positive <- zp > 0 & is.finite(zp)
  nonnegative <- zp >= 0 & is.finite(zp)
  safe <- function(val, ok = is.finite(val)) { val[!ok] <- NA_real_; val }
  list(
    `Cubic (x^3)` = safe(z^3),
    `Square (x^2)` = safe(z^2),
    `Identity (x)` = safe(z),
    `Square root (sqrt x)` = safe(sqrt(pmax(zp, 0)), nonnegative),
    `Log (ln x)` = safe(log(zp), positive),
    `1 / sqrt(x)` = safe(1 / sqrt(zp), positive),
    `1 / x` = safe(1 / zp, positive),
    `1 / x^2` = safe(1 / (zp^2), positive),
    `1 / x^3` = safe(1 / (zp^3), positive),
    shift = add
  )
}

.r4vn_varform_stats <- function(z, label, digits = 3, p_digits = 3) {
  z <- z[is.finite(z)]
  if (!length(z)) {
    return(data.frame(Transformation = label, n = 0L, Mean = NA, SD = NA, Median = NA,
                      Min = NA, Max = NA, Skewness = NA, Kurtosis = NA,
                      Shapiro.W = NA, Shapiro.p = NA, stringsAsFactors = FALSE, check.names = FALSE))
  }
  sk <- .r4vn_skew_kurt(z)
  sh <- if (length(z) >= 3L && length(z) <= 5000L && length(unique(z)) > 1L) tryCatch(stats::shapiro.test(z), error = function(e) NULL) else NULL
  data.frame(
    Transformation = label, n = length(z), Mean = .r4vn_num(mean(z), digits),
    SD = .r4vn_num(stats::sd(z), digits), Median = .r4vn_num(stats::median(z), digits),
    Min = .r4vn_num(min(z), digits), Max = .r4vn_num(max(z), digits),
    Skewness = .r4vn_num(sk[["skewness"]], digits), Kurtosis = .r4vn_num(sk[["kurtosis"]], digits),
    Shapiro.W = if (is.null(sh)) NA else .r4vn_num(unname(sh$statistic), digits),
    Shapiro.p = if (is.null(sh)) NA else .r4vn_p(sh$p.value, p_digits),
    stringsAsFactors = FALSE, check.names = FALSE
  )
}

.r4vn_varform_draw <- function(trans, variable, context = NULL, bins = "Sturges",
                               color = NULL, palette = "default", alpha = 0.8,
                               normal_color = "black", normal_lty = 1, normal_lwd = 2,
                               theme = "journal", size = 10) {
  old <- graphics::par(no.readonly = TRUE)
  on.exit(graphics::par(old), add = TRUE)
  graphics::par(mfrow = c(3, 3), mar = c(3.1, 3.4, 3.1, 0.8), mgp = c(2.1, 0.55, 0), tcl = -0.2,
                cex.axis = size / 11, cex.lab = size / 11, cex.main = size / 11)
  cols <- .r4vn_colors(9, color, palette, alpha)
  nm <- setdiff(names(trans), "shift")
  for (i in seq_along(nm)) {
    z <- trans[[nm[i]]]; z <- z[is.finite(z)]
    if (length(z) < 2L || length(unique(z)) < 2L) {
      graphics::plot.new(); graphics::title(main = nm[i]); graphics::text(0.5, 0.5, "Not estimable", cex = 0.9)
      next
    }
    h <- graphics::hist(z, breaks = bins, probability = TRUE, col = cols[i], border = "white",
                        main = nm[i], xlab = "", ylab = "Density")
    xx <- seq(min(h$breaks), max(h$breaks), length.out = 300)
    graphics::lines(xx, stats::dnorm(xx, mean(z), stats::sd(z)), col = normal_color,
                    lty = normal_lty, lwd = normal_lwd)
    sk <- .r4vn_skew_kurt(z)
    graphics::mtext(sprintf("n=%d; skew=%s; kurt=%s", length(z), .r4vn_num(sk[["skewness"]], 2), .r4vn_num(sk[["kurtosis"]], 2)),
                    side = 3, line = 0.15, cex = 0.7)
  }
  main <- paste0("Transformation ladder: ", variable)
  if (!is.null(context) && nzchar(context)) main <- paste0(main, " | ", context)
  graphics::mtext(main, outer = FALSE, side = 3, line = -1.1, adj = 0, cex = 0.8, font = 2)
}

#' Explore transformations of continuous variables
#' @usage varform(x = NULL, vars = NULL, data = NULL, by = NULL, shift = c("none", "auto"), bins = "Sturges", color = NULL, palette = "default", alpha = 0.8, normal_color = "black", normal_lty = 1, normal_lwd = 2, theme = "journal", size = 10, file = NULL, width = 9, height = 8, dpi = 300, show = TRUE, bg = "white", digits = 3, p_digits = 3)
#'
#' @description
#' `varform()` is an R4VN transformation diagnostic inspired by Stata's
#' `gladder`. It compares the nine Tukey ladder transformations using
#' histograms with fitted normal curves and returns descriptive distribution
#' statistics for every transformation.
#'
#' @param x One numeric variable. Can be omitted when `vars` is supplied.
#' @param vars Optional `vars(...)` selection of several numeric variables.
#' @param data Data frame. If `NULL`, the active R4VN data frame is used.
#' @param by Optional grouping specification. `by = vars(region, sex)` treats
#'   `region` as an outer stratum and `sex` as the innermost grouping variable.
#' @param shift `"none"` keeps the original scale; transformations requiring
#'   positive values are unavailable when the data are nonpositive. `"auto"`
#'   adds the smallest constant required to make the transformed base positive.
#' @param bins Histogram break specification passed to `hist()`.
#' @param color,palette,alpha Histogram appearance.
#' @param normal_color,normal_lty,normal_lwd Appearance of the fitted normal curve.
#' @param theme,size Graph theme/font controls.
#' @param file Optional graph file. With several graphs, R4VN appends a sequence
#'   number to the file name.
#' @param width,height,dpi,show,bg Standard graph export/display options.
#' @param digits,p_digits Formatting digits.
#' @return An `r4vn_stat` object containing transformation diagnostics; its raw
#'   component also stores the transformed values. Graphs are drawn as a 3 x 3
#'   transformation ladder.
#' @export
#' @examples
#' d <- data.frame(weight = c(29,26,13,23,23,25,17,22,17,19,12,26,30,30,18,14,12,26,17,18))
#' varform(weight, data = d)
#' varform(vars = vars(weight), data = d, shift = "auto")
varform <- function(x = NULL, vars = NULL, data = NULL, by = NULL,
                    shift = c("none", "auto"), bins = "Sturges",
                    color = NULL, palette = "default", alpha = 0.8,
                    normal_color = "black", normal_lty = 1, normal_lwd = 2,
                    theme = "journal", size = 10, file = NULL, width = 9, height = 8,
                    dpi = 300, show = TRUE, bg = "white", digits = 3, p_digits = 3) {
  call <- match.call(); env <- parent.frame(); d <- .r4vn_stat_data(data); shift <- match.arg(shift)
  x_expr <- substitute(x); vars_expr <- substitute(vars); by_expr <- substitute(by)
  selected <- character()
  if (!missing(vars) && !.r4vn_expr_is_null(vars_expr)) selected <- .r4vn_vars_spec_names(vars_expr, d, env, "vars", numeric_only = TRUE)
  if (!missing(x) && !.r4vn_expr_is_null(x_expr)) selected <- unique(c(.r4vn_resolve_name_spec(x_expr, d, env, "x", multiple = FALSE), selected))
  if (!length(selected)) stop("Supply `x` or `vars = vars(...)`.", call. = FALSE)
  bad <- selected[!vapply(d[selected], is.numeric, logical(1))]
  if (length(bad)) stop("`varform()` requires numeric variables: ", paste(bad, collapse = ", "), call. = FALSE)
  bspec <- .r4vn_by_spec(by_expr, d, env, allow_null = TRUE)
  splitvars <- bspec$all
  ids <- .r4vn_strata_indices(d, splitvars)
  if (!length(splitvars)) ids <- list(Overall = seq_len(nrow(d)))

  results <- list(); labels <- character(); raw <- list(); graph_index <- 0L
  for (nm in selected) {
    for (idx in ids) {
      graph_index <- graph_index + 1L
      context <- if (length(splitvars)) .r4vn_stratum_label(d, splitvars, idx) else "Overall"
      trans <- .r4vn_varform_transformations(d[[nm]][idx], shift)
      sh <- trans$shift; trans$shift <- NULL
      tab <- do.call(rbind, Map(function(z, lab) .r4vn_varform_stats(z, lab, digits, p_digits), trans, names(trans)))
      tab$Variable <- nm; tab$Stratum <- context; tab$Shift.added <- .r4vn_num(sh, digits)
      tab <- tab[c("Variable", "Stratum", "Transformation", "n", "Mean", "SD", "Median", "Min", "Max", "Skewness", "Kurtosis", "Shapiro.W", "Shapiro.p", "Shift.added")]
      results[[length(results) + 1L]] <- tab
      labels <- c(labels, paste(nm, context, sep = " | "))
      raw[[labels[length(labels)]]] <- trans
      target <- file
      if (!is.null(file) && (length(selected) * length(ids) > 1L)) {
        ext <- tools::file_ext(file); stem <- sub(paste0("\\.", ext, "$"), "", file, ignore.case = TRUE)
        target <- paste0(stem, "-", graph_index, if (nzchar(ext)) paste0(".", ext) else "")
      }
      draw <- local({ tt <- trans; vv <- nm; cc <- context; function() .r4vn_varform_draw(tt, vv, cc, bins, color, palette, alpha, normal_color, normal_lty, normal_lwd, theme, size) })
      .r4vn_render(draw, file = target, show = show, width = width, height = height, dpi = dpi, bg = bg)
    }
  }
  table <- do.call(rbind, results); rownames(table) <- NULL
  note <- "The nine panels follow Tukey's transformation ladder. Transformations involving log or reciprocal powers require positive values unless shift = 'auto'."
  .r4vn_show(.r4vn_result("Variable transformation diagnostics", list("Transformation diagnostics" = table),
    notes = note, raw = list(transformations = raw, variables = selected, by = bspec), call = call), show = show, console = FALSE)
}