Evaluating PCA on public data

Andrew Stein

This article runs synpmx_pca() over the same eight public datasets that Evaluating AVATAR on public data runs synpmx_avatar() over, so the two generators can be read against each other on identical studies.

The eight split three ways, and the split is the finding rather than a presentational choice. synpmx_pca() builds every feature on a declared nominal grid and refuses without one, so what separates these datasets is where that grid comes from. Two ship it: case1_pkpd and mad carry a NOMTIME column. Five can supply it — warfarin and wbcSim because their recorded observation times already are the protocol’s, and theo_md, nimoData and mavoglurant because the protocol can be written down from the dose times, the declared occasion, or an occasion-relative clock. One cannot: pheno_sd is routine neonatal care with no protocol behind it, and the generator refuses.

Seven runs follow, then three cross-dataset tables and one finding — the copy check that fails on two of the nine, and what it does and does not mean.

Plotting and reporting helpers used throughout this vignette

# Observation rows only, source beside synthetic on a shared y axis, one row per
# endpoint. The same shape as the AVATAR vignette's figure, reduced to what the
# PCA output needs: this generator writes onto a nominal grid, so there is no
# time-after-dose view to switch to.
# `endpoints` draws one endpoint of a study that has several. A log axis is the
# reason it exists: `facet_grid()` puts one scale on every panel, so a study
# whose concentration wants a log axis and whose PD endpoint goes negative
# cannot be one figure.
overlay_plot <- function(source, synthetic, roles, title,
                         y_label = "DV", log_y = FALSE, alpha = 0.3,
                         max_time = Inf, endpoints = NULL) {
  frame <- function(data, label) {
    observed <- as.character(data[[roles$evid]]) %in% c("0", "0.0") &
      !is.na(data[[roles$dv]])
    if (!is.null(roles$mdv)) {
      observed <- observed & as.character(data[[roles$mdv]]) %in% c("0", "0.0")
    }
    observed <- observed & as.numeric(data[[roles$time]]) <= max_time
    data.frame(
      dataset = factor(label, levels = c("Source", "Synthetic")),
      subject = as.character(data[[roles$id]][observed]),
      # A study whose clock resets inside an occasion needs one line per
      # occasion; joining them draws a line from the end of the first period
      # back to the start of the second.
      occasion = if (is.null(roles$occasion)) "1" else
        as.character(data[[roles$occasion]][observed]),
      time = as.numeric(data[[roles$time]][observed]),
      dv = as.numeric(data[[roles$dv]][observed]),
      endpoint = if (is.null(roles$dvid)) "DV" else
        as.character(data[[roles$dvid]][observed]),
      stringsAsFactors = FALSE
    )
  }
  plotted <- rbind(frame(source, "Source"), frame(synthetic, "Synthetic"))
  if (!is.null(endpoints)) {
    plotted <- plotted[plotted$endpoint %in% endpoints, , drop = FALSE]
  }
  figure <- ggplot2::ggplot(
    plotted,
    ggplot2::aes(time, dv,
                 group = interaction(dataset, subject, occasion),
                 colour = dataset)
  ) +
    ggplot2::geom_line(alpha = alpha) +
    ggplot2::geom_point(alpha = alpha, size = 0.7) +
    # Endpoint down the side, source and synthetic across, so each endpoint has
    # its own y scale while source and synthetic share one. The transpose lets
    # the two panels drift onto different axes and squeezes every endpoint onto
    # a single scale, which is the arrangement the eye cannot compare.
    (if (length(unique(plotted$endpoint)) > 1L) {
      ggplot2::facet_grid(endpoint ~ dataset, scales = "free_y", switch = "y")
    } else {
      ggplot2::facet_wrap(~dataset)
    }) +
    ggplot2::scale_colour_manual(values = comparison_colours) +
    ggplot2::labs(x = "Time (hours)", y = y_label, title = title) +
    ggplot2::theme_minimal() +
    ggplot2::theme(legend.position = "none", strip.placement = "outside")
  if (isTRUE(log_y)) figure + ggplot2::scale_y_log10() else figure
}

overlay_height <- function(data, roles, per_endpoint = 2.2, minimum = 3.4) {
  observed <- as.character(data[[roles$evid]]) %in% c("0", "0.0")
  endpoints <- if (is.null(roles$dvid)) 1L else
    length(unique(data[[roles$dvid]][observed]))
  max(minimum, per_endpoint * endpoints)
}

# Every dataset is run the same way and its numbers collected in one place, so
# the cross-dataset tables at the end cannot drift from the sections above them.
runs <- list()
pca_run <- function(label, source, roles, seed) {
  summary <- synpmx_pca_summarize(source, roles, seed = seed)
  synthetic <- synpmx_pca_generate(summary, seed = seed)
  runs[[label]] <<- list(label = label, source = source, roles = roles,
                         summary = summary, synthetic = synthetic,
                         card = synpmx_scorecard(source, synthetic, roles))
  invisible(runs[[label]])
}

# Nearest-neighbour snapping onto a stated design grid. Used where a study's
# recorded times have to be placed on the protocol's planned ones; which grid to
# snap to is a statement about the study and is written out at each call.
snap_to <- function(x, grid) {
  grid[max.col(-abs(outer(x, grid, "-")), ties.method = "first")]
}

The datasets

Every dataset here is public package data from nlmixr2data or xgxr, both Suggests of this package. The last column is what this vignette is about.

Dataset	Package	Patients	Endpoints	Where its nominal grid comes from
case1_pkpd	xgxr	180	2	Declared. Ships a NOMTIME column
mad	xgxr	60	5	Declared. Ships a NOMTIME column
warfarin	nlmixr2data	32	2	Already on it. Its 16 observation times are the protocol’s
wbcSim	nlmixr2data	45	1	Already on it, apart from a few extended-follow-up patients
theo_md	nlmixr2data	12	1	Constructed from the Q24H dose times and a design grid
nimoData	nlmixr2data	12	1	Constructed from the declared occasion and time after dose
mavoglurant	nlmixr2data	120	1	Constructed by snapping an occasion-relative clock
pheno_sd	nlmixr2data	59	1	None. Routine care, individualised dosing — refused
mixroute_sim	synpmx	90	1	Declared. Ships an NTIME column, and an ADM for the route
onc_sim	synpmx	200	2	Declared. Scans and troughs are planned, so NTIME is the recorded day

The last two are simulated by synpmx rather than public package data, and carry what the eight above do not: a declared administration column, and a daily regimen written with ADDL/II. See ?mixroute_sim and ?onc_sim.

Shared workflow

Every example follows the same five steps:

1. declare the column meanings with pmx_roles(), including nominal_time;
2. summarize with synpmx_pca_summarize() and generate with synpmx_pca_generate();
3. plot the real and synthetic data;
4. plot the observation and baseline covariate distributions, source against synthetic;
5. report the scorecard.

Eight scorecard rows read not applicable on every card below, which the cards count in their tally and do not list. B1a, B1b and C2 read a run record that synpmx_avatar() writes and this generator does not. B4a asks whether a generated set of observation times copies a real one, and B2 whether a synthetic patient stands out from its stratum; both are disclosure questions where a real patient’s record was reused, and this generator draws each visit and each profile from a model of the cohort instead. E1, E2 and E3 score a fitted population model, and this generator fits none. That is a limitation of those rows rather than a property of any dataset here, and it is the same on every card. B4b asks the copy question directly on the finished tables and is the row that carries the claim; the section at the end reports what it found.

The nominal grid is the requirement

synpmx_pca() places dose rows and observation rows on the same declared grid, and every feature it models is a cell on it. A grid inferred from elapsed time puts a sample on the wrong side of a dose as soon as doses were recorded as actuals, and deciding where a trough belongs is a statement about the protocol rather than something an algorithm can derive. So the function refuses without nominal_time instead of guessing:

data("pheno_sd", package = "nlmixr2data")
pheno_roles <- pmx_roles(
  id = "ID", time = "TIME", dv = "DV", amt = "AMT", evid = "EVID",
  mdv = "MDV", covariates = c("WT", "APGR")
)
synpmx_pca_summarize(pheno_sd, pheno_roles, seed = 111)
#> Error:
#> ! `synpmx_pca()` requires `nominal_time` in `pmx_roles()`.
#>   The nominal grid is the axis every feature sits on, and inferring it from
#>   recorded times is a statement about the protocol that only you can make.
#>   Fix: Add the protocol's planned times as a column and declare it.

This is the sharpest difference between the two generators. synpmx_avatar() derives a grid from the recorded times when none is declared, and on pheno_sd that derivation does real work. synpmx_pca() has no such fallback, and the nine runs below are the nine studies where a grid could honestly be supplied.

case1_pkpd: a declared nominal grid

180 patients across six treatment arms, two endpoints keyed by a character NAME column, a baseline weight and a CENS column. The industry shape, and the one the generator is easiest on.

case1_pkpd <- as.data.frame(
  get(utils::data(list = "case1_pkpd", package = "xgxr"))
)
# That study's CENS is meaningful only for PK: the PD effect is signed, so a
# left-censored PD row would report a value above the uncensored ones.
case1_pkpd$CENS <- ifelse(case1_pkpd$NAME == "PD - Continuous", 0,
                          case1_pkpd$CENS)
case1_roles <- pmx_roles(
  id = "ID", time = "TIME", dv = "LIDV", cens = "CENS", amt = "AMT",
  evid = "EVID", cmt = "CMT", dvid = "NAME", nominal_time = "NOMTIME",
  strata = c("TRTACT", "DOSE"), covariates = "WEIGHTB", keep = "STUDY"
)
case1 <- pca_run("case1_pkpd", case1_pkpd, case1_roles, seed = 808)

compare_pmx_distributions(case1$source, case1$synthetic, case1_roles)

synpmx_scorecard_datatable(case1$card, report = "minimal")

Nothing fails. The figure above is on a linear axis over the whole study, which is the structural view: both endpoints present, on the same ranges, over the same 2040 hours. It is the wrong axis for reading a PK concentration, and the censoring behaviour it cannot show is worth a number instead — the fraction of PK observations reported at the assay limit, per arm:

blq_share <- function(data, label) {
  rows <- data$EVID == 0 & data$NAME == "PK Concentration" & !is.na(data$LIDV)
  out <- aggregate(list(pct = data$CENS[rows] == 1),
                   list(arm = data$TRTACT[rows]),
                   function(x) round(100 * mean(x)))
  stats::setNames(out, c("arm", label))
}
knitr::kable(
  merge(blq_share(case1$source, "source_pct_blq"),
        blq_share(case1$synthetic, "synthetic_pct_blq"), by = "arm"),
  caption = "Percentage of PK observations below the limit of quantification."
)

Percentage of PK observations below the limit of quantification.

arm	source_pct_blq	synthetic_pct_blq
10 mg	70	67
100 mg	11	11
3 mg	95	92
30 mg	50	40
300 mg	6	4

The censored share tracks the source arm by arm — within a few points on four of the five arms, and 50% against 39% on the 30 mg arm, which is the one sitting nearest the limit and so the one a small shift in the drawn values moves most. That is the assay limit being reapplied at emit: the drawn value is the latent one, and where the source declared a lower limit of quantification the reported value and CENS are rebuilt from it together. vignette("pca-demo") shows the same thing as a log-scale figure over the day-1 and last-dose windows.

A5a reads review, at 30.7 -> 29 observations per patient, source against synthetic. Every patient was seen at every visit the grid models, so both endpoints’ attendance rates are 1 and the synthetic count is the grid’s. The source also holds two pre-dose concentration records per treated patient, at -24 and -0.1 hours, which carry no value and are not modelled.

mad: five endpoints, including ordinal, count and binary

60 subjects in a multiple-ascending-dose study with a declared NOMTIME and five endpoints — PK concentration, continuous PD, and ordinal, count and binary PD. An endpoint silently vanishing during generation is invisible with fewer than three, which is what this dataset is here to catch.

mad <- as.data.frame(get(utils::data(list = "mad", package = "xgxr")))
mad_roles <- pmx_roles(
  id = "ID", time = "TIME", dv = "LIDV", amt = "AMT", evid = "EVID",
  cmt = "CMT", dvid = "NAME", mdv = "MDV", nominal_time = "NOMTIME",
  strata = c("TRTACT", "DOSE"), covariates = c("WEIGHTB", "SEX")
)
mad_run <- pca_run("mad", mad, mad_roles, seed = 909)

synpmx_scorecard_datatable(mad_run$card, report = "minimal")

Nothing fails, and A3 reads 5 of 5: every endpoint survives. A6 is the check that the discrete endpoints kept their scale, and it reads the finished table rather than trusting the mechanism.

discrete <- c("PD - Binary", "PD - Count", "PD - Ordinal")
values_taken <- function(data, endpoint) {
  values <- data$LIDV[data$NAME == endpoint & data$EVID == 0]
  values <- values[!is.na(values)]
  sprintf("%d distinct, %.2f to %.2f", length(unique(values)),
          min(values), max(values))
}
knitr::kable(
  data.frame(
    endpoint = discrete,
    source = vapply(discrete, values_taken, character(1), data = mad),
    synthetic = vapply(discrete, values_taken, character(1),
                       data = mad_run$synthetic),
    row.names = NULL
  ),
  caption = "Values taken by the three discrete endpoints."
)

Values taken by the three discrete endpoints.

endpoint	source	synthetic
PD - Binary	2 distinct, 0.00 to 1.00	2 distinct, 0.00 to 1.00
PD - Count	20 distinct, 0.00 to 19.00	25 distinct, 0.00 to 42.00
PD - Ordinal	3 distinct, 1.00 to 3.00	3 distinct, 1.00 to 3.00

A drawn feature is a real number, so without the snap a binary endpoint would come back as a spread of values around zero and one. Step 7 puts each drawn value back onto the source’s level set, and the cost is stated rather than hidden: a generated value for a discrete endpoint is a source value, because a 0/1 endpoint has no third value to emit.

warfarin: the observation times are already the protocol’s

32 subjects, a single dose, a PK endpoint (cp) and a PD one (pca) on different time courses. Its 515 rows carry only 16 distinct observation times — 0.5, 1, 1.5, 2, 3, 6, 9, 12, 24, 36, 48, 72, 96, 120, 144 and a baseline — which is a protocol grid already written down.

data("warfarin", package = "nlmixr2data")
warfarin <- as.data.frame(warfarin)
sort(unique(warfarin$time[warfarin$evid == 0]))
#>  [1]   0.0   0.5   1.0   1.5   2.0   3.0   6.0   9.0  12.0  24.0  36.0  48.0
#> [13]  72.0  96.0 120.0 144.0

Declaring nominal_time is then a statement that these recorded times are the planned ones, which for this study is true. No arithmetic is involved, and that is the point: the column is a declaration, not a derivation.

warfarin$ntime <- warfarin$time
warfarin_roles <- pmx_roles(
  id = "id", time = "time", nominal_time = "ntime", dv = "dv", amt = "amt",
  evid = "evid", dvid = "dvid", covariates = c("wt", "age", "sex")
)
warfarin_run <- pca_run("warfarin", warfarin, warfarin_roles, seed = 404)

compare_pmx_distributions(warfarin_run$source, warfarin_run$synthetic,
                          warfarin_roles)

synpmx_scorecard_datatable(warfarin_run$card, report = "minimal")

wbcSim: infusions and a delayed response

45 subjects with infusion start/stop pairs and a delayed white-blood-cell decline, nadir and recovery. Its observation times are a grid too, apart from a few patients followed far longer than the rest.

data("wbcSim", package = "nlmixr2data")
wbcSim <- as.data.frame(wbcSim)
wbcSim$NTIME <- wbcSim$TIME
wbc_roles <- pmx_roles(
  id = "ID", time = "TIME", nominal_time = "NTIME", dv = "DV", amt = "AMT",
  evid = "EVID", cmt = "CMT", rate = "RATE"
)
wbc_run <- pca_run("wbcSim", wbcSim, wbc_roles, seed = 505)

synpmx_scorecard_datatable(wbc_run$card, report = "minimal")

The figure is cut at 700 h because a handful of source patients are followed to 4580 h and the rest are not. Those late visits are exactly what the column guard removes: a grid cell fewer than min_column_patients patients reached describes those patients and nobody else, so it is dropped rather than modelled.

features <- as.data.frame(pca_features(wbc_run$summary))
c(source_times = length(unique(wbcSim$TIME[wbcSim$EVID == 0])),
  modelled_cells = sum(features$kind == "endpoint_cell"),
  last_modelled = max(features$time, na.rm = TRUE))
#>   source_times modelled_cells  last_modelled 
#>             38             21            648

A5a and A5b read 3.91 -> 3.11 and 1.16 -> 1 for the same reason: dropping the long tail shortens both the observation record and the dosing. This is the cost of the guard rather than a defect, and it is the number to judge — a study whose late follow-up matters is one to run at a lower min_column_patients and re-read this row.

theo_md: seven oral doses on a constructible grid

12 subjects, seven doses exactly 24 h apart, with dense concentration sampling around the first and last dose. The doses are on a clean grid and the samples are not, so the grid can be written down: the dose interval a sample falls in, plus its time after that dose snapped to the design times.

data("theo_md", package = "nlmixr2data")
theo_md <- as.data.frame(theo_md)
sort(unique(theo_md$TIME[theo_md$EVID != 0]))
#> [1]   0  24  48  72  96 120 144

theo_doses <- seq(0, 144, by = 24)          # the protocol's Q24H schedule
theo_design <- c(0, 0.25, 0.5, 1, 2, 3, 4, 5, 7, 9, 12, 24)  # planned samples
occasion <- pmax(1L, findInterval(theo_md$TIME, theo_doses))
theo_md$NTIME <- ifelse(
  theo_md$EVID == 0,
  theo_doses[occasion] +
    snap_to(theo_md$TIME - theo_doses[occasion], theo_design),
  theo_md$TIME
)
c(recorded_times = length(unique(theo_md$TIME[theo_md$EVID == 0])),
  nominal_times = length(unique(theo_md$NTIME[theo_md$EVID == 0])))
#> recorded_times  nominal_times 
#>            156             26

156 recorded times become 26 nominal ones. The design grid above is a statement about the study — it says which samples the protocol planned — and it is written in the dataset where a reader can check it rather than inside a function they cannot.

theo_roles <- pmx_roles(
  id = "ID", time = "TIME", nominal_time = "NTIME", dv = "DV", amt = "AMT",
  evid = "EVID", cmt = "CMT", covariates = "WT"
)
theo_run <- pca_run("theo_md", theo_md, theo_roles, seed = 303)

synpmx_scorecard_datatable(theo_run$card, report = "minimal")

Nothing fails, on twelve subjects. The cap holds the basis to two components here — a fifth of the cohort — which is what stops a basis fitted on that few subjects from interpolating them.

c(subjects = theo_run$summary$n_source,
  components = theo_run$summary$basis$k,
  features = length(theo_run$summary$basis$columns))
#>   subjects components   features 
#>         12          2         27

Two components against 27 features is the honest reading of twelve patients, and it is why the synthetic profiles are visibly smoother than the source ones.

nimoData: ten infusions recorded at the times they happened

12 subjects, ten roughly weekly infusions, with the dose times recorded as actuals — 165.70, 167.24 and 168.05 for what the protocol called one weekly visit. The dataset declares an occasion and a time after dose, so the protocol can be written down from what it records. This is the construction worked out in the AVATAR evaluation, reused unchanged.

data("nimoData", package = "nlmixr2data")
nimoData <- as.data.frame(nimoData)
nimo_interval <- 168                    # the protocol's weekly infusion
last_occasion <- ave(nimoData$OCC, nimoData$ID, FUN = max)

nominal_tad <- round(nimoData$TAD / 24) * 24
# A sample at the end of its interval is the NEXT infusion's pre-dose trough,
# not that infusion's time-zero sample. Only where another infusion follows:
# after the last one, a long time after dose is real follow-up.
pre_dose <- nimoData$EVID == 0 & nominal_tad >= nimo_interval &
  nimoData$OCC < last_occasion
nominal_tad[pre_dose] <- nimo_interval - 1

nimoData$NTIME <- (nimoData$OCC - 1) * nimo_interval +
  ifelse(nimoData$EVID == 0, nominal_tad, 0)

The trough is the decision the arithmetic will not make. Rounding time after dose to the nearest day would send a sample drawn just before the next infusion onto that infusion’s own nominal time, where it stops being a trough and becomes a time-zero sample; the line above gives it a slot of its own instead.

nimo_roles <- pmx_roles(
  id = "ID", time = "TIME", nominal_time = "NTIME", dv = "DV", amt = "AMT",
  evid = "EVID", rate = "RATE", mdv = "MDV", tad = "TAD", occasion = "OCC",
  covariates = c("BSA", "AGE", "HGT"), keep = "DOS"
)
nimo_run <- pca_run("nimoData", nimoData, nimo_roles, seed = 606)

#> Warning in transformation$transform(x): NaNs produced
#> Warning in ggplot2::scale_y_log10(): log-10 transformation
#> introduced infinite values.
#> Warning in transformation$transform(x): NaNs produced
#> Warning in ggplot2::scale_y_log10(): log-10 transformation
#> introduced infinite values.
#> Warning: Removed 1 row containing missing values or values outside the scale range
#> (`geom_line()`).
#> Warning: Removed 1 row containing missing values or values outside the scale range
#> (`geom_point()`).

synpmx_scorecard_datatable(nimo_run$card, report = "minimal")

The dosing comes through in full, which is the contrast worth drawing with the AVATAR run on this study. There, reaching the B1b guarantee cost the dosing: 10 doses per patient became 1.6, and every avatar’s course was truncated. Here the doses are the arm’s planned schedule and every generated patient receives all ten of them.

knitr::kable(pca_dosing(nimo_run$summary),
             caption = "The planned schedule, which every generated patient receives.",
             digits = 1)

The planned schedule, which every generated patient receives.

arm	cycle	time	planned_amt
all	1	0	100
all	2	168	100
all	3	336	100
all	4	504	100
all	5	672	100
all	6	840	100
all	7	1008	100
all	8	1176	100
all	9	1344	100
all	10	1512	100

The cost is on the other side and is stated in vignette("pca-algorithm"): there is one schedule where the source had twelve. C2 is the row that would say so, and it is one of the rows that cannot read this generator’s output. pca_dosing() above is what to read in its place.

mavoglurant: an occasion-reset clock

120 subjects in one- and two-period profiles, with TIME resetting inside OCC so it is already dose-relative. There is nothing to align, only a grid to snap to: 117 distinct recorded times across the two periods against a protocol that planned a dozen or so.

data("mavoglurant", package = "nlmixr2data")
mavoglurant <- as.data.frame(mavoglurant)
mavo_design <- c(0, 0.25, 0.5, 0.75, 1, 1.5, 2, 3, 4, 6, 8, 10, 12, 24, 36, 48)
mavoglurant$NTIME <- ifelse(mavoglurant$EVID == 0,
                            snap_to(mavoglurant$TIME, mavo_design),
                            mavoglurant$TIME)
mavo_roles <- pmx_roles(
  id = "ID", time = "TIME", nominal_time = "NTIME", dv = "DV", amt = "AMT",
  evid = "EVID", cmt = "CMT", rate = "RATE", mdv = "MDV", occasion = "OCC",
  keep = "DOSE", covariates = c("AGE", "SEX", "WT", "HT")
)
mavo_run <- pca_run("mavoglurant", mavoglurant, mavo_roles, seed = 707)

synpmx_scorecard_datatable(mavo_run$card, report = "minimal")

This is the most expensive construction in the vignette, and A5a says so: 20.2 -> 11 observations per patient, source against synthetic. Snapping 117 distinct recorded times onto 15 nominal ones collapses several real samples into one cell, and a cell holds one value per subject.

c(source_rows = nrow(mavoglurant),
  synthetic_rows = nrow(mavo_run$synthetic),
  recorded_times = length(unique(mavoglurant$TIME[mavoglurant$EVID == 0])),
  nominal_times = length(unique(mavoglurant$NTIME[mavoglurant$EVID == 0])))
#>    source_rows synthetic_rows recorded_times  nominal_times 
#>           2678           1440            117             15

That is a property of the grid supplied, not of the generator: a wider design grid keeps more of the samples, at the price of cells fewer patients share. It is the one dataset here where the construction is worth revisiting before using the output, and A5a is where a reader would notice.

pheno_sd: no protocol to write down

59 neonates on phenobarbital in routine care, with individualised dosing and sparse irregular sampling: 155 observations at 118 distinct times. There is no protocol grid because there was no protocol, and the refusal at the top of this vignette is the whole of this section.

c(observations = sum(pheno_sd$EVID == 0),
  distinct_times = length(unique(pheno_sd$TIME[pheno_sd$EVID == 0])),
  patients = length(unique(pheno_sd$ID)))
#>   observations distinct_times       patients 
#>            155            118             59

Snapping these onto a day grid would produce a column, and it would be an invention rather than a declaration — nobody planned those visits, so there are no planned times to recover. synpmx_avatar() is the generator for this study; its derived grid does real work there, and the AVATAR evaluation runs it.

mixroute_sim: a grid that is already the protocol’s

Ninety patients in three arms — intravenous, subcutaneous, and intravenous then subcutaneous. Its NTIME is the protocol’s planned time and ships with the dataset, so nothing has to be constructed. What it is here to test is the route: PCA does not model a route, and until it carried one, a declared ADM was dropped from the output entirely and the generated table could not be validated at all.

mixroute_roles <- pmx_roles(
  id = "ID", time = "TIME", nominal_time = "NTIME", dv = "DV", amt = "AMT",
  evid = "EVID", cmt = "CMT", adm = "ADM",
  routes = c("1" = "iv", "2" = "extravascular"),
  cens = "CENS", strata = "ARM", covariates = "WT"
)
mixroute_run <- pca_run("mixroute_sim", mixroute_sim, mixroute_roles, seed = 808)

synpmx_scorecard_datatable(mixroute_run$card, report = "minimal")

Nothing fails. The route comes back because each generated dose enters the compartment its own route doses into — .draw_schedule() carries the route of every planned dose, and the administration id is written from it:

dosed <- function(data) data[data$EVID == 1L, ]
list(source = table(dosed(mixroute_sim)$ARM, dosed(mixroute_sim)$ADM),
     synthetic = table(dosed(mixroute_run$synthetic)$ARM,
                       dosed(mixroute_run$synthetic)$ADM))
#> $source
#>             
#>               1  2
#>   IV only    90  0
#>   IV then SC 30 60
#>   SC only     0 90
#> 
#> $synthetic
#>             
#>               1  2
#>   IV only    90  0
#>   IV then SC 30 60
#>   SC only     0 90

The arm dosed both ways keeps both, because the route is a property of each planned cycle in the arm’s dosing model rather than of the arm.

onc_sim: a dose that changes for reasons the patient’s data explains

Two hundred patients, a year of daily everolimus, tumour size beside trough concentrations. Most placebo patients cross over on progression and about a quarter of the treated arm reduce their dose, so the dose history is neither constant nor protocol-driven — it is a response to what happened to that patient.

onc_roles <- pmx_roles(
  id = "ID", time = "TIME", nominal_time = "NTIME", dv = "DV", amt = "AMT",
  evid = "EVID", cmt = "CMT", dvid = "NAME", addl = "ADDL", ii = "II",
  cens = "CENS", strata = "ARM", covariates = c("BSLD", "AGE", "SEX"),
  keep = "CROSSOVER"
)
onc_run <- pca_run("onc_sim", onc_sim, onc_roles, seed = 808)

synpmx_scorecard_datatable(onc_run$card, report = "minimal")

Nothing fails, and one row is worth reading against the source rather than against a threshold. A5b, doses per patient, reads 327 -> 266. That is counted on the doses, not on the records they are written in: both tables are expanded before the card is scored, so ADDL and II cannot flatter or distort it. The records themselves tell a different story:

records <- function(data) {
  d <- data[data$EVID == 1L, ]
  c(records = nrow(d), doses = sum(d$ADDL + 1L))
}
onc_records <- rbind(source = records(onc_sim),
                     synthetic = records(onc_run$synthetic))
onc_records
#>           records doses
#> source        322 71180
#> synthetic    3283 60488

The doses are close and the encoding is not: 3,283 records where the source has 322. The dosing model draws a reduction, a skip or a discontinuation, at rates pooled over both arms, independently at each planned cycle, and this study’s cycle is one day, so a synthetic patient’s dose flickers between levels where a real patient held one dose for months and changed it once. The regimen is close to right on average and wrong in its runs, and because pmx_compress_doses() writes each constant run as its own record, the record count is what makes it visible — which is exactly why A5b counts doses instead, and why the two numbers have to be read together.

That is a fair description of what the model is: per-cycle hazards, with no memory of how long a patient has been at a dose. On a study dosed every three weeks it is a good description of dose reduction. On one dosed daily it is not, and this is the dataset that says so.

What the nine runs held

The trial summary is the whole of what leaves the source, so its size across seven studies is worth one table.

inventory <- do.call(rbind, lapply(runs, function(entry) {
  report <- as.data.frame(pca_report(entry$summary))
  basis <- entry$summary$basis
  features <- as.data.frame(pca_features(entry$summary))
  sparsest_cell <- min(features$patients, na.rm = TRUE)
  data.frame(
    Dataset = entry$label,
    Patients = basis$n_source,
    Arms = length(entry$summary$arms$arms),
    Features = length(basis$columns),
    Components = basis$k,
    `Numbers held` = sum(report$numbers),
    `Fewest patients` = min(report$min_patients, na.rm = TRUE),
    `Set by` = if (sparsest_cell <= basis$min_group) "sparsest grid cell" else
      "smallest arm",
    check.names = FALSE, stringsAsFactors = FALSE
  )
}))
knitr::kable(
  inventory, row.names = FALSE,
  caption = paste(
    "What each run read out of its source. `Numbers held` is every quantity in",
    "the trial summary; `Fewest patients` is the smallest number of patients",
    "standing behind any one of them, and `Set by` names which quantity that is."
  )
)

What each run read out of its source. Numbers held is every quantity in the trial summary; Fewest patients is the smallest number of patients standing behind any one of them, and Set by names which quantity that is.

Dataset	Patients	Arms	Features	Components	Numbers held	Fewest patients	Set by
case1_pkpd	180	6	34	11	1925	30	smallest arm
mad	60	6	67	12	1715	10	smallest arm
warfarin	32	1	25	6	292	5	sparsest grid cell
wbcSim	45	1	21	4	195	6	sparsest grid cell
theo_md	12	1	27	2	185	12	sparsest grid cell
nimoData	12	1	41	2	272	9	sparsest grid cell
mavoglurant	120	1	18	7	254	19	sparsest grid cell
mixroute_sim	90	3	18	5	269	30	smallest arm
onc_sim	200	2	31	9	2222	51	sparsest grid cell

Two things to read across the rows. Components track the cohort, not the grid. The retained count is capped at a fifth of the subjects, so nimoData gets two components out of 41 features while case1_pkpd gets eleven out of 34: the wider grid is the one with fewer components. Twelve patients do not support more than two directions however many visits they were seen at, and that cap is what stops a basis from interpolating the subjects it was built from.

And which quantity is most exposed depends on the study’s shape. On the two six-arm studies it is the smallest treatment arm — 30 patients on case1_pkpd and 10 on mad — because the score means and covariances are fitted per arm. On the single-arm studies there is no small arm, so the most exposed quantity is a grid cell: warfarin’s sparsest modelled cell is held by 5 patients and wbcSim’s by 6. min_column_patients is the floor under that number, and it is the setting to raise if the exposure matters more than the late timepoints.

What the scorecards said

verdicts <- do.call(rbind, lapply(runs, function(entry) {
  card <- as.data.frame(entry$card)
  counts <- table(factor(card$verdict,
                         levels = c("pass", "review", "FAIL",
                                    "not applicable")))
  data.frame(Dataset = entry$label, as.list(counts),
             Failing = paste(card$check[card$verdict == "FAIL"],
                             collapse = ", "),
             check.names = FALSE, stringsAsFactors = FALSE)
}))
knitr::kable(verdicts, row.names = FALSE,
             caption = "Scorecard verdicts across the nine runs.")

Scorecard verdicts across the nine runs.

Dataset	pass	review	FAIL	not applicable	Failing
case1_pkpd	11	2	0	8
mad	12	1	0	8
warfarin	12	1	0	8
wbcSim	10	3	0	8
theo_md	12	1	0	8
nimoData	11	2	0	8
mavoglurant	10	3	0	8
mixroute_sim	12	1	0	8
onc_sim	11	2	0	8

not applicable is the same eight rows on every card, and it is a gap rather than a result. B1a, B1b and C2 read a run record this generator does not write; all three questions are answerable from the two tables, and until they are computed that way the card is silent exactly where this generator’s largest loss is — C2, the variety of dose schedules. B2 and B4a ask whether a real patient’s profile or visit set was reused, which this generator does not do at all, and E1 to E3 score a population model it does not fit.

D1 is review on all seven, and it is the row that describes the method rather than any study. It reports the furthest standard-deviation ratio over every numeric variable — endpoints and covariates alike — and a basis of a few components reproduces less spread than the data had. case1_pkpd is furthest on its PK endpoint at 0.96 of the source spread; nimoData, on twelve subjects and two components, is furthest on the BSA covariate at 0.23. On four of the seven the variable flagged is a covariate rather than an endpoint — wt, WT, BSA and HT — which is worth knowing before reading the row as a statement about the endpoints. A continuous covariate is one number per subject and carries no trajectory for the components to reconstruct it from, so it is the feature the basis compresses hardest.

A5a is the row that carries the cost of each grid, and it is review on exactly the three studies where something was merged or dropped: case1_pkpd (two pre-dose records without a value), wbcSim (the column guard) and mavoglurant (the snapping). The four studies whose grids cost nothing read pass on it.

What B4b found

B4b asks whether any generated subject’s vector of measured values equals a real patient’s that no other patient shares: the copy question, asked of what somebody’s assay said.

copies <- do.call(rbind, lapply(runs, function(entry) {
  card <- as.data.frame(entry$card)
  observations <- entry$source[[entry$roles$evid]] == 0
  data.frame(
    Dataset = entry$label,
    `Median observations` = as.numeric(stats::median(table(
      entry$source[[entry$roles$id]][observations]
    ))),
    `DV vectors copied (B4b)` = card$result[card$check == "B4b"],
    check.names = FALSE, stringsAsFactors = FALSE
  )
}))
knitr::kable(copies, row.names = FALSE,
             caption = "B4b across the nine runs, beside how long a source patient's record is.")

B4b across the nine runs, beside how long a source patient’s record is.

Dataset	Median observations	DV vectors copied (B4b)
case1_pkpd	35	0
mad	68	0
warfarin	13	0
wbcSim	4	0
theo_md	22	0
nimoData	27	0
mavoglurant	24	0
mixroute_sim	17	0
onc_sim	16	0

B4a, the same question asked of the list of observation times, does not apply to this generator. Attendance is drawn independently per visit at one rate per endpoint, onto the declared grid, so a generated list that equals a real patient’s is a coincidence with a computable chance rather than a copy. On a study whose patients’ visit lists are nearly all unique, which short records make likely, that chance is high, and a row firing on it would measure the width of the grid rather than a disclosure. The scorecard gives the measurements behind that ruling, under B4a.

What remains is that the visit draw has no floor. Nothing stops a generated set from landing on one that a single real patient held. synpmx_avatar() has such a floor: min_pattern_share draws each avatar’s attendance from patterns at least that many subjects share. On a study with genuinely distinctive attendance, a coincidence with a real patient’s visits is possible here and no row reports it.

What is preserved, and what is not

Preserved, on the evidence above: the endpoint set (A3 on mad), the level sets of discrete endpoints less any level fewer than three patients recorded, the arms and their sizes, the assay limit and the censored fraction per arm, the covariate–profile relationship that the shared feature matrix carries, and the planned dose schedule in full — including nimoData’s ten infusions, which the AVATAR run on the same study truncated.

Not preserved, and each is a decision rather than a defect:

• Visit-to-visit measurement noise. The model holds a handful of components on a visit grid, so the generated profiles are smoother than the source ones. Visible in every figure above and most obvious on theo_md, where twelve subjects buy two components.
• The variety of dose schedules. One schedule per arm, against one per patient on a study recording actuals. C2 would report it and reads not applicable; pca_dosing() reports the same quantity per arm.
• When patients missed visits, and differences between arms in attendance and dose changes. Attendance is one rate per endpoint and the dose changes three rates pooled over the arms, so a late visit few patients reached is drawn as often as an early visit of the same endpoint, and an arm whose patients dropped out or reduced more often than the others’ takes the study’s rates. onc_sim is the worked case.
• Anything the grid construction merged. mavoglurant is the worked case: 117 recorded times onto 15 nominal ones costs nine observations per patient.
• Late follow-up held by few patients, which the column guard drops. wbcSim is the worked case.
• Any study without a protocol grid. pheno_sd is refused, and that is the boundary of this generator rather than a gap in it.

Where to go next

• vignette("pca-demo") — one study end to end.
• vignette("pca-fingerprint") — every quantity the fingerprint holds.
• vignette("pca-algorithm") — how each is estimated, and what it costs.
• Evaluating AVATAR on public data — the same ten datasets under synpmx_avatar().