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, a planned dose schedule per arm, and dose-change rates and one attendance rate per endpoint pooled over the arms. 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: a planned dose schedule per arm, rates for reduction, interruption and discontinuation pooled over the arms, the visits each arm has, and one attendance rate per endpoint.
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 dropped. Its mean score vector would be those patients and its covariance would be noise around them, so synpmx_pca() warns, names the short arms and leaves them out rather than generating from one or two people. The synthetic data then has no such arm. min_arm_patients sets the floor: pool the arm, drop the column from strata, or lower the floor to keep them.

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.

A regimen recorded as one record plus ADDL and II is written out by pmx_expand_doses() before the features are built, so the dose cells count every dose the patient received; the generated study is folded back by pmx_compress_doses() on the way out and leaves in the encoding its source used (REV-051).

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 that at least min_category_patients patients hold, 3 by default. A rarer level gets no column, so it is not in the basis and cannot be generated, and the patients holding it are missing on that covariate. A covariate left with one level, which no column can carry because it does not vary, is written as that level for every generated subject.

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 its planned schedule, which is the protocol; the dose-change rates and the visit model are pooled over the arms. All of it is a summary of the study rather than a fact about a patient.

The dosing model

Each arm carries a planned schedule, and three rates, pooled over the arms, say how patients departed from their schedules. 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 every arm, so an arm whose patients reduced or stopped more often than the others’ takes the study’s rates:

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)

A rate rests on at least min_arm_patients patients. A rate whose events fewer than min_arm_patients patients had, across every arm, is their dosing history, and is zero.

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 and to one significant figure, that several patients of some arm dropped to, one ladder for the study: 0.74, 0.75 and 0.78 are one level, 0.8.

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 for the study, a planned schedule 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. On public studies, a generated subject’s list of observation times matches one a single real patient held at about the rate chance predicts: those studies’ visit lists are nearly all unique, and attendance drawn per visit lands on one of them by coincidence. Because a match there is a coincidence rather than a copy, scorecard row B4a reads not applicable for this generator, and the scorecard vignette gives the measurements. B4b, the value side, carries the copy claim, and it is zero on every public study in Evaluating the PCA generator on public data.

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 one attendance rate per endpoint, the share of the endpoint’s scheduled visits at which a patient was observed, at every visit of it an arm has; an arm has a visit where at least min_arm_patients of its patients were observed there. Where fewer than min_arm_patients patients missed any of an endpoint’s visits its rate is 1: “one patient missed a visit” discloses as much as “one patient came”, and the threshold rule of statistical disclosure control asks for enough patients on both sides. A rate per endpoint loses when in the study visits were missed, and keeps each endpoint’s number of observations per patient on average.

Attendance is drawn independently per visit at its endpoint’s rate, 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. B4b compares the finished tables directly and is the row that carries the copy claim; B4a reads not applicable, because attendance drawn per visit can match a real patient’s list only by coincidence.

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; pooled over the arms, the dose ladder and three rates
Visit model	The visits each arm has, and one attendance rate per endpoint
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, a kept one that fewer than `min_category_patients` of the arm’s patients hold written as missing. No factor level or discrete endpoint level fewer than that many patients hold is stored
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 ones the source held. Both keep only the levels at least min_category_patients patients hold, because a level held by one patient is a level that identifies them, and the B5 row of synpmx_scorecard() checks the output for one.

Subject identifiers are minted rather than derived. No source label is kept, including the levels of an identifier stored as a factor. The only thing kept about a numeric identifier is the largest one rounded up to a power of ten, so that a synthetic ID cannot collide with a real one and the largest real ID is not stored.

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, rounded down to 1, 2 or 5 times a power of ten so that the summary does not carry one patient’s measurement, 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. Thirteen of its twenty-one rows read the source and synthetic tables, compute the same way for any generator, and mean the same thing.

Eight read not applicable. B1a, B1b and C2 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. B2 and B4a ask whether a real patient’s profile or visit set was reused, and this generator draws both from a model of the cohort. E1, E2 and E3 score a fitted population model, and this generator fits none.

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 one attendance rate per endpoint and one shared schedule per arm. B4b asks the copy question directly, on the finished tables, and is the row that carries 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, or, where no amount is shared by min_arm_patients patients, the arm’s mean to two significant figures, so that variation is lost even though reductions relative to it are kept. The mode of amounts each held by one patient would be one patient’s dose.
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, a per-arm quantity by one arm, and a covariate level by at least min_category_patients patients.

Andrew Stein