The synpmx PCA Algorithm

This article is the full specification of synpmx_pca(), which generates a synthetic pharmacometric (PMX) dataset from a fitted model rather than from real patients’ values. Use of the algorithm is demonstrated in vignette("pca-demo").

synpmx_pca() makes no formal privacy claim. It does not estimate a pharmacokinetic (PK), pharmacodynamic (PD) or nonlinear mixed-effects model, so it reports no clearance, no volume and no between-subject variability on either. Use the synthetic datasets to develop and debug analysis code. Do not use them to estimate parameters.

What the principal components are used for

The components are a basis. Each subject’s trajectories become a short vector of scores, the scores become the response of a fitted linear model, and new subjects come from drawing that model’s residual.

That is what makes the output fully synthetic. Nothing selects a real patient, nothing weights one, and no number a patient measured reaches the generated table. What leaves the source data is a mean, a scale, a set of loadings, one mean score vector per arm, a residual covariance, and a dosing and visit model per arm. pca_report() inventories all of it.

Overview of Algorithm

Declare the columns. The roles variable says which column is the id, time, dose, measurement, arm and covariate. Only named columns survive.
Build one feature vector per subject. Place every row on the declared nominal grid, align each subject to their own first dose, and read each endpoint’s trajectory onto that grid, beside the baseline covariates.
Fill, standardize and decompose. Replace censored values with a draw inside the censoring region, median-fill the cells a subject has no observation in, standardize every column, and take the principal components.
Model the scores against the arm. Each arm gets its own mean score vector and its own residual covariance.
Draw new subjects by adding a fresh residual to an arm’s mean, and invert the standardization to get a feature vector back.
Fit a dosing model and a visit model per arm: a planned dose schedule with rates for reduction, interruption and discontinuation, and the probability of a visit at each nominal time.
Place the drawn values on that skeleton, snapping a discrete endpoint back onto its levels and putting the assay limit back where the source had one.
Construct the full synthetic dataset with new subject IDs and every column back in the shape the source had, and validate it.

Steps 1 to 4 and Step 6 read the source. Steps 5, 7 and 8 read only the summaries those produced. The section below says why that split is enforced rather than intended.

Step 1: Declare the Meaning of the Columns

The package never guesses PMX roles from column names. Declare them with pmx_roles(), as in vignette("pca-demo").

strata is load-bearing. It names the arms, and the arms are the groups Step 4 fits against, so a column added to strata adds a group whose mean and covariance are estimated separately.

An arm of fewer than three patients is refused. Its mean score vector would be those patients and its covariance would be noise around them. synpmx_pca() stops and names the short arms rather than generating from one or two people. min_arm_patients sets the floor. Pool the arm, drop the column from strata, or exclude those patients before calling.

Step 2: One Feature Vector per Subject

`nominal_time` is required

The grid is the model’s axis. Every feature is a cell on it, and dose rows and observation rows are placed on it together, which is what keeps the two in step. A grid inferred from elapsed time puts a sample on the wrong side of a dose as soon as the doses were recorded as actuals rather than as planned times, and choosing where a trough sample belongs is a statement about the protocol rather than something an algorithm can derive.

So synpmx_pca() requires nominal_time and refuses without it. A study that does not carry one needs the protocol’s planned times added as a column first. Rows where it is missing are named in the error rather than filled.

Each subject is then aligned to their own first dose, so that time zero means the same thing for everybody.

For each endpoint, the grid is every nominal time in the study. There is no cap on its width: the number of components is bounded by the subject count however many columns there are, so extra columns cost nothing but arithmetic.

Endpoint values are log-transformed where the endpoint is positive. Covariates join the same matrix, continuous ones as themselves and categorical ones as one indicator column per level.

Covariates and trajectories share one matrix on purpose. A single draw in Step 5 therefore carries the relationship between a covariate and the profile, which a generator drawing covariates separately and merging them on afterwards cannot do.

The minimum patients per cell

A grid cell only a few patients ever reached is filled from the median for everybody else. Such a cell describes those few patients and nobody else, so it is dropped rather than modelled. The floor is min_column_patients, defaulting to the larger of 3 and a tenth of the cohort.

This is the only place the algorithm defends the late timepoints. Follow-up thins out towards the end of a study, and the last visits are where a cell held by one patient is most likely.

Step 3: Fill, Standardize and Decompose

Missing cells take the column median, every column is centered and divided by its standard deviation, and prcomp() gives the loadings and the scores.

Median filling shrinks a subject towards the middle wherever they have no observation. That is accepted rather than corrected: the cells where it happens most are the ones Step 2 has already dropped, and an iterative fill would add machinery for an effect the column guard has bounded.

How many components

The retained count is the smallest reaching pca_variance of the total, capped at a fifth of the cohort. The cap is the part that matters. With many components against few subjects the basis interpolates the source, and a drawn subject can land on a real one; the variance rule alone does not prevent that, because a small study reaches any variance target quickly.

Step 4: Model the Scores Against the Arm

Scores become the response of a linear model with the arm as the predictor. Every arm gets its own mean score vector, and its own residual covariance shrunk towards the pooled one in proportion to its size.

dose_term = "log" is the alternative: one regression of the scores on log1p() of each subject’s total dose, which spends a single coefficient rather than a mean per arm and generates at doses the study never ran.

"factor" is the default, and the reason is measurable. "log" assumes the dose-response is log-linear and that every arm has the same between-subject spread. A lower limit of quantification (LLOQ) breaks both: the low arms sit on the floor, so the relationship is flat and then rising rather than straight, and the arms on the floor have far less spread than the top arm, so one pooled covariance smears them upward and pulls the top arm in. On a six-arm study whose source PK medians ran from the LLOQ to 0.53, "log" returned a range of 0.10 to 0.18 and "factor" recovered the arms. Use "log" where the arms are genuinely log-linear and there are too many dose levels to spend a mean on each.

Step 5: Draw and Invert

A new subject’s scores are their arm’s mean plus a fresh multivariate normal residual drawn through the Cholesky factor of that arm’s covariance. Multiplying by the loadings, then by the column scales, then adding the column centers, gives a feature vector in the original units.

The residual draw is what makes this generation rather than blending. Nothing selects a real subject, and nothing weights one.

Step 6: The Dosing Model and the Visit Model

Each generated subject is assigned an arm, keeping each arm’s share of the cohort. The arm carries two models, and both are summaries of the arm rather than facts about a patient.

The dosing model

Each arm carries a planned schedule and three rates that say how its patients departed from it. pca_dosing() reports the first and pca_dose_rates() the second.

The planned schedule is one row per cycle: the nominal dose time, and the amount the arm was planned to receive there. The amount is read per cycle from the patients still on their starting dose, so a protocol-prescribed escalation is part of the plan, while accumulated reductions do not drag the plan down with them.

The three rates are discrete-time hazards, pooled over the arm:

Rate	What it is
`discontinuation`	P(this cycle is the last, given still on treatment)
`interruption`	P(this cycle is skipped, given still on treatment)
`reduction`	P(drop one level, given still on and a lower level exists)

The dose ladder is built from within-patient decreases. A study dosed by body weight gives every patient a different amount, so the pooled spread of amounts scatters continuously and any threshold on it would invent a ladder out of ordinary between-patient variation. A reduction is something else: the same patient receiving less than they did at the previous cycle. The ladder is the set of fractions, of a patient’s own starting dose, that several patients dropped to.

Why there is no detector

A study where nobody reduces, skips or stops early has all three rates at zero and a single level, and the model then reproduces the planned schedule exactly for every patient. Nothing decides which kind of study this is, because the rates already say it. On a fixed-schedule study the generated dose rows, the doses per patient and the scorecard’s A5b row are identical to what the modal schedule produced.

Where the rates are not zero — oncology being the usual case — generated patients differ from one another in the same way and to the same degree the source patients did. Relative dose intensity, time on treatment and the fraction of patients ever reduced are often what such a dataset exists to analyse, and a generator handing everyone the modal schedule has removed all three.

What this does and does not protect

No patient’s own sequence of reductions is copied. What leaves the source is three probabilities and a ladder per arm, and the planned grid stops at the last cycle min_arm_patients patients reached, so a treatment duration only one patient had cannot be generated — the same reasoning as SIM-047 on the observation side.

What it does not give is a guarantee of non-coincidence. Simulating a schedule per patient means a generated sequence can land on a real patient’s by chance, exactly as the visit draw can. Scorecard row B1b asks that question and cannot answer it here, for the reason given under “Reading the Scorecard on This Output”. Read distinct in pca_dose_rates() against the arm size instead: a study whose schedules are nearly all unique is one where a coincidence is least likely to be meaningful, and a study with few distinct schedules is one where it is not a disclosure at all.

That coincidence has been measured on the visit side, and it is not hypothetical. Over seven public studies, B4a — which asks the same question of the finished tables — finds a generated subject’s list of observation times matching a singleton real one on two of them, reproducibly across seeds. The matches are between short lists on sparsely sampled studies rather than between distinctive schedules, and B4b, the value side, is zero everywhere. The numbers and the reasoning are in Evaluating the PCA generator on public data, along with what a floor on the visit draw would cost.

The visit model

pca_visits() reports it: one row per arm, endpoint and modelled nominal time, giving the probability that a generated subject in that arm has an observation there. It is the fraction of the arm’s patients who did.

Attendance is drawn independently per visit from that probability, so no real patient’s set of attended visits is reused.

Why the timing survives

Dose rows and observation rows are placed on the same declared nominal grid, which is what Step 2 requires nominal_time for. The consequence is that the event structure making a dataset fittable comes through at its source values: doses per patient, the time of the last dose, the end of follow-up, the baseline sample before the first dose, and the distribution of time after dose.

What does not come through is the variety of schedules, as above.

B1a, B1b and C2 cannot be scored on this output, for the reason set out under “Reading the Scorecard on This Output” below. B4a and B4b compare the finished tables directly and are the rows that carry the copy claim.

The Model Is the Only Thing Generation Reads

Fitting produces one object, and generation is a function of that object and a subject count. No patient row is in scope while a synthetic dataset is being written, and that is enforced rather than intended: the generator’s only arguments are the model and the count, and a regression test regenerates a dataset from the model alone and requires it to match.

The object holds, and this is all of it:

Group	What it is
Basis	Column means, column standard deviations, and the component loadings
Scores	One mean score vector per arm, and one residual covariance per arm
Dosing model	Per arm, the planned cycles, the dose ladder, and three rates
Visit model	Per arm, endpoint and time, the probability of an observation
Schema	Column order, column types as empty prototypes, the discrete endpoints’ level sets, the log-or-identity choice and assay limit per endpoint, and one value per arm for the strata and kept columns
Arms	The number of patients in each

pca_report() prints the inventory with the count of numbers in each group and, for each, the smallest number of patients standing behind any one of them. vignette("pca-fingerprint") walks through every group in turn, on a worked study.

Two entries deserve naming because they are not averages. The level set of a discrete endpoint is read from real values: a 0/1 endpoint has no third value to emit, so a generated value is a source value. A categorical covariate’s levels are likewise the ones the source held, which is what the B5 row of synpmx_scorecard() measures — a level held by one patient is a level that identifies them.

Subject identifiers are minted rather than derived. The only thing the model keeps about the source’s identifiers is the largest one, so that a synthetic ID cannot collide with a real one.

Step 7: Place the Values

Each attended visit takes its drawn value for that endpoint, inverted out of the log transform where one was applied, and snapped back onto the source’s levels where the endpoint is discrete, so that a 0/1 endpoint cannot come back as a number between its levels.

Categorical covariates come out of the draw as one number per level, and the level with the largest is the one written.

The assay limit

The drawn value is the latent one: what the subject would have measured with no assay limit. Where the source declared censoring, the reported value, CENS and LIMIT are then reconstructed from it together, so a value below the lower limit of quantification (LLOQ) is reported on the boundary and flagged, as the source reports it.

This runs in both directions and both are needed.

Censored source values are replaced by a draw inside the censoring region before the basis is fitted, back in Step 3. The log transform’s offset is read from the reported values rather than from those draws, and that distinction is load-bearing rather than a detail. The offset is half the smallest positive value and exists to keep log() finite near zero; .impute_censored() draws uniformly on the raw scale, so a value censored at a limit of 0.05 can land three orders of magnitude below it and take the offset with it. Every censored observation then occupies several log units of what is purely uniform noise, the per-arm residual covariance is fitted on that inflated spread, and exponentiating a wide Gaussian has a heavy upper tail — so an arm that is entirely below quantification comes back with detectable concentrations that rise with time.

Not imputing at all is worse, and was measured. With every censored value sitting exactly at the limit, the arm’s mean score is the boundary itself and a symmetric residual puts half the draws above it: on a six-arm study the lowest arm returned 42% below the limit against a source 88%. The draw is what moves the arm’s centre to where the latent values are. A grid cell where most subjects sit exactly at the limit has almost no variance, so a basis fitted on the reported values would learn the assay rather than the patients, and the generated data would inherit a floor slightly above the real one. This is the only random step in summarizing, which is why synpmx_pca_summarize() takes a seed of its own.

Reapplying the boundary at emit is what makes a low-dose arm come back flat. Without it every value below the limit is emitted as itself, and an arm that the study recorded as a straight line at the LLOQ comes back as a spread of small numbers — visibly wrong on a log-scale concentration plot, and wrong in a way that would carry into any BLQ handling the dataset is used to develop.

Step 8: Construct and Validate

New subject IDs, the strata and keep columns carried from the assigned arm, every column back in the source’s order and type, and validate_pmx() over the result.

Reading the Scorecard on This Output

synpmx_scorecard() was built to judge a synthetic dataset rather than a generator, and most of it holds here. Fifteen of its eighteen rows read the source and synthetic tables, compute the same way for any generator, and mean the same thing.

Three do not: B1a, B1b and C2. They read a run record that synpmx_scorecard() expects to find on the generated table as a pmx_settings attribute, rather than reading the two tables, and this table carries no such record. All three report not applicable.

For B1a and B1b the guarantee is structural rather than unrecorded. No individual’s visit set and no individual’s dose schedule exists in the model to be copied, because the model holds a per-visit probability and one shared schedule per arm. B4a and B4b ask the copy question directly, on the finished tables, and are the rows that carry the claim.

C2 is the one to be careful about, because a blank row is at its most misleading exactly where the loss is largest. C2 asks how many of the study’s distinct dose-time schedules survived, and this algorithm generates one schedule per arm. On a study whose recorded dose times are actuals, that is a large loss and the card says nothing about it. pca_dosing() reports the same quantity as distinct and share, per arm, and is what to read in its place.

What This Algorithm Does Not Preserve

Within-visit detail. The model holds a handful of components on a visit grid. Measurement noise from one visit to the next is not in them, so the generated profiles are smoother than the source.
Anything the components discarded. Reconstructing the source subjects from their own retained scores shows how much that is, per endpoint. It is the diagnostic that says whether the retained count is large enough.
Interpretable pharmacology. A score standard deviation is a variance decomposition, not a parameter. There is no clearance here to report a coefficient of variation on. A loading plotted against time can be described — flat and positive is overall magnitude, crossing zero separates early from late — but naming a mechanism for it is a claim the fit does not support.
Doses and times the study never ran, under the default dose_term. "log" extrapolates; "factor" does not.
Amounts that vary continuously between patients. Body-weight or body-surface-area dosing gives every patient their own milligrams; the planned amount here is the arm’s modal one, so that variation is lost even though reductions relative to it are kept.
Within-subject dose escalation. A subject’s total dose is one number and an arm’s schedule is one schedule, so a study escalating within a patient is represented by its arms rather than by the escalation.

What Leaves the Source Data

The inventory is above, under “The Model Is the Only Thing Generation Reads”, and pca_report() prints it for a given run. The column to read is the smallest number of patients standing behind a quantity. A grid cell mean is backed by most of the cohort and a per-arm quantity by one arm, while a rare covariate level can be backed by a single patient.

Andrew Stein