---
title: "Creating an exposure-response dataset"
format:
  html:
    toc: true
vignette: >
  %\VignetteIndexEntry{Creating an exposure-response dataset}
  %\VignetteEngine{quarto::html}
  %\VignetteEncoding{UTF-8}
---

```{r}
#| label: setup
#| include: false
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
```

## Introduction

An exposure-response (ER) dataset links a measure of drug **exposure** (for example a model-derived
`AUC` or `CMAX`) to one or more **response** endpoints, so that the relation between the two can be
analysed. Typically the exposure metrics and the response data live in two separate datasets: the
exposure metrics are derived by a modeler from a pharmacokinetic model, while the endpoint data are
prepared by a programmer. Building the ER dataset is the step that brings the two together and gets
everything ready for analysis.

This vignette focuses on the parts of that step that are *specific* to an ER dataset:

- combining the exposure and response data into one dataset, and
- preparing for the analysis (inspecting and lumping sparse categories, imputing
  missing values and flagging outliers).

The parts of the workflow that are the same as for any other dataset -- reading in source data,
attaching meta data through attributes, tabulating and documenting the result, and writing the
analysis file -- are described in the `vignette("getting_started")` and `vignette("example_study")`
vignettes and are therefore not repeated here. To keep this example self-contained we simulate two small
source datasets rather than read them from file.

```{r}
#| message: false
library(dplyr)
library(amp.dm)
```

## Example data

The **exposure** dataset has one row per subject. Besides the exposure metrics (`AUC`, `CMAX`) it
carries the subject-level covariates that the ER analysis may need (`RACE`, `SEX`, `AGE`, `WT`, `ALT`).

```{r}
set.seed(1234)
n <- 30
exposure <- tibble(
  STDY = 1L,
  ID   = 1:n,
  TRT  = rep(1:3, each = 10),                  # dose group
  RACE = sample(rep(1:4, c(14, 10, 4, 2)))     # 4 = a rare category
) |>
  mutate(
    SEX  = rbinom(n(), 1, 0.45),
    AGE  = round(rnorm(n(), 55, 12)),
    WT   = round(rnorm(n(), 78, 14), 1),
    AUC  = round(TRT * rlnorm(n(), log(450), 0.25)),
    CMAX = round(AUC / rnorm(n(), 9, 1.2), 1),
    ALT  = round(runif(n(), 18, 40)),
    ALT  = replace(ALT, ID %in% c(5, 22), c(120, 150)),   # two subjects with elevated ALT
    AGE  = replace(AGE, ID %in% c(15, 16), NA)            # two subjects with a missing age
  ) |>
  select(STDY, ID, TRT, AUC, CMAX, RACE, SEX, AGE, WT, ALT)
head(exposure)
```

The **response** dataset has one row per subject per endpoint. Exposure-response endpoints are often
**binary**, so here `DV` is a 0/1 indicator. Every subject has two endpoints (`EP` 1: an efficacy
response, `EP` 2: a safety event), each with its own response rate.

```{r}
response <- tibble(
  STDY = 1L,
  ID   = rep(1:n, each = 2),
  EP   = rep(1:2, times = n)
) |>
  mutate(DV = rbinom(n(), 1, if_else(EP == 1, 0.5, 0.25)))   # binary response
head(response)
```

## Combining exposure and response

The two datasets are combined with `left_joinr()`. Keeping the response data on the left means every
endpoint record is retained and the matching exposure and covariates are attached to it. `left_joinr()`
is a thin wrapper around `dplyr::left_join()` that additionally **logs** the join, which is worth doing
here: joining exposure onto repeated endpoint records is a classic place to accidentally drop subjects
or create a cartesian product, and the message tells you immediately how many records went in and came
out.

```{r}
er <- left_joinr(response, exposure, by = c("STDY", "ID"),
                 comment = "Combine response with exposure")
head(er)
```

Notice that the exposure metrics and covariates are now **repeated** for each of a subject's endpoint
records (`AUC`, `RACE`, ... are identical on the two rows of every `ID`). This repetition is the one
thing to keep in mind while preparing the covariates below.

## Preparing the covariates

Because each subject contributes more than one row, any summary of a covariate over the combined
dataset counts every subject as many times as it has endpoints. The natural fix is to prepare the
covariates **within each endpoint**: `dplyr::group_by(EP)` splits the data so that each subject appears
exactly once per group. This is also what you want when the endpoints are measured in different subsets
of subjects, so that lumping and imputation follow each endpoint's own data.

### Inspecting the categories

`check_cat()` gives a quick overview of a categorical covariate. Running it on the per-subject
`exposure` data shows the true subject counts and highlights that races `3` and `4` are sparse.

```{r}
check_cat(exposure$RACE)
```

### Lumping sparse categories

Sparse categories often cannot support a covariate effect and are usually lumped together. `num_lump()`
does this; a category is lumped when it stays below the `min` count (or `prop` proportion) you specify.

Running it over the whole combined dataset is misleading, because every subject is counted once per
endpoint: race `3` has 4 subjects but 8 records, so with `min = 5` it is *not* lumped:

```{r}
table(num_lump(er$RACE, lumpcat = 99, min = 5))
```

Lumping **within each endpoint** solves this. Grouped by `EP` each subject appears once, so `min` acts
on the subject counts and both sparse races move into the new category `99`:

```{r}
er <- er |>
  group_by(EP) |>
  mutate(RACE_L = num_lump(RACE, lumpcat = 99, min = 5)) |>
  ungroup()
table(er$RACE_L)
```

`lumpcat` controls *where* the lumped subjects go. Instead of a new category you can send them to the
largest existing group, or map specific categories yourself -- again within each endpoint:

```{r}
# move the sparse races into the largest group
er2 <- er |> group_by(EP) |> mutate(RACE_L = num_lump(RACE, lumpcat = "largest", min = 5))
table(er2$RACE_L)
# or map specific categories by hand (3 and 4 both become 2)
er2 <- er |> group_by(EP) |> mutate(RACE_L = num_lump(RACE, lumpcat = c("3" = 2, "4" = 2), min = 5))
table(er2$RACE_L)
```

### Imputing missing covariates

Missing covariate values can be filled in with `impute_covar()`, which replaces `NA` with a statistic
(the median by default, or the largest group for a categorical covariate). Done within each endpoint,
the statistic is again computed over subjects rather than records. The two subjects whose age was
missing receive the median age:

```{r}
filter(er, ID%in%15:16) |> pull(AGE)            # missing before
er <- er |> group_by(EP) |> mutate(AGE = impute_covar(AGE, type = "median")) |> ungroup()
filter(er, ID%in%15:16) |> pull(AGE)            # imputed after
```

### Flagging outliers

Finally, `flag_outliers()` marks records whose value is a boxplot outlier. This does not remove
anything; it adds a `0/1` flag you can use later to review or exclude the points. The two subjects with
an elevated `ALT` are picked up (two endpoint records each, so four rows):

```{r}
er$ALT_OUT <- flag_outliers(er$ALT)
table(er$ALT_OUT)
```

## The analysis-ready dataset

The result is a single dataset with one row per subject and endpoint, each response value sitting next
to the subject's exposure and cleaned-up covariates -- ready for the exposure-response analysis.

```{r}
head(er)
```

From here the dataset is finished like any other: attach meta data with `attr_xls()` / `attr_add()`,
review it with `define_tbl()`, `stats_df()` and `counts_df()`, and write it out with `output_data()`.
Those steps are covered in the `vignette("getting_started")` and `vignette("example_study")` vignettes.