---
title: "Comparative Statistics and Binding External Results"
output:
  rmarkdown::html_vignette:
    toc: true
vignette: >
  %\VignetteIndexEntry{Comparative Statistics and Binding External Results}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
library(tplyr2)
library(knitr)
data(tplyr_adsl, package = "tplyr2")
```

# Where statistics come from in tplyr2

tplyr2 is a *grammar for building clinical summary tables*, not a statistical
modeling engine. Every number in a table arrives through one of four routes,
and knowing which to reach for is most of the work:

| You need | Use | Notes |
|---|---|---|
| A descriptive summary (mean, SD, quantiles, ...) | **`group_desc()`** built-in stats | See `vignette("desc")` |
| Counts / incidence, `n (%)` | **`group_count()`** built-in stats | See `vignette("count")` |
| A test or interval computed **across the treatment columns** *from the cell counts* | **`risk_diff`**, **`assoc_test()`**, single-proportion CIs | This vignette |
| A statistic computed by a **model fit** (MMRM, ANCOVA, Cox, logistic, ...) | fit it **externally** and **bind** the formatted result on | This vignette |
| A fully custom in-spec statistic computed *within* one column | **`group_analyze()`** | See `vignette("analyze")` |

The dividing line is simple: if the statistic can be built from the assembled
cell counts (or the raw rows) *across arms*, tplyr2 can compute it natively. If
it requires fitting a model, tplyr2's job is to build the descriptive block and
give you a clean seam to attach the model result. This vignette covers the last
two rows of that table.

# Native comparisons across the treatment columns

These are computed by the layer itself, so the comparison shares one source of
truth with the `n (%)` / summary block beside it.

## Risk difference

`risk_diff` emits pairwise, per-level difference-in-proportion columns
(`rdiff1`, `rdiff2`, ... — one per comparison, a value on every target-level
row) via `stats::prop.test()`. It has its own vignette,
`vignette("riskdiff")`; reach for it when you want an asymptotic risk
difference with a confidence interval.

## Association tests with `assoc_test()`

`assoc_test()` attaches an **arbitrary** test — you supply the function, so
Fisher's exact, a chi-square, or a `coin::cmh_test()` all drop in. It has two
modes.

**Omnibus mode** (the default) runs your function once per `by` group over the
**raw source rows** for that group (all treatment columns at once), and lands a
single value on the group's first output row. It works on **count**, **shift**,
and **desc** layers (a desc layer is shown below). This is the shape for a
per-analyte test that collapses across arms:

```{r assoc-omnibus}
fisher_p <- function(.data) {
  suppressWarnings(fisher.test(table(.data$TRT01P, .data$SEX))$p.value)
}

spec <- tplyr_spec(
  cols = "TRT01P",
  layers = tplyr_layers(
    group_count("SEX",
      settings = layer_settings(
        format_strings = list(n_counts = f_str("xx (xx.x%)", "n", "pct")),
        assoc_test = assoc_test(fn = fisher_p,
                                format = f_str("x.xxx", "p"),
                                label = "Fisher p")))))

b <- tplyr_build(spec, tplyr_adsl)
kable(as_display(b))
```

The p-value sits in the trailing `pval1` column on the first row. Your `fn`
receives the raw rows, so any test that works on a data frame works here.

Omnibus mode also works on a **`group_desc`** layer — the natural home for a
**continuous** comparison across arms (ANOVA, Kruskal-Wallis, a t-test). The
contract is identical (one p per `by` group, on that group's first statistic
row), so a demographics table gets its continuous p-values the same way its
categorical ones do, sharing a single `pval1` column:

```{r assoc-desc}
spec_age <- tplyr_spec(
  cols = "TRT01P",
  layers = tplyr_layers(
    group_desc("AGE",
      settings = layer_settings(
        format_strings = list(Mean = f_str("xx.x", "mean"), SD = f_str("xx.xx", "sd")),
        assoc_test = assoc_test(
          fn = function(.data) anova(lm(AGE ~ TRT01P, .data))[["Pr(>F)"]][1],
          format = f_str("x.xxx", "p"), label = "ANOVA p")))))

kable(as_display(tplyr_build(spec_age, tplyr_adsl)))
```

**Pairwise / per-level mode** (new in 0.2.0, count layers only) is switched on by
supplying `comparisons`. It compares a `reference` arm to each other arm, emits one
`pval<k>` column per comparison with a value on **every** target-level row, and
hands your `fn` a ready-made incidence 2x2 per (level, comparison) built from
the assembled counts *and* population denominators:

```{r assoc-pairwise, eval = FALSE}
group_count("AEDECOD",
  settings = layer_settings(
    distinct_by = "USUBJID",
    stat_columns = list("n" = f_str("xx (xx.x%)", "distinct_n", "distinct_pct")),
    assoc_test = assoc_test(
      fn          = function(m) fisher.test(m)$p.value,  # m = 2x2 count matrix
      reference   = "Placebo",
      comparisons = c("Low", "High"),
      format      = f_str("x.xxx", "p"),
      label       = c("Placebo vs Low", "Placebo vs High"))))
#>   rowlabel1       res1       res2       res3 pval1 pval2
#> 1  HEADACHE  3 (50.0%)  1 (20.0%)  2 (50.0%) 0.524 1.000
#> 2    NAUSEA  2 (33.3%)  3 (60.0%)  1 (25.0%) 0.524 1.000
```

The `fn` receives
`matrix(c(n_ref, n_cmp, N_ref - n_ref, N_cmp - n_cmp), nrow = 2)` — rows are the
(reference, comparison) arm, columns are (event, no-event) — using distinct
counts and denominators when `distinct_by` is set. Because it is just a matrix
in and a scalar out, any 2x2 test (Fisher, chi-square, relative risk, ...)
plugs in.

These snippets show only the count layer, to keep the focus on `assoc_test()`.
In a real adverse-event build the enclosing `tplyr_spec()` sets `pop_data()`, so
the displayed incidence *and* the `N_ref`/`N_cmp` cells of the 2x2 use the safety
population rather than only the subjects who had events (see `vignette("count")`
and `vignette("adverse-events")`, which builds this pattern end to end). The
`stat_columns` shown here is also optional — `format_strings` works identically;
pairwise mode does not require `stat_columns`.

**Nested layers** — the AE-by-SOC/PT case — work the same way: give
`group_count()` a nested target such as `c("AEBODSYS", "AEDECOD")` and the
`pval` columns land on **every** row of **every** level. Each preferred-term row
and each system-organ-class subtotal row gets its own 2x2, built from that row's
counts (the subtotal uses the SOC-level "any event in that class" subject count)
and the population denominators — the standard per-row Fisher column. When the
layer carries a total row, the grand-total ("any event anywhere") p-value is
computed too; pass `assoc_test(..., total_row = FALSE)` to leave it blank.
Missing rows are always blank.

**Returning the finished display string.** `fn` may return a **character**
string instead of a number, and it is dropped into the cell verbatim (`format`
then applies only to numeric returns). The function that computes the test can
therefore also encode the display conventions clinical AE tables use — a
significance flag, a `>.99`/`<.0001` ceiling, trailing-space alignment, or a
`"NE"` sentinel — right where the raw p-value is in hand (`NA`, numeric or
character, still renders a blank):

```{r assoc-char, eval = FALSE}
ae_pval <- function(m) {
  if (sum(m[, 1]) == 0) return(NA_character_)          # both arms zero -> blank
  p <- fisher.test(m)$p.value
  d <- formatC(round(p, 3), format = "f", digits = 3)
  if (p > .99) ">.99" else if (p < .15) paste0(d, "*") else paste0(d, " ")
}
group_count("AEDECOD",
  settings = layer_settings(
    distinct_by = "USUBJID",
    stat_columns = list("n" = f_str("xx (xx.x%)", "distinct_n", "distinct_pct")),
    assoc_test = assoc_test(fn = ae_pval,
                            reference = "Placebo", comparisons = c("Low", "High"))))
```

**Returning several statistics.** A p-value is not the only thing a test produces. When `format` references more than one variable, `fn` returns a numeric *vector* matching it (mapped positionally), so an effect size and its confidence interval land in a single cell. For example, an odds ratio with a 95% CI straight from `fisher.test()`:

```{r assoc-multi, eval = FALSE}
or_ci <- function(m) {
  ft <- fisher.test(m)
  c(ft$estimate, ft$conf.int[1], ft$conf.int[2])   # OR, lower, upper
}
group_count("AEDECOD",
  settings = layer_settings(
    distinct_by = "USUBJID",
    stat_columns = list("n" = f_str("xx (xx.x%)", "distinct_n", "distinct_pct")),
    assoc_test = assoc_test(fn = or_ci,
                            reference = "Placebo", comparisons = c("Low", "High"),
                            format = f_str("xx.xx (xx.xx, xx.xx)", "or", "lo", "hi"),
                            label = "OR (95% CI)")))
```

Each `pval` column then reads like `1.85 (1.10, 3.02)` -- the odds ratio and its interval in one cell. Any statistical procedure that emits a small tuple fits this pattern: an estimate with a p-value, a hazard ratio with an interval, and so on. The values map to the `format` variables positionally, so their names are free; an all-`NA` return (or one whose length does not match `format`) blanks the cell.

## Single-proportion confidence intervals

Count layers can attach an exact or score confidence interval for each cell's
proportion (new in 0.2.0). Set `ci_method` / `ci_level` and reference the
`ci_lower` / `ci_upper` (or `distinct_ci_lower` / `distinct_ci_upper`) keywords
in a format string:

```{r prop-ci, eval = FALSE}
group_count("AEDECOD",
  settings = layer_settings(
    distinct_by = "USUBJID",
    ci_method   = "clopper_pearson",     # also: wilson, wald, agresti_coull, jeffreys
    format_strings = list(
      n_counts = f_str("xx (xx.x%) [xx.x, xx.x]",
                       "distinct_n", "distinct_pct",
                       "distinct_ci_lower", "distinct_ci_upper"))))
#>   rowlabel1                       res1
#> 1  HEADACHE  12 (30.0%) [16.6, 46.5]
```

`clopper_pearson` matches SAS `PROC FREQ ... EXACT` (and `stats::binom.test`);
`wilson` matches `stats::prop.test(correct = FALSE)`.

# Bringing your own statistics

Model-based results — MMRM/ANCOVA p-values, LS-means, hazard ratios, odds
ratios, confidence-interval strings — are computed by a dedicated package
(`emmeans`, `mmrm`, `survival`, `rbmi`, ...) and **bound onto** the assembled
table. The workflow is always the same five steps:

1. **Build** the descriptive block with tplyr2 and pull the display frame with
   `as_display()` (drops the internal `ord_*` columns; the build is already
   ordered).
2. **Compute** the statistic with whatever package you like.
3. **Format** the result into character cells with `apply_formats()` (the same
   `f_str()` spec, rounding, and alignment as the rest of the table).
4. **Align** it to the table — as a new *column* (join on the row label) or a
   new *row* (`rbind`), matching the `rowlabel*` / `res*` shape.
5. **Hand off** to your table-rendering package.

## Attach a statistic as a column

Here we compute a per-level p-value ourselves and bind it beside the `n (%)`
block. (We use base R so the vignette stays dependency-free; swap in your own
model.)

```{r bind-column}
# 1. descriptive block
spec <- tplyr_spec(
  cols = "TRT01P",
  layers = tplyr_layers(
    group_count("SEX",
      settings = layer_settings(
        format_strings = list(n_counts = f_str("xx (xx.x%)", "n", "pct"))))))
disp <- as_display(tplyr_build(spec, tplyr_adsl))

# 2. compute a statistic per row (here: chi-square of SEX x TRT for each SEX
#    level vs the rest). Substitute any model here.
pval <- vapply(disp$rowlabel1, function(lvl) {
  tab <- table(tplyr_adsl$TRT01P, tplyr_adsl$SEX == lvl)
  suppressWarnings(chisq.test(tab)$p.value)
}, numeric(1))

# 3. format with the SAME machinery as the table
disp$pval <- apply_formats(f_str("x.xxx", "p"), pval)

# 4. it is already aligned row-for-row (as_display() preserves build order)
kable(disp)
```

Formatting the external value through `apply_formats()` — rather than
`sprintf()` or `format()` — keeps rounding (including the IBM half-up option)
and width identical to the layer-computed cells.

## Attach a statistic as a row

For a p-value that belongs *under* a descriptive block (e.g. an ANCOVA of a
continuous endpoint), format one cell and `rbind` a labelled row:

```{r bind-row}
spec <- tplyr_spec(
  cols = "TRT01P",
  layers = tplyr_layers(
    group_desc("AGE",
      settings = layer_settings(
        format_strings = list(n = f_str("xx", "n"),
                              "Mean (SD)" = f_str("xx.x (xx.xx)", "mean", "sd"))))))
disp <- as_display(tplyr_build(spec, tplyr_adsl))

# an overall test across arms (substitute lm()/emmeans()/mmrm() as needed)
p <- anova(lm(AGE ~ TRT01P, tplyr_adsl))[["Pr(>F)"]][1]

# build a matching row: label + one p-value cell in the first result column,
# blank in the rest, then bind it on
res_cols <- grep("^res", names(disp), value = TRUE)
p_row <- disp[1, ]                                   # a row of the right shape
p_row$rowlabel1 <- "p-value"
p_row[res_cols] <- ""
p_row[[res_cols[1]]] <- apply_formats(f_str("x.xxx", "p"), p)

out <- rbind(disp, p_row)
rownames(out) <- NULL
kable(out)
```

## Missing or non-estimable results

When a model does not converge or a cell is not estimable, the statistic is
`NA`. `apply_formats()` gains an `na` argument (new in 0.2.0) so a missing value
collapses to a truly empty string rather than a run of spaces — which keeps
`is.na()` / `nzchar()` checks and downstream trimming honest:

```{r na-arg, eval = FALSE}
apply_formats(f_str("xx.x", "x"), c(2.3, NA, 12.7), na = "")
#> [1] " 2.3" ""     "12.7"

apply_formats(f_str("xx.x", "x"), NA_real_, na = "NE")   # or any placeholder
#> [1] "NE"
```

An optional `width =` (with `pad = "right"`/`"left"`) pads the token to a fixed
total width for monospace alignment when the renderer does not own it.

## `group_analyze()` for in-spec custom statistics

`group_analyze()` runs a function you supply *inside the spec*, so a custom
statistic travels with the layer, ordering, and metadata. It is the right tool
for a bespoke *descriptive* statistic (a trimmed mean, a custom responder
definition, a formatted median with a non-standard CI). See
`vignette("analyze")` for the full contract.

One boundary is worth stating here: `analyze_fn` is called **once per
`cols` x `by` combination**, so it only ever sees a *single* treatment column
at a time. It therefore cannot compute a statistic that compares *across* arms
— for that, use `assoc_test()` (which sees all columns) or the external-binding
pattern above.

# Handing off to a table package

`as_display()` returns the display columns only (`rowlabel*`, `res*`, `rdiff*`,
`pval*`), already ordered, ready for `clinify`, `flextable`, `gt`, or
`huxtable`. Pass `labels = TRUE` to rename the result columns to their
column-group headers, and see `vignette("post_processing")` for row masking and
label collapsing.

```{r as-display-labels}
b <- tplyr_build(
  tplyr_spec(cols = "TRT01P",
             layers = tplyr_layers(group_count("SEX",
               settings = layer_settings(
                 format_strings = list(n_counts = f_str("xx (xx.x%)", "n", "pct")))))),
  tplyr_adsl)
kable(as_display(b, labels = TRUE))
```
