# Celeus - reproduce.R (run from this package folder; bundled data in data/) # Test: Descriptive table (Table 1) | engine 0.41.7 | catalog c40-dc08 # Research context: not provided for this run. # Reproduces the recorded result deterministically. set.seed(20261101) .data_path <- "data/table-1-participant-characteristics.csv" if (!file.exists(.data_path)) stop("Run reproduce.R from the package folder (the bundled data is in data/).") data <- readr::read_csv(.data_path, show_col_types = FALSE) # needs the readr package clean_cells <- function (df) { ctrl_collapse <- "[\r\n\t\v\f]+" zero_width_chars <- intToUtf8(c(8203L, 65279L), multiple = TRUE) uni_space_chars <- intToUtf8(c(160L, 8239L, 8199L), multiple = TRUE) clean_text <- function(x) { out <- x for (ch in zero_width_chars) out <- gsub(ch, "", out, fixed = TRUE) out <- gsub(ctrl_collapse, " ", out) for (ch in uni_space_chars) out <- gsub(ch, " ", out, fixed = TRUE) trimws(gsub(" +", " ", out)) } changed_cols <- character(0) cells_cleaned <- 0L retyped <- character(0) for (j in seq_along(df)) { col <- df[[j]] if (!is.character(col)) next cleaned <- clean_text(col) n_changed <- sum(!is.na(col) & !is.na(cleaned) & col != cleaned, na.rm = TRUE) if (n_changed > 0L) { changed_cols <- c(changed_cols, names(df)[j]) cells_cleaned <- cells_cleaned + n_changed } df[[j]] <- cleaned non_blank <- cleaned[!is.na(cleaned) & nzchar(cleaned)] if (length(non_blank) && all(grepl("^[+-]?(0|[1-9][0-9]*)([.][0-9]+)?$", non_blank))) { parsed <- suppressWarnings(as.numeric(cleaned)) was_blank <- is.na(cleaned) | !nzchar(cleaned) lost <- (is.na(parsed) | is.infinite(parsed)) & !was_blank if (!any(lost)) { finite <- parsed[!is.na(parsed)] all_whole <- all(finite == floor(finite)) in_range <- all(abs(finite) <= .Machine$integer.max) df[[j]] <- if (all_whole && in_range) as.integer(parsed) else parsed retyped <- c(retyped, names(df)[j]) } } } names(df) <- clean_text(names(df)) if (length(changed_cols)) { attr(df, "ingest_clean_columns") <- changed_cols attr(df, "ingest_clean_cells_cleaned") <- cells_cleaned } if (length(retyped)) attr(df, "ingest_clean_retyped") <- retyped df } data <- clean_cells(data) # --- The descriptive table, emitted verbatim so this script needs no package --------------- # These are the engine's own functions, copied in so the table reproduces from base R alone. CELEUS <- list(shapiro_min_n = 3L, shapiro_max_n = 5000L, table_one_welch_min_stratum_n = 10L, assumption_check_alpha = 0.05, expected_count_min = 5L, default_p_adjust = "holm", max_levels_report = 20L, high_cardinality_ratio = 0.5, small_cell_threshold = 11L, censored_pct_heavy = 0.5) `%||%` <- function (a, b) if (is.null(a) || length(a) == 0) b else a `.fmt` <- function (v) { if (is.null(v) || length(v) == 0) return("NA") v <- suppressWarnings(as.numeric(v[1])) if (is.na(v)) "NA" else format(signif(v, 4)) } `.fmt_pval` <- function (p) { if (is.null(p) || length(p) == 0 || is.na(suppressWarnings(as.numeric(p[1])))) return("NA") p <- as.numeric(p[1]) if (p < 0.001) "< 0.001" else format(signif(p, 3)) } `.fmt_p` <- function (p) { v <- .fmt_pval(p) if (identical(v, "NA")) "NA" else if (startsWith(v, "<")) paste("p", v) else paste("p =", v) } `.col_type` <- function (x) { if (inherits(x, c("Date", "POSIXct", "POSIXt"))) return("date") if (is.factor(x)) return("factor") if (is.logical(x)) return("logical") if (is.integer(x)) return("integer") if (is.numeric(x)) return("numeric") "character" } `.ID_NAME_HINTS` <- c("id", "caseid", "casenum", "caseno", "recordid", "recordno", "rownum", "seqno", "serial") `.ID_NAME_HINT_RE` <- "^(.*[._ -])?(id|caseid|casenum|caseno|recordid|recordno|rownum|seqno|serial)$" `.is_id_named` <- function (name) { grepl(.ID_NAME_HINT_RE, tolower(name)) || .phi_norm(name) %in% .ID_ROW_HINT_TOKENS_FOLDED } `.role_hint` <- function (type, n_distinct, ratio, name = "") { if (n_distinct <= 1) return("constant") if (type %in% c("numeric", "integer")) { if (n_distinct == 2) return("binary") if (ratio == 1 && .is_id_named(name)) return("identifier_or_text") return("continuous") } if (type == "date") return("date") high_card <- (n_distinct > CELEUS$max_levels_report) || (ratio > CELEUS$high_cardinality_ratio && n_distinct > CELEUS$small_cell_threshold) if (type != "logical" && high_card) return("identifier_or_text") if (n_distinct == 2) "binary" else "categorical" } `check_normality` <- function (x) { x <- x[!is.na(x)] n <- length(x) if (n < CELEUS$shapiro_min_n || n > CELEUS$shapiro_max_n) { return(list(id = ASSUMPTION_NORMALITY, check = "normality (Shapiro-Wilk)", statistic = NA_real_, p_value = NA_real_, passed = NA, note = sprintf(ASSUMPTION_NOTES[["shapiro_range"]], n, CELEUS$shapiro_min_n, CELEUS$shapiro_max_n))) } if (stats::sd(x) == 0) { return(list(id = ASSUMPTION_NORMALITY, check = "normality (Shapiro-Wilk)", statistic = NA_real_, p_value = NA_real_, passed = NA, note = ASSUMPTION_NOTES[["no_variance"]])) } w <- stats::shapiro.test(x) list(id = ASSUMPTION_NORMALITY, check = "normality (Shapiro-Wilk)", statistic = unname(w$statistic), p_value = w$p.value, passed = w$p.value > CELEUS$assumption_check_alpha, note = if (w$p.value > CELEUS$assumption_check_alpha) ASSUMPTION_NOTES[["normal"]] else ASSUMPTION_NOTES[["non_normal"]]) } `ASSUMPTION_NORMALITY` <- "normality" `ASSUMPTION_NOTES` <- c(shapiro_range = "n=%d outside Shapiro-Wilk valid range (%d-%d)", no_variance = "all values identical (no variance) - not evaluable", not_computed = "could not be computed", calibration_ok = "no evidence of poor calibration (observed event rates track the predicted risks)", calibration_poor = "evidence of poor calibration: observed event rates diverge from the predicted risks - the odds ratios may still be usable, but a predicted probability from this model is not", calibration_sparse = "- caution: at least one risk group expects fewer than %d events or fewer than %d non-events, where the chi-squared approximation behind this p-value is unreliable", coxzph_failed = "cox.zph could not be computed", sphericity_corrected = "F, df and p use the Greenhouse-Geisser correction, so a sphericity violation is already accounted for", normal = "no evidence against normality", non_normal = "evidence of non-normality", context_rank = "- normality is not assumed by this rank-based test", context_signed_rank = "- the signed-rank test assumes the differences are symmetric, not normal", equal_variances = "no evidence of unequal variances", unequal_variances = "evidence of unequal variances", expected_ok = "all expected counts >= %d", expected_low = "min expected count %.2f < %d; prefer Fisher's exact", expected_low_exact = "min expected count %.2f < %d; the exact test does not rely on this approximation", discordant_ok = "%d discordant pairs (the test uses only these pairs)", discordant_low = "only %d discordant pairs (< %d); with so few, this test cannot reach significance however the pairs split", discordant_none = "no discordant pairs; the data carry no information about a change, and the p-value is reported as 1", discordant_approximation = "- this method uses a chi-squared approximation, and passing this check says nothing about how accurate it is here; the mid-p and exact methods use the binomial distribution of the discordant pairs instead (mid-p is not an exact test and its error rate can slightly exceed the stated level)", slopes_common = "no evidence a covariate's effect differs by group; the common-slope ANCOVA model is appropriate", slopes_differ = "evidence that a covariate's effect differs by group - the common-slope adjusted comparison may not be valid; interpret it with caution", dispersion_ok = "residual deviance/df = %.2f; no strong evidence of overdispersion", dispersion_high = "residual deviance/df = %.2f (> %.1f); evidence of overdispersion - a negative-binomial model is preferable to Poisson", dispersion_high_handled = "residual deviance/df of a Poisson fit = %.2f (> %.1f); evidence of overdispersion, which this negative-binomial model accounts for by estimating a dispersion parameter", prop_odds_ok = "no evidence against the proportional-odds assumption (the predictors' effects look consistent across cut points)", prop_odds_violated = "evidence the proportional-odds assumption is violated (at least one predictor's effect differs across cut points) - interpret the single odds ratio with caution", n_ok = "adequate n", n_small = "small sample; interpret with caution", events_ok = "adequate events", events_few = "few events; the analysis may be underpowered", epv_ok = "at or above the rule-of-thumb floor (a convention, not a guarantee of a stable model)", epv_few = "below the rule-of-thumb floor; the model may be overfit", calibration_too_few = "too few distinct predicted risks to assess calibration", calibration_degenerate = "the model's predicted risks are degenerate in at least one risk group, so Hosmer-Lemeshow is undefined", ph_ok = "no evidence against proportional hazards", ph_violated = "evidence the proportional-hazards assumption is violated" ) `.ID_ROW_HINT_TOKENS_FOLDED` <- c("hastano", "dosyano", "fallnr", "patientennummer") `.PHI_FOLD_MAP` <- c(Ç = "C", ç = "c", Ğ = "G", ğ = "g", İ = "I", ı = "i", Ş = "S", ş = "s", ß = "ss", À = "A", Á = "A", Â = "A", Ã = "A", Ä = "A", Å = "A", à = "a", á = "a", â = "a", ã = "a", ä = "a", å = "a", È = "E", É = "E", Ê = "E", Ë = "E", è = "e", é = "e", ê = "e", ë = "e", Ì = "I", Í = "I", Î = "I", Ï = "I", ì = "i", í = "i", î = "i", ï = "i", Ñ = "N", ñ = "n", Ò = "O", Ó = "O", Ô = "O", Õ = "O", Ö = "O", ò = "o", ó = "o", ô = "o", õ = "o", ö = "o", Ù = "U", Ú = "U", Û = "U", Ü = "U", ù = "u", ú = "u", û = "u", ü = "u", Ý = "Y", ý = "y", ÿ = "y") `.phi_fold` <- function (x) { for (from in names(.PHI_FOLD_MAP)) x <- gsub(from, .PHI_FOLD_MAP[[from]], x, fixed = TRUE) x } `.phi_norm` <- function (x) gsub("[^a-z0-9]", "", tolower(.phi_fold(x))) `.CLIN_CI_LEVEL` <- 0.95 `PROPORTION_INTERVAL_WILSON` <- "wilson" `PROPORTION_INTERVAL_CLOPPER_PEARSON` <- "clopper_pearson" `PROPORTION_INTERVAL_METHODS` <- c("wilson", "clopper_pearson") `.wilson_ci` <- function (k, n, level = .CLIN_CI_LEVEL) { if (is.null(n) || is.na(n) || n <= 0 || is.null(k) || is.na(k)) return(c(NA_real_, NA_real_)) z <- stats::qnorm(1 - (1 - level)/2) denom <- n + z^2 centre <- (k + z^2/2)/denom half <- z/denom * sqrt(k * (n - k)/n + z^2/4) c(max(0, centre - half), min(1, centre + half)) } `.binom_ci` <- function (x, n, level = .CLIN_CI_LEVEL) { if (is.na(n) || n <= 0) return(c(NA_real_, NA_real_)) as.numeric(stats::binom.test(x, n, conf.level = level)$conf.int) } `.proportion_ci` <- function (k, n, method = PROPORTION_INTERVAL_WILSON, level = .CLIN_CI_LEVEL) { if (identical(method, PROPORTION_INTERVAL_CLOPPER_PEARSON)) .binom_ci(k, n, level) else .wilson_ci(k, n, level) } `valid_p_adjust` <- function (test_id) { if (identical(test_id, "kruskal_wallis") || identical(test_id, "peto_peto")) c("holm", "BH", "bonferroni", "none") else if (identical(test_id, "table_one_full")) c("holm", "BH", "bonferroni", "none") else "single-step" } `P_ADJUST_METHODS` <- function () valid_p_adjust("kruskal_wallis") `.adjust_rows` <- function (p, method = CELEUS$default_p_adjust) { if (!method %in% P_ADJUST_METHODS()) stop("unsupported p-adjust method: ", method, call. = FALSE) p <- as.numeric(p) k <- sum(!is.na(p)) if (k == 0L) return(rep(NA_real_, length(p))) stats::p.adjust(p, method = method, n = k) } `.P_ADJUST_DISPLAY_NAMES` <- c(holm = "Holm", BH = "Benjamini-Hochberg", bonferroni = "Bonferroni" ) `test_catalog` <- function () { list(welch_t = list(label = "Welch's t-test", label_key = "test.welch_t.label", kind = "two_group_continuous", when = "Two independent groups, continuous outcome, unequal variances (default two-group test)."), student_t = list(label = "Student's t-test", label_key = "test.student_t.label", kind = "two_group_continuous", when = "Two independent groups, continuous outcome, equal variances assumed."), paired_t = list(label = "Paired t-test", label_key = "test.paired_t.label", kind = "paired_continuous", when = "Two paired/repeated continuous measurements, normal differences."), mann_whitney = list(label = "Mann-Whitney U test", label_key = "test.mann_whitney.label", kind = "two_group_continuous", when = "Two independent groups, continuous/ordinal outcome, non-normal."), wilcoxon_signed_rank = list(label = "Wilcoxon signed-rank test", label_key = "test.wilcoxon_signed_rank.label", kind = "paired_continuous", when = "Two paired measurements, non-normal differences."), one_way_anova = list(label = "One-way ANOVA", label_key = "test.one_way_anova.label", kind = "multi_group_continuous", when = "3+ independent groups, continuous outcome, equal variances."), welch_anova = list(label = "Welch's ANOVA", label_key = "test.welch_anova.label", kind = "multi_group_continuous", when = "3+ independent groups, continuous outcome, unequal variances."), kruskal_wallis = list(label = "Kruskal-Wallis test", label_key = "test.kruskal_wallis.label", kind = "multi_group_continuous", when = "3+ independent groups, continuous/ordinal outcome, non-normal."), peto_peto = list(label = "Peto-Peto generalized Wilcoxon test", label_key = "test.peto_peto.label", kind = "censored_group_comparison", when = "Compare a censored continuous outcome across 2+ groups when detection limits overlap observed values (not rank-determinable)."), censored_regression = list(label = "Censored regression (Tobit-style)", label_key = "test.censored_regression.label", kind = "regression_censored", when = "Model a censored continuous outcome from predictors without substituting detection limits (lognormal or gaussian error distribution)."), chi_square = list(label = "Chi-squared test of association", label_key = "test.chi_square.label", kind = "association", when = "Association between two categorical variables, adequate expected counts."), fisher_exact = list(label = "Fisher's exact test", label_key = "test.fisher_exact.label", kind = "association", when = "Association between two categorical variables, small expected counts."), pearson = list(label = "Pearson correlation", label_key = "test.pearson.label", kind = "correlation", when = "Linear association between two continuous variables."), spearman = list(label = "Spearman correlation", label_key = "test.spearman.label", kind = "correlation", when = "Monotonic association between two continuous/ordinal variables."), linear_regression = list(label = "Linear regression", label_key = "test.linear_regression.label", kind = "regression_lm", when = "Model a continuous outcome from predictors, adjusting for covariates."), logistic_regression = list(label = "Logistic regression", label_key = "test.logistic_regression.label", kind = "regression_glm", when = "Model a binary outcome from predictors, adjusting for covariates."), firth_logistic = list(label = "Firth penalised logistic regression", label_key = "test.firth_logistic.label", kind = "regression_binary", when = "Model a binary outcome when ordinary logistic regression breaks down: a predictor that perfectly (or almost perfectly) separates the outcome, or very few events per predictor. Firth's penalised likelihood returns FINITE odds ratios with profile-likelihood confidence intervals where the ordinary fit runs off to infinity."), kaplan_meier = list(label = "Kaplan-Meier (log-rank test)", label_key = "test.kaplan_meier.label", kind = "survival_km", when = "Time-to-event outcome (a time column + a 0/1 event indicator): compare survival between groups (KM curves + log-rank test)."), cox_ph = list(label = "Cox proportional-hazards regression", label_key = "test.cox_ph.label", kind = "survival_cox", when = "Time-to-event outcome (time + 0/1 event): hazard ratios for predictors, adjusting for covariates."), aft_weibull = list(label = "Weibull accelerated failure time regression", label_key = "test.aft_weibull.label", kind = "survival_cox", when = "Time-to-event outcome (time + 0/1 event): time ratios for predictors (how much they stretch or shrink expected survival time) - an accelerated failure time model that reports time ratios instead of Cox's hazard ratios and assumes Weibull-distributed event times."), mixed_model = list(label = "Linear mixed-effects model (random intercept)", label_key = "test.mixed_model.label", kind = "mixed_model", label_random_slope = "Linear mixed-effects model (random intercept and slope)", label_random_slope_key = "test.mixed_model.label_random_slope", when = "Continuous outcome with repeated or clustered measurements (multiple rows per subject or site): fixed-effect predictors plus a per-cluster random intercept."), mixed_logistic = list(label = "Mixed-effects logistic regression (random intercept)", label_key = "test.mixed_logistic.label", kind = "mixed_glm_binomial", label_random_slope = "Mixed-effects logistic regression (random intercept and slope)", label_random_slope_key = "test.mixed_logistic.label_random_slope", when = "Binary outcome with repeated or clustered measurements: odds ratios for predictors plus a per-cluster random intercept."), mixed_poisson = list(label = "Mixed-effects Poisson regression (random intercept)", label_key = "test.mixed_poisson.label", kind = "mixed_glm_poisson", label_random_slope = "Mixed-effects Poisson regression (random intercept and slope)", label_random_slope_key = "test.mixed_poisson.label_random_slope", when = "Count outcome with repeated or clustered measurements: rate ratios for predictors plus a per-cluster random intercept."), rm_anova = list(label = "Repeated-measures ANOVA (Greenhouse-Geisser)", label_key = "test.rm_anova.label", kind = "rm_anova", when = "Continuous outcome measured on the SAME subjects across the levels of a within-subject factor (e.g. pre / post / follow-up). Classic repeated-measures ANOVA with a sphericity correction; for unbalanced or missing data prefer the linear mixed model."), mcnemar = list(label = "McNemar's test", label_key = "test.mcnemar.label", kind = "paired_association", when = "Two paired binary measurements on the SAME subjects (e.g. before/after, two raters on one occasion): tests whether the proportion changes."), ancova = list(label = "ANCOVA (analysis of covariance)", label_key = "test.ancova.label", kind = "ancova", when = "Compare a continuous outcome across 2+ groups while adjusting for one or more covariates (e.g. a baseline/pre-treatment value, age, sex)."), poisson_regression = list(label = "Poisson regression", label_key = "test.poisson_regression.label", kind = "regression_count", when = "Model a COUNT outcome (non-negative integers, e.g. number of events) from predictors; assumes the variance equals the mean."), negbin_regression = list(label = "Negative-binomial regression", label_key = "test.negbin_regression.label", kind = "regression_count", when = "Model a COUNT outcome when its variance exceeds its mean (overdispersion), relaxing Poisson's equal-mean-variance assumption."), ordinal_logistic = list(label = "Ordinal logistic regression", label_key = "test.ordinal_logistic.label", kind = "regression_ordinal", when = "Model an ORDERED categorical outcome (3+ ordered levels, e.g. low/medium/high) from predictors via proportional-odds logistic regression."), roc_auc = list(label = "ROC curve / AUC", label_key = "test.roc_auc.label", kind = "roc", when = "How well a CONTINUOUS score/marker discriminates a BINARY gold-standard outcome: the ROC curve, area under the curve (AUC) with a DeLong confidence interval, and the Youden-optimal cutoff with its sensitivity/specificity."), diagnostic_accuracy = list(label = "Diagnostic accuracy (2x2)", label_key = "test.diagnostic_accuracy.label", kind = "diagnostic", when = "Agreement of a BINARY test against a BINARY gold standard: sensitivity, specificity, PPV, and NPV (each with a confidence interval), the 2x2 counts, and disease prevalence."), bland_altman = list(label = "Bland-Altman agreement", label_key = "test.bland_altman.label", kind = "bland_altman", when = "Agreement between two NUMERIC measurements of the same quantity (method comparison): bias (mean difference) and 95% limits of agreement, each with a confidence interval."), roc_compare = list(label = "Compare two ROC curves (DeLong)", label_key = "test.roc_compare.label", kind = "roc_comparison", when = "Which of TWO continuous scores/markers better discriminates the SAME binary gold-standard outcome: each marker's ROC AUC with a DeLong confidence interval, and DeLong's test for the difference between the two correlated (paired) AUCs."), two_way_anova = list(label = "Two-way ANOVA (factorial)", label_key = "test.two_way_anova.label", kind = "two_way_anova", when = "A continuous outcome and TWO categorical factors (2+ levels each): tests each factor's main effect and whether the two factors INTERACT (whether one factor's effect on the outcome depends on the level of the other)."), multinomial_logistic = list(label = "Multinomial logistic regression", label_key = "test.multinomial_logistic.label", kind = "regression_multinomial", when = "Model an UNORDERED categorical outcome with 3+ levels (e.g. treatment A/B/C, no natural order) from predictors, reporting an odds ratio per predictor for each outcome level vs a reference level."), factors_continuous = list(label = "Factors associated with outcome (continuous outcome)", label_key = "test.factors_continuous.label", kind = "factor_screen_lm", when = "Which of 2-30 candidate factors are associated with a CONTINUOUS outcome: the full multivariable model with every candidate is the primary result, with per-factor tests (Benjamini-Hochberg adjusted), correlations and collinearity, and an exploratory LASSO layer reporting how often each factor is selected across 200 bootstrap resamples."), factors_binary = list(label = "Factors associated with outcome (binary outcome)", label_key = "test.factors_binary.label", kind = "factor_screen_binary", when = "Which of 2-30 candidate factors are associated with a BINARY outcome: the full multivariable logistic model with every candidate is the primary result, with per-factor likelihood-ratio tests (Benjamini-Hochberg adjusted), correlations and collinearity, and an exploratory LASSO layer reporting how often each factor is selected across 200 bootstrap resamples."), factors_survival = list(label = "Factors associated with outcome (time-to-event outcome)", label_key = "test.factors_survival.label", kind = "factor_screen_survival", when = "Which of 2-30 candidate factors are associated with a TIME-TO-EVENT outcome (time + 0/1 event): the full Cox model with every candidate is the primary result, with per-factor likelihood-ratio tests (Benjamini-Hochberg adjusted), correlations and collinearity, and an exploratory LASSO layer reporting how often each factor is selected across 200 bootstrap resamples."), table_one_full = list(label = "Descriptive table (Table 1)", label_key = "test.table_one_full.label", kind = "descriptive_table", when = "Describe the whole dataset: every variable summarised on its own terms (median [Q1, Q3] or mean (SD), chosen per variable), counts with confidence intervals for proportions, and a missingness row - optionally split by one grouping variable. No hypothesis is tested."), descriptive_summary = list(label = "Descriptive summary (mean, SD, SE, quartiles)", label_key = "test.descriptive_summary.label", kind = "descriptive_summary", when = "Summarise one or more chosen NUMERIC columns as numbers: n, missing, mean, standard deviation, variance, standard error of the mean with a t-based confidence interval, median, quartiles (R quantile type 7), IQR and range - optionally split by one grouping variable. Each variable uses its own non-missing values. No hypothesis is tested.")) } `.TABLE_ONE_OVERALL_LABEL` <- "Overall" `.TABLE_ONE_OVERALL_COLLISION_LABEL` <- "Overall (all)" `.TABLE_ONE_DESCRIPTIVE_MEAN_SD` <- "mean_sd" `.TABLE_ONE_DESCRIPTIVE_MEDIAN_IQR` <- "median_iqr" `.TABLE_ONE_LEVEL_LABEL_MEAN_SD` <- "mean (SD)" `.TABLE_ONE_LEVEL_LABEL_MEDIAN_Q1Q3` <- "median [Q1, Q3]" `.TABLE_ONE_RANGE_COLUMN` <- "Range (min-max)" `.TABLE_ONE_MISSING_LEVEL` <- "missing, n (%)" `.TABLE_ONE_REF_SUFFIX` <- "_ref" `.TABLE_ONE_REF_LOW` <- "L" `.TABLE_ONE_REF_NORMAL` <- "N" `.TABLE_ONE_REF_HIGH` <- "H" `.TABLE_ONE_REF_LEVELS` <- c("L", "N", "H") `.table_one_ref_flags` <- function (x) toupper(trimws(as.character(x))) `.TABLE_ONE_MIN_STRATA` <- 2L `.TABLE_ONE_MIN_STRATUM_N` <- 2L `.TABLE_ONE_MIN_CATEGORIES` <- 2L `.P_ADJUST_NONE` <- "none" `.TABLE_ONE_REF_ABNORMAL` <- "abnormal" `.TABLE_ONE_REF_NORMAL_CLASS` <- "normal" `.TABLE_ONE_REF_LEVEL_OUTSIDE` <- "outside reference range, n (%)" `.TABLE_ONE_REF_LEVEL_LOW` <- "below reference range (L), n (%)" `.TABLE_ONE_REF_LEVEL_HIGH` <- "above reference range (H), n (%)" `.TABLE_ONE_EXCLUDED_ROLES` <- c("identifier_or_text", "date") `.TABLE_ONE_ROLE_CONTINUOUS` <- "continuous" `.TABLE_ONE_PCT_DIGITS` <- 1L `.table_one_pct` <- function (x) { if (is.null(x) || length(x) == 0 || is.na(x)) return("NA") sprintf(paste0("%.", .TABLE_ONE_PCT_DIGITS, "f"), x) } `.TABLE_ONE_TEST_COLUMN` <- "Test" `.TABLE_ONE_P_COLUMN` <- "p" `.table_one_p_adjusted_column` <- function (method) { sprintf("p (%s)", unname(.P_ADJUST_DISPLAY_NAMES[method]) %||% method) } `table_one_multiplicity_sentence` <- function (k) { sprintf("k = %d baseline comparisons (one per variable); exploratory, hypothesis-generating", as.integer(k)) } `.table_one_role` <- function (x, name) { non_na <- x[!is.na(x)] n_distinct <- length(unique(non_na)) ratio <- if (length(non_na) == 0) 0 else n_distinct/length(non_na) .role_hint(.col_type(x), n_distinct, ratio, name) } `.table_one_ref_companion_for` <- function (df, name) { if (!grepl(paste0(.TABLE_ONE_REF_SUFFIX, "$"), name)) return(NULL) base <- substr(name, 1, nchar(name) - nchar(.TABLE_ONE_REF_SUFFIX)) if (!nzchar(base) || !base %in% names(df) || !is.numeric(df[[base]])) return(NULL) vals <- .table_one_ref_flags(df[[name]][!is.na(df[[name]])]) if (!length(vals)) return(base) if (!all(vals %in% .TABLE_ONE_REF_LEVELS)) return(NULL) base } `.table_one_level_order` <- function (x) { if (is.factor(x)) { obs <- unique(as.character(x[!is.na(x)])) return(levels(x)[levels(x) %in% obs]) } if (is.logical(x)) { lv <- c("FALSE", "TRUE") return(lv[lv %in% unique(as.character(x[!is.na(x)]))]) } unique(as.character(x[!is.na(x)])) } `.table_one_numeric_descriptive` <- function (x) { xx <- x[!is.na(x)] n <- length(xx) if (n < CELEUS$shapiro_min_n || n > CELEUS$shapiro_max_n) return(.TABLE_ONE_DESCRIPTIVE_MEDIAN_IQR) passed <- tryCatch(check_normality(xx)$passed, error = function(e) NA) if (isTRUE(passed)) .TABLE_ONE_DESCRIPTIVE_MEAN_SD else .TABLE_ONE_DESCRIPTIVE_MEDIAN_IQR } `.table_one_num_cell` <- function (x, descriptive) { xx <- x[!is.na(x)] if (!length(xx)) return("NA") if (identical(descriptive, .TABLE_ONE_DESCRIPTIVE_MEAN_SD)) return(sprintf("%s (%s)", .fmt(mean(xx)), .fmt(stats::sd(xx)))) q <- stats::quantile(xx, probs = c(0.25, 0.75), names = FALSE) sprintf(.TABLE_ONE_MEDIAN_Q1Q3_FORMAT, .fmt(stats::median(xx)), .fmt(q[1]), .fmt(q[2])) } `.table_one_count_cell` <- function (k, n, with_interval = FALSE, interval = PROPORTION_INTERVAL_WILSON) { if (is.na(n) || n <= 0) return("NA") pct <- .table_one_pct(100 * k/n) if (!with_interval) return(sprintf("%d (%s%%)", k, pct)) ci <- .proportion_ci(k, n, interval) sprintf("%d (%s%%, 95%% CI %s-%s)", k, pct, .table_one_pct(100 * ci[1]), .table_one_pct(100 * ci[2])) } `.table_one_compare_numeric` <- function (values, g) { ok <- !is.na(values) & !is.na(g) d <- data.frame(value = as.numeric(values[ok]), grp = droplevels(as.factor(g[ok]))) if (nlevels(d$grp) < .TABLE_ONE_MIN_STRATA) return(NULL) sizes <- as.integer(table(d$grp)) if (any(sizes < .TABLE_ONE_MIN_STRATUM_N)) return(NULL) safe_p <- function(expr) tryCatch(suppressWarnings(expr), error = function(e) NA_real_) if (nlevels(d$grp) == 2L) { if (any(sizes < CELEUS$table_one_welch_min_stratum_n)) return(list(test_id = "mann_whitney", p = safe_p(stats::wilcox.test(value ~ grp, data = d)$p.value))) return(list(test_id = "welch_t", p = safe_p(stats::t.test(value ~ grp, data = d)$p.value))) } list(test_id = "kruskal_wallis", p = safe_p(stats::kruskal.test(value ~ grp, data = d)$p.value)) } `.table_one_compare_categories` <- function (x, g) { ok <- !is.na(x) & !is.na(g) if (!any(ok)) return(NULL) tbl <- table(droplevels(as.factor(x[ok])), droplevels(as.factor(g[ok]))) if (nrow(tbl) < .TABLE_ONE_MIN_CATEGORIES || ncol(tbl) < .TABLE_ONE_MIN_STRATA) return(NULL) chi_p <- function() tryCatch(suppressWarnings(stats::chisq.test(tbl)$p.value), error = function(e) NA_real_) minexp <- tryCatch(min(suppressWarnings(stats::chisq.test(tbl)$expected)), error = function(e) NA_real_) if (!is.na(minexp) && minexp >= CELEUS$expected_count_min) return(list(test_id = "chi_square", p = chi_p(), note = "")) exact <- tryCatch(list(p = stats::fisher.test(tbl, workspace = .TABLE_ONE_FISHER_WORKSPACE)$p.value), error = function(e) list(error = conditionMessage(e))) if (!is.null(exact$p)) return(list(test_id = "fisher_exact", p = exact$p, note = "")) if (!grepl(FISHER_WORKSPACE_ERROR_RE, exact$error)) return(NULL) list(test_id = "chi_square", p = chi_p(), note = .TABLE_ONE_FISHER_FALLBACK_NOTE) } `.TABLE_ONE_FISHER_WORKSPACE` <- 200000L `FISHER_WORKSPACE_ERROR_RE` <- "FEXACT|workspace|LDKEY|LDSTP" `.TABLE_ONE_FISHER_FALLBACK_NOTE` <- "exact test not computable for a table this size; chi-squared used" `.TABLE_ONE_MEDIAN_Q1Q3_FORMAT` <- "%s [%s, %s]" `.TABLE_ONE_RANGE_FORMAT` <- "%s to %s" `.is_censor_companion` <- function (name) grepl("^\\.\\.censor\\.\\.", name, perl = TRUE) `.censor_col` <- function (col) paste0("..censor..", col) `.has_censor_companion` <- function (df, col) .censor_col(col) %in% names(df) `.fmt_censored_descriptive` <- function (value, censor, n_below = NULL, n_total = NULL) { qs <- .km_quantiles(value, censor, probs = c(q1 = 0.25, median = 0.5, q3 = 0.75)) heavy <- isTRUE(attr(qs, "heavy")) if (isTRUE(attr(qs, "two_sided"))) { return(list(cell = NA_character_, heavy = heavy, two_sided = TRUE, qs = qs)) } med <- .fmt_km_quantile(qs$median) q1 <- .fmt_km_quantile(qs$q1) q3 <- .fmt_km_quantile(qs$q3) cell <- if (heavy || is.na(med) || med == "NA") { NA_character_ } else { sprintf("%s [%s, %s] (KM)", med, q1, q3) } list(cell = cell, heavy = heavy, two_sided = FALSE, qs = qs, n_below = n_below, n_total = n_total) } `.fmt_km_quantile` <- function (x) { if (is.null(x) || (length(x) == 1L && is.na(x))) return(NA_character_) if (is.character(x)) return(x) format(as.numeric(x), trim = TRUE, scientific = FALSE, digits = 6) } `.km_quantiles` <- function (value, censor, probs = c(0.25, 0.5, 0.75)) { stopifnot(is.numeric(value), length(value) == length(censor)) pn <- names(probs) probs <- as.numeric(probs) names(probs) <- if (!is.null(pn) && length(pn) == length(probs) && any(nzchar(pn))) pn else paste0("p", probs) out <- vector("list", length(probs)) names(out) <- names(probs) ok <- !is.na(value) & !is.na(censor) if (!any(ok)) { for (i in seq_along(out)) out[[i]] <- NA_real_ return(out) } v <- value[ok] c <- as.integer(censor[ok]) n_cens <- sum(c != 0L) pct <- n_cens/length(v) heavy <- pct > CELEUS$censored_pct_heavy left_only <- any(c == -1L) && !any(c == 1L) right_only <- any(c == 1L) && !any(c == -1L) both <- any(c == -1L) && any(c == 1L) if (both) { for (i in seq_along(out)) out[[i]] <- NA_real_ attr(out, "two_sided") <- TRUE attr(out, "heavy") <- heavy return(out) } if (left_only) { fl <- .flip(v, c) sf <- survival::survfit(survival::Surv(fl$time, fl$event) ~ 1) q_flip <- tryCatch({ qq <- stats::quantile(sf, probs = 1 - probs, conf.int = FALSE) as.numeric(if (is.list(qq) && !is.null(qq$quantile)) qq$quantile else qq) }, error = function(e) rep(NA_real_, length(probs))) for (i in seq_along(probs)) { if (heavy && identical(as.numeric(probs[i]), 0.5)) { out[[i]] <- NA_real_ } else if (is.na(q_flip[i])) { L <- max(v[c == -1L]) out[[i]] <- paste0("<", format(L, trim = TRUE, scientific = FALSE, digits = 15)) } else { out[[i]] <- fl$M - q_flip[i] } } } else if (right_only) { sf <- survival::survfit(survival::Surv(v, c == 0L) ~ 1) q <- tryCatch({ qq <- stats::quantile(sf, probs = probs, conf.int = FALSE) as.numeric(if (is.list(qq) && !is.null(qq$quantile)) qq$quantile else qq) }, error = function(e) rep(NA_real_, length(probs))) for (i in seq_along(probs)) { if (heavy && identical(as.numeric(probs[i]), 0.5)) { out[[i]] <- NA_real_ } else if (is.na(q[i])) { L <- min(v[c == 1L]) out[[i]] <- paste0(">", format(L, trim = TRUE, scientific = FALSE, digits = 15)) } else { out[[i]] <- q[i] } } } else { qs <- as.numeric(stats::quantile(v, probs = probs, names = FALSE, type = 7)) for (i in seq_along(probs)) out[[i]] <- qs[i] } attr(out, "heavy") <- heavy attr(out, "pct_censored") <- pct out } `.flip` <- function (value, censor) { ok <- !is.na(value) M <- if (any(ok)) max(value[ok]) else NA_real_ list(time = M - value, event = as.integer(!is.na(censor) & censor == 0L), M = M) } `.fmt_lod_counts` <- function (censor, n_total) { if (n_total <= 0L) return("0 below LOD") n_below <- sum(!is.na(censor) & censor == -1L) n_above <- sum(!is.na(censor) & censor == 1L) below <- sprintf("%d below LOD (%s%%)", n_below, .fmt(100 * n_below/n_total)) if (n_above == 0L) return(below) paste(below, sprintf("%d above LOD (%s%%)", n_above, .fmt(100 * n_above/n_total)), sep = "; ") } `.censored_detected_range` <- function (value, censor) { obs <- !is.na(value) & !is.na(censor) & censor == 0L if (!any(obs)) return(list(min = NA_real_, max = NA_real_)) list(min = min(value[obs]), max = max(value[obs])) } `.censored_two_sided` <- function (value, censor) { ok <- !is.na(value) & !is.na(censor) any(censor[ok] == -1L) && any(censor[ok] == 1L) } `build_descriptive_table` <- function (df, spec) { interval <- spec$interval %||% PROPORTION_INTERVAL_WILSON grp <- spec$group has_group <- !is.null(grp) && length(grp) == 1L && grp %in% names(df) if (has_group) { lv <- .table_one_level_order(df[[grp]]) idx <- lapply(lv, function(l) which(!is.na(df[[grp]]) & as.character(df[[grp]]) == l)) groups <- stats::setNames(idx, lv) } else { groups <- stats::setNames(list(seq_len(nrow(df))), .TABLE_ONE_OVERALL_LABEL) } overall_label <- if (has_group && .TABLE_ONE_OVERALL_LABEL %in% names(groups)) .TABLE_ONE_OVERALL_COLLISION_LABEL else .TABLE_ONE_OVERALL_LABEL ref_of <- stats::setNames(lapply(names(df), function(n) .table_one_ref_companion_for(df, n)), names(df)) ref_cols <- names(df)[!vapply(ref_of, is.null, logical(1))] vars <- setdiff(names(df), c(grp %||% character(0), ref_cols)) vars <- Filter(function(v) !.is_censor_companion(v), vars) vars <- Filter(function(v) !(.table_one_role(df[[v]], v) %in% .TABLE_ONE_EXCLUDED_ROLES), vars) if (!length(vars)) return(NULL) gvec <- if (has_group) as.character(df[[grp]]) else NULL rows <- list() comparisons <- list() add_row <- function(variable, level, cells, range_cell = "", comparison = NULL) { r <- c(variable = variable, level = level, cells) r[[.TABLE_ONE_RANGE_COLUMN]] <- range_cell rows[[length(rows) + 1L]] <<- r comparisons[length(rows)] <<- list(if (has_group) comparison else NULL) } group_cells <- function(f, overall_f) { if (!has_group) return(stats::setNames(overall_f(), .TABLE_ONE_OVERALL_LABEL)) c(vapply(groups, f, character(1)), stats::setNames(overall_f(), overall_label)) } for (v in vars) { col <- df[[v]] if (.has_censor_companion(df, v)) { censor <- df[[.censor_col(v)]] cell_one <- function(ix) { desc <- .fmt_censored_descriptive(col[ix], censor[ix]) n_tot <- sum(!is.na(col[ix]) & !is.na(censor[ix])) lod <- .fmt_lod_counts(censor[ix], n_tot) if (is.na(desc$cell) || !nzchar(desc$cell)) lod else paste(desc$cell, lod, sep = "; ") } add_row(v, "median [Q1, Q3] (KM)", group_cells(cell_one, function() cell_one(seq_len(nrow(df)))), "", NULL) next } if (identical(.table_one_role(col, v), .TABLE_ONE_ROLE_CONTINUOUS)) { conv <- .table_one_numeric_descriptive(col) label <- if (identical(conv, .TABLE_ONE_DESCRIPTIVE_MEAN_SD)) .TABLE_ONE_LEVEL_LABEL_MEAN_SD else .TABLE_ONE_LEVEL_LABEL_MEDIAN_Q1Q3 xx <- col[!is.na(col)] range_cell <- if (length(xx)) sprintf(.TABLE_ONE_RANGE_FORMAT, .fmt(min(xx)), .fmt(max(xx))) else "NA" add_row(v, label, group_cells(function(ix) .table_one_num_cell(col[ix], conv), function() .table_one_num_cell(col, conv)), range_cell, if (has_group) .table_one_compare_numeric(col, gvec)) } else { var_comparison <- if (has_group) .table_one_compare_categories(as.character(col), gvec) first <- TRUE for (l in .table_one_level_order(col)) { local({ lvl <- l add_row(v, lvl, group_cells(function(ix) { sub <- as.character(col[ix]) .table_one_count_cell(sum(sub == lvl, na.rm = TRUE), sum(!is.na(sub))) }, function() { sub <- as.character(col) .table_one_count_cell(sum(sub == lvl, na.rm = TRUE), sum(!is.na(sub)), with_interval = TRUE, interval = interval) }), comparison = if (first) var_comparison) }) first <- FALSE } } rc <- names(ref_of)[vapply(ref_of, function(b) identical(b, v), logical(1))] for (r in rc) { flags <- .table_one_ref_flags(df[[r]]) if (!any(!is.na(flags))) next abnormal_row <- function(level_label, want, comparison = NULL) add_row(v, level_label, group_cells(function(ix) { f <- flags[ix] .table_one_count_cell(sum(f %in% want, na.rm = TRUE), sum(!is.na(f))) }, function() .table_one_count_cell(sum(flags %in% want, na.rm = TRUE), sum(!is.na(flags)), with_interval = TRUE, interval = interval)), comparison = comparison) abnormal <- c(.TABLE_ONE_REF_LOW, .TABLE_ONE_REF_HIGH) ref_comparison <- if (has_group) .table_one_compare_categories(ifelse(is.na(flags), NA, ifelse(flags %in% abnormal, .TABLE_ONE_REF_ABNORMAL, .TABLE_ONE_REF_NORMAL_CLASS)), gvec) abnormal_row(.TABLE_ONE_REF_LEVEL_OUTSIDE, abnormal, ref_comparison) abnormal_row(.TABLE_ONE_REF_LEVEL_LOW, .TABLE_ONE_REF_LOW) abnormal_row(.TABLE_ONE_REF_LEVEL_HIGH, .TABLE_ONE_REF_HIGH) } add_row(v, .TABLE_ONE_MISSING_LEVEL, group_cells(function(ix) .table_one_count_cell(sum(is.na(col[ix])), length(ix)), function() .table_one_count_cell(sum(is.na(col)), length(col)))) } if (!length(rows)) return(NULL) out <- as.data.frame(do.call(rbind, rows), stringsAsFactors = FALSE, check.names = FALSE) rownames(out) <- NULL if (!has_group) return(out) method <- spec$p_adjust %||% CELEUS$default_p_adjust p_raw <- vapply(comparisons, function(c) if (is.null(c)) NA_real_ else as.numeric(c$p), numeric(1)) labels <- vapply(comparisons, function(c) { if (is.null(c)) return("") lbl <- test_catalog()[[c$test_id]]$label %||% c$test_id if (nzchar(c$note %||% "")) sprintf("%s (%s)", lbl, c$note) else lbl }, character(1)) p_cell <- function(p) vapply(p, function(x) if (is.na(x)) "" else .fmt_p(x), character(1)) out[[.TABLE_ONE_TEST_COLUMN]] <- labels out[[.TABLE_ONE_P_COLUMN]] <- p_cell(p_raw) if (!identical(method, .P_ADJUST_NONE)) out[[.table_one_p_adjusted_column(method)]] <- p_cell(.adjust_rows(p_raw, method)) attr(out, "n_comparisons") <- sum(!is.na(p_raw)) attr(out, "multiplicity_note") <- table_one_multiplicity_sentence(sum(!is.na(p_raw))) out } # Complete-case analysis: drop rows missing any analysis value. .acols <- c("smoke") data <- data[stats::complete.cases(data[, .acols, drop = FALSE]), , drop = FALSE] result <- build_descriptive_table(data, list(group = "smoke")) print(result)