---
title: "Risk Difference"
output:
  rmarkdown::html_vignette:
    toc: true
vignette: >
  %\VignetteIndexEntry{Risk Difference}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

## Introduction

In clinical trials, adverse event tables typically display the frequency of events by treatment group. While these counts are informative on their own, regulatory reviewers and study teams often want to see a direct statistical comparison between groups -- specifically, the difference in proportions and a confidence interval around that difference. This quantity is the risk difference.

The risk difference answers a straightforward question: how much more (or less) likely is an event in one group compared to another? A risk difference of 10% with a 95% confidence interval of (2%, 18%) tells you the event rate was 10 percentage points higher in the treatment group, and you can be reasonably confident the true difference falls between 2 and 18 percentage points.

tplyr2 computes risk differences using `stats::prop.test()` with no continuity correction, producing an asymptotic Wald-type confidence interval. This is configured entirely through the `risk_diff` parameter in `layer_settings()`, and the results appear as additional columns in the output alongside the standard count summaries.

Because the proportions being compared are incidence rates, the denominators must be the treatment-arm populations, not just the subjects who had events. Every example below therefore supplies population data via `pop_data()` -- exactly as an adverse event count layer would (see `vignette("count")` and `vignette("denom")`). Without it, `prop.test()` would be handed the wrong denominators and every risk difference would be biased.

## Basic Risk Difference

To add a risk difference to a count layer, pass a `risk_diff` list inside `layer_settings()`. At minimum, you need to specify which two treatment levels to compare via the `comparisons` parameter.

```{r basic}
spec <- tplyr_spec(
  cols = "TRTA",
  pop_data = pop_data(cols = c("TRTA" = "TRT01A")),
  layers = tplyr_layers(
    group_count("AEDECOD",
      settings = layer_settings(
        distinct_by = "USUBJID",
        format_strings = list(
          n_counts = f_str("xxx (xx.x%)", "distinct_n", "distinct_pct")
        ),
        risk_diff = list(
          comparisons = list(c("Xanomeline High Dose", "Placebo"))
        )
      )
    )
  )
)

result <- tplyr_build(spec, tplyr_adae, pop_data = tplyr_adsl)
kable(head(result[, c("rowlabel1", "res1", "res2", "res3", "rdiff1")], 10))
```

The risk difference appears in the `rdiff1` column. Each comparison pair you specify produces one `rdiff` column, numbered sequentially. The column carries a label attribute indicating which groups were compared.

```{r label}
attr(result$rdiff1, "label")
```

Note the order of the comparison pair matters: the first element is the treatment group and the second is the reference. The risk difference is computed as the treatment proportion minus the reference proportion, expressed as a percentage.

## Configuring Comparisons

You are not limited to a single comparison. When your study has multiple active dose groups, you often want to compare each one against placebo. Pass multiple pairs in the `comparisons` list.

```{r multiple}
spec <- tplyr_spec(
  cols = "TRTA",
  pop_data = pop_data(cols = c("TRTA" = "TRT01A")),
  layers = tplyr_layers(
    group_count("AEDECOD",
      settings = layer_settings(
        distinct_by = "USUBJID",
        format_strings = list(
          n_counts = f_str("xxx (xx.x%)", "distinct_n", "distinct_pct")
        ),
        risk_diff = list(
          comparisons = list(
            c("Xanomeline High Dose", "Placebo"),
            c("Xanomeline Low Dose", "Placebo")
          ),
          format = f_str("xx.x (xx.x, xx.x)", "rdiff", "lower", "upper")
        )
      )
    )
  )
)

result <- tplyr_build(spec, tplyr_adae, pop_data = tplyr_adsl)
kable(head(result[, c("rowlabel1", "res1", "res2", "res3", "rdiff1", "rdiff2")], 8))
```

The first comparison (`Xanomeline High Dose` vs `Placebo`) goes into `rdiff1`, and the second (`Xanomeline Low Dose` vs `Placebo`) goes into `rdiff2`. Each column gets its own label.

```{r multiple_labels}
attr(result$rdiff1, "label")
attr(result$rdiff2, "label")
```

## Formatting Output

The risk difference format is controlled through the `format` parameter, which takes an `f_str()` object just like the count format strings. Four variables are available for use in the format:

- `rdiff` -- the risk difference as a percentage
- `lower` -- the lower bound of the confidence interval
- `upper` -- the upper bound of the confidence interval
- `p_value` -- the p-value from `prop.test()`

When no format is specified, tplyr2 uses the default: `f_str("xx.x (xx.x, xx.x)", "rdiff", "lower", "upper")`.

You can include the p-value in the formatted string if your table requires it.

```{r format_pvalue}
spec <- tplyr_spec(
  cols = "TRTA",
  pop_data = pop_data(cols = c("TRTA" = "TRT01A")),
  layers = tplyr_layers(
    group_count("AEDECOD",
      settings = layer_settings(
        distinct_by = "USUBJID",
        format_strings = list(
          n_counts = f_str("xxx (xx.x%)", "distinct_n", "distinct_pct")
        ),
        risk_diff = list(
          comparisons = list(c("Xanomeline High Dose", "Placebo")),
          format = f_str("xx.x (xx.x, xx.x) [x.xxxx]", "rdiff", "lower", "upper", "p_value")
        )
      )
    )
  )
)

result <- tplyr_build(spec, tplyr_adae, pop_data = tplyr_adsl)
kable(head(result[, c("rowlabel1", "res1", "res2", "res3", "rdiff1")], 8))
```

The `x` characters in the format string control field width, just as they do for count format strings. Each `x` reserves one character position, so `xx.x` gives one decimal place with room for a two-digit integer part. You can adjust the precision to match your table shell requirements.

## Confidence Interval Level

By default, tplyr2 computes a 95% confidence interval. You can change this with the `ci` parameter.

```{r ci90}
spec <- tplyr_spec(
  cols = "TRTA",
  pop_data = pop_data(cols = c("TRTA" = "TRT01A")),
  layers = tplyr_layers(
    group_count("AEDECOD",
      settings = layer_settings(
        distinct_by = "USUBJID",
        format_strings = list(
          n_counts = f_str("xxx (xx.x%)", "distinct_n", "distinct_pct")
        ),
        risk_diff = list(
          comparisons = list(c("Xanomeline High Dose", "Placebo")),
          ci = 0.90,
          format = f_str("xx.x (xx.x, xx.x)", "rdiff", "lower", "upper")
        )
      )
    )
  )
)

result <- tplyr_build(spec, tplyr_adae, pop_data = tplyr_adsl)
kable(head(result[, c("rowlabel1", "rdiff1")], 8))
```

A 90% interval is narrower than a 95% interval for the same data. The confidence level you choose should match your study's statistical analysis plan.

## Risk Difference with Distinct Counts

Risk difference calculations naturally work with the `distinct_by` setting. When `distinct_by` is specified, the proportions used for the risk difference are based on distinct subject counts rather than event counts. This is almost always what you want for adverse event tables, where a single subject can contribute multiple events.

```{r distinct}
spec <- tplyr_spec(
  cols = "TRTA",
  pop_data = pop_data(cols = c("TRTA" = "TRT01A")),
  layers = tplyr_layers(
    group_count("AEDECOD",
      settings = layer_settings(
        distinct_by = "USUBJID",
        format_strings = list(
          n_counts = f_str("xxx (xx.x%) [xxx]", "distinct_n", "distinct_pct", "n")
        ),
        risk_diff = list(
          comparisons = list(c("Xanomeline High Dose", "Placebo")),
          format = f_str("xx.x (xx.x, xx.x)", "rdiff", "lower", "upper")
        )
      )
    )
  )
)

result <- tplyr_build(spec, tplyr_adae, pop_data = tplyr_adsl)
kable(head(result[, c("rowlabel1", "res1", "res2", "res3", "rdiff1")], 8))
```

In this output, the count columns show distinct subjects, their percentage, and the total event count in brackets. The risk difference is computed from the distinct subject proportions, which is the clinically meaningful comparison.

## Interaction with Special Rows

An important detail: risk differences are computed before special rows (total rows and missing rows) are appended. This means total and missing rows will have empty risk difference values, which is the expected behavior -- a risk difference on a total row is not meaningful.

```{r special_rows}
spec <- tplyr_spec(
  cols = "TRTA",
  pop_data = pop_data(cols = c("TRTA" = "TRT01A")),
  layers = tplyr_layers(
    group_count("AEDECOD",
      settings = layer_settings(
        distinct_by = "USUBJID",
        format_strings = list(
          n_counts = f_str("xxx (xx.x%)", "distinct_n", "distinct_pct")
        ),
        total_row = TRUE,
        total_row_label = "Any adverse event",
        risk_diff = list(
          comparisons = list(c("Xanomeline High Dose", "Placebo")),
          format = f_str("xx.x (xx.x, xx.x)", "rdiff", "lower", "upper")
        )
      )
    )
  )
)

result <- tplyr_build(spec, tplyr_adae, pop_data = tplyr_adsl)
# Show the last few rows including the total row
tail_rows <- tail(result[, c("rowlabel1", "res1", "res2", "res3", "rdiff1")], 5)
kable(tail_rows)
```

The "Any adverse event" row has an empty `rdiff1` value because it is a total row added after the risk difference computation.

## Extracting Raw Numbers

The formatted risk difference strings are useful for display, but sometimes you need the underlying numeric values for further analysis or custom formatting. The `tplyr_numeric_data()` function gives you access to the raw counts that were used to compute the risk differences.

```{r numeric_data}
spec <- tplyr_spec(
  cols = "TRTA",
  pop_data = pop_data(cols = c("TRTA" = "TRT01A")),
  layers = tplyr_layers(
    group_count("AEDECOD",
      settings = layer_settings(
        distinct_by = "USUBJID",
        format_strings = list(
          n_counts = f_str("xxx (xx.x%)", "distinct_n", "distinct_pct")
        ),
        risk_diff = list(
          comparisons = list(c("Xanomeline High Dose", "Placebo")),
          format = f_str("xx.x (xx.x, xx.x)", "rdiff", "lower", "upper")
        )
      )
    )
  )
)

result <- tplyr_build(spec, tplyr_adae, pop_data = tplyr_adsl)
nd <- tplyr_numeric_data(result, layer = 1)
kable(head(nd, 10))
```

This data.frame contains the raw counts, percentages, and denominators per treatment group and preferred term. These are the values that feed into the `prop.test()` calls.

You can also extract numeric values directly from the formatted risk difference strings using `str_extract_num()`. This function pulls the nth number from a formatted string.

```{r extract_num}
# Extract the risk difference value (1st number)
result$rdiff_value <- str_extract_num(result$rdiff1, 1)

# Extract the lower CI bound (2nd number)
result$rdiff_lower <- str_extract_num(result$rdiff1, 2)

# Extract the upper CI bound (3rd number)
result$rdiff_upper <- str_extract_num(result$rdiff1, 3)

kable(head(result[, c("rowlabel1", "rdiff1", "rdiff_value", "rdiff_lower", "rdiff_upper")], 8))
```

This approach is useful when you need numeric risk difference values for downstream tasks like sorting, filtering, or creating forest plots.

## Risk Difference Versus Association Tests

A risk difference gives an *effect size* -- the difference in incidence with a confidence interval. When you instead (or additionally) want a *p-value* per row from a formal test, use `assoc_test()` in its pairwise mode. It sits in the same place as `risk_diff` (one `pval` column per comparison, a value on every target-level row) but delegates the test to a function you supply -- typically `fisher.test()` for adverse event incidence. See `vignette("binding-statistics")` for the full treatment.

`risk_diff` applies to **single-level** count layers, like the examples above. For a **nested** SOC/PT layer, pairwise `assoc_test()` is the tool that produces a per-row comparison on every level; `vignette("adverse-events")` shows that pattern.

Do not confuse the two confidence-interval controls: `risk_diff`'s `ci` argument sets the coverage of the *between-arm difference* interval shown here, whereas the `ci_method`/`ci_level` layer settings (and the `ci_lower`/`ci_upper` `f_str` keywords) produce a *single-proportion* interval for each individual cell. They are independent features -- see `vignette("denom")` for the single-proportion interval.

## Summary

Risk difference in tplyr2 is configured entirely through the `risk_diff` parameter in `layer_settings()`. The key points to remember:

- **Comparisons** are specified as pairs of treatment levels, with the first element as the treatment and the second as the reference.
- **Multiple comparisons** each produce a separate `rdiff` column (`rdiff1`, `rdiff2`, etc.).
- **Formatting** uses `f_str()` with the variables `rdiff`, `lower`, `upper`, and `p_value`.
- **Confidence level** defaults to 0.95 and is adjustable via the `ci` parameter.
- Risk differences are computed **before** special rows (total, missing), so those rows have empty risk difference values.
- Raw count data is accessible through `tplyr_numeric_data()`, and formatted values can be parsed with `str_extract_num()`.

## Single-Level Layers Only

`risk_diff` works on single-level count layers. Attaching it to a **nested** count
layer -- the `group_count(c("AEBODSYS", "AEDECOD"))` shape most adverse event tables
use -- is a hard error rather than a silent partial result:

```{r nested-error, error = TRUE}
spec <- tplyr_spec(
  cols = "TRTA",
  layers = tplyr_layers(
    group_count(c("AEBODSYS", "AEDECOD"),
      settings = layer_settings(
        distinct_by = "USUBJID",
        risk_diff = list(
          comparisons = list(c("Xanomeline High Dose", "Placebo"))
        )
      )
    )
  )
)

tplyr_build(spec, tplyr_adae)
```

For a per-level comparison column on a nested layer, use a pairwise `assoc_test()`
instead. It emits one `pval` column per comparison with a value on every preferred-term
row *and* every body-system subtotal row, and it accepts any 2x2 test you supply. See
`vignette("binding-statistics")`, which also covers association tests on single-level
and descriptive layers and the pattern for binding externally fitted model results onto
a built table.