The synpmx PMX Model Algorithm

Andrew Stein

This article is the full specification of synpmx_model(), which fits a small set of linear pharmacokinetic (PK) models to a study and generates new subjects by simulating from the fit. Use of the algorithm is demonstrated in vignette("pmxmodel-demo").

The aim is a simple model that roughly reproduces the properties of the study, not one that characterises it. The synthetic profiles it generates are for prototyping analysis code. The goodness of fit of the selected model is not assessed: model_report() says what the search had to go on, not whether the model describes the source. The candidate set is the package’s linear models and the covariate model is allometric scaling or nothing, which is too little to answer a scientific question, so the fitted parameters are not estimates to report, even though the printed object looks exactly like the output of a real population analysis.

No formal privacy guarantee is offered. The protections on what leaves the study, and the checks that measure them, are set out in Privacy protections in the PMX model generator; every number that leaves is listed in vignette("pmxmodel-fingerprint").

Three Exported Functions

Function	Reads patient data	Needs nlmixr2	Returns
synpmx_model_estimate(data, roles, ...)	yes	yes	pmx_fitted_model
synpmx_model_generate(fitted_model, n_subjects, seed)	no	no	data frame
synpmx_model(data, roles, n_subjects, seed, ...)	yes	yes	data frame

model_report() inventories everything the fitted object carries, model_candidates() returns the attempted models and their acceptance results, and model_parameters() returns the fixed effects, the between-subject covariance matrix and the residual error. model_fingerprint() returns the part of the fit that generation reads, which is what synpmx_model_generate() attaches to its output, and model_privacy_checks() runs five checks on it (Step 5).

nlmixr2est is an imported dependency, used by synpmx_model_estimate(), and the candidate set is exactly the models the generator can evaluate in closed form. A candidate the fitter could estimate and the generator could not simulate would be a model that fits and then generates nothing, so the two lists are one list. The vignettes read a stored synpmx_model_estimate() fit, so R CMD check and the pkgdown site never compile a model.

The Arguments

Beyond the roles and the seed:

Argument	Default	Effect
pk	NULL	One of the built-in models, forcing it and skipping the search; or several, which is how a search of your own choosing is asked for.
pd	NULL	Named vector of PD shapes per endpoint. Skips that search.
pd_by_arm	FALSE	Fit each PD endpoint’s shape, and read each binary or ordinal endpoint’s level frequencies, per arm rather than once over the pooled cohort.
endpoint_roles	NULL	Declares concentration endpoints with pk, or a PD-only study with c(pd = "response"). More than one concentration may be named.
start_param	NULL	Starting values for the population fit, keyed by endpoint where more than one concentration is fitted. The escape hatch where the non-compartmental read of the median profile starts the optimizer somewhere it cannot move from.
covariate_effects	"none"	"none" puts no covariate in the structural model; "auto" applies allometric scaling on clearance and volume where a weight-like covariate is declared.
min_subjects	20L	Warn and fit anyway below this cohort size.
min_arm_patients	3L	Warn and drop the patients in any arm below this, as synpmx_pca_summarize() does.
min_time_bins	6L	Warn and fit anyway below this many distinct nominal times after a dose; no post-dose observation at all still refuses.
max_fit_subjects	60L	Fit the population model to this many subjects, drawn in proportion to the arms. The dosing, visit and covariate models read every subject.
estimation	"focei"	Passed to nlmixr2.
drop_influential	TRUE	Leave out of the PK, PD and covariate estimates any patient who moves a released estimate too far, and estimate again (Step 5). FALSE keeps every patient and only reports the reading.

The default tries two compartments first, with a one-compartment fallback only if the first fit fails acceptance checks. An ambiguous administration route gets this sequence for each possible route. A named pk model is fitted alone; a vector of names requests an Akaike information criterion (AIC) comparison of acceptable fits.

Two Time Axes, and nominal_time Is Required

Estimation reads recorded times. A population PK fit is a statement about time after the dose that was actually given, so the subject-level dosing records go to the fitter as recorded, reductions and skipped cycles included. Fitting against the planned schedule where a patient’s dose was reduced would push the drop in concentration that followed into clearance, and the model would report a population that eliminates the drug faster than the real one.

Everything else reads the nominal grid. The dosing model, the visit model, and every signal in Steps 1 and 2 sit on the protocol’s planned times. On recorded times a “median profile” is one point per time, because no two patients share a clock reading, and a study sampled at six times looks like a study sampled at two hundred. A shape, a peak position and a count of sampling times are all statements about the protocol.

Inferring that grid is a statement about the protocol only the caller can make, so nominal_time is required rather than derived.

Repeated doses are written out before anything reads them. A regimen recorded as one record plus ADDL and II — “and 364 more like it, every day” — is expanded by pmx_expand_doses() on the way in, so that every dose the patient received is a row to the dosing model, the derived time after dose and the estimation table alike. Before this, every reader counted rows: onc_sim showed the fitter 237 dose records standing for 71,180 doses and got back a clearance six orders of magnitude off, from a fit that had moved and so tripped no gate (REV-051). The fitter is handed the compact form again, because nlmixr2 reads it natively and is fastest on it, and the generated study is folded back by pmx_compress_doses() on the way out, so it comes back in the encoding its source used.

Overview of Algorithm

1. Classify the endpoints. Four signals decide which endpoint or endpoints are drug concentrations. Everything else continuous becomes a pharmacodynamic (PD) endpoint.
2. Detect the design. The route of administration determines the candidate models. Sampling richness is reported alongside the fit.
3. Estimate the parameters. Fit each surviving PK model as a population nonlinear mixed-effects model through nlmixr2, screening before fallback, once per concentration endpoint. No covariate enters the structural model unless covariate_effects = "auto" asks for it. Fit each PD endpoint’s time course by least squares.
4. 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.
5. Store the fit: the structural model, the parameters, the arm models and the schema, and no per-subject quantity.
6. Generate new subjects. Draw covariates and random effects, draw the dose schedule, evaluate the profile against that schedule, add residual error and put the assay limit back.

Steps 1 to 4 read the source. Step 6 reads only what Step 5 stored, which is why generation runs on base R and needs neither nlmixr2 nor a compiler.

Step 1: Identify the Drug Concentration Endpoint

pmx_endpoint_types() already infers continuous, integer, binary and ordinal from the values, and that inference is reused unchanged. What this generator adds is a second classification on top of it. Four signals, computed for every endpoint that is a time course:

1. Compartment. The cmt role puts the endpoint in the compartment the doses go to, or one above a dosing compartment nobody observes, which is the depot-and-central convention.
2. Post-dose only. Observations at or before the first recorded dose are absent, zero, or at or below the censoring limit, in most subjects. Recorded observation and dose times are compared on the same clock. Later troughs do not count as baseline.
3. Shape. The median profile within the richest dose interval rises to a maximum and declines without rising again.
4. Dose proportionality. Between the highest and lowest dose level, the ratio of median subject maxima is within a factor of two of the dose ratio.

Signals 2 and 4 are required; 1 and 3 break ties. Being absent before the first dose and scaling with the dose are properties only a drug concentration has. A compartment number is a convention the dataset’s author chose, and a rise-and-fall shape is one many biomarkers also have.

Signal 4 is not computable on a single-level study, nor on one dosed by body weight. A study that gives every patient a slightly different amount has no dose levels: reading each amount as its own level compares a median of one patient against a median of one patient over a dose ratio near 1, which almost any endpoint passes. Where the distinct amounts outnumber half the cohort the signal reports “not computable” and signals 1 to 3 decide.

Where every continuous endpoint fails the baseline signal, the study is inferred as PD-only. Otherwise, if no endpoint passes or two pass and neither tie-break separates them, the function errors and names endpoint_roles as the way through. That is the inference-versus-declaration fork, answered the way dose_covariate answers it: infer where the data settles the question, offer the declaration as an override.

A study can contain only PD endpoints

endpoint_roles = c(pd = "response") declares a response-only study. endpoint_roles = list(pk = character()) makes the same declaration for all observed endpoints. The pipeline fits continuous responses using the PD shapes in Step 3 and discrete responses using visit frequencies. It skips compartment selection, population PK estimation, and the PK acceptance checks. It also skips the post-dose coverage gate: dose records need not exist, although the nominal observation grid remains required.

A positive baseline at the first dose time identifies wbcSim as PD-only without an override. Missing doses likewise provide no evidence of PK. An immediate post-bolus sample sharing the first dose timestamp, a pretreated patient, or an endogenous concentration can contradict this baseline heuristic; use endpoint_roles = c(pk = "response") when the assay is known to measure a concentration. The available continuous PD shapes remain constant, linear and exponential in study time, with no exposure term or cell-turnover model.

A PD-only fitted object has structural = NULL, parameters = NULL, and an empty pk_models list. model_candidates() returns an empty PK candidate table; the PD candidates are in model_report(fit)$pd[[endpoint]]$candidates, and model_parameters(fit)$pd returns the PD fits. Generation still carries the declared dosing, observation, covariate and arm summaries.

More than one endpoint can be a concentration

A study measuring two drugs, or a parent and its metabolite, has two drug concentrations, and each is fitted its own population model — its own structural model, its own parameters, its own residual error. No correlation between their random effects is estimated, which is a limitation rather than a claim: a patient with a high parent clearance is given an independently drawn metabolite clearance.

Which doses drive which endpoint is a declaration. By default every dose record drives every concentration: one administration, two analytes, which is what a parent and its metabolite are. Two different drugs given together are the opposite case — each model has to see its own drug’s doses and not the other’s, or a clearance and a volume come back having absorbed an input the patient never received of that drug.

The data cannot settle it. A metabolite has no dose records of its own, so “this endpoint’s doses are the ones marked 2” and “this endpoint has no doses” are the same table, and guessing wrong leaves a model with no input at all. So it is declared, on the argument routes is declared on:

pmx_roles(
  adm            = "ADM",
  routes         = c("1" = "extravascular", "2" = "extravascular"),
  dose_endpoints = c("1" = "drug A PK", "2" = "drug B PK")
)

In a NONMEM dataset the compartment already does this job, so declare the column twice — cmt = "CMT", adm = "CMT" — and key the mapping on compartment numbers.

Declared, each drug gets its own dose schedule per arm — its own amounts, its own interval, its own ladder of reductions — and each endpoint’s fit and generated profile read only that drug’s doses. Including the dosing history each observation is read against: time after dose, the dose interval and the first dose amount follow the endpoint’s own drug, which is what the starting values, the route reading, the sampling richness and the dose-proportionality signal are all computed from. Without that half, a fit table split correctly still starts from a volume scaled by the wrong drug’s milligrams. The generated study then writes both drugs’ dose records back with their own adm value and compartment.

Measured on a simulated combination — 100 mg weekly at cl 5, v 50, ka 1.2 beside 900 mg fortnightly at cl 20, v 200, ka 0.5, sixteen subjects — the second drug comes back cl 14.4, v 3.26, ka 9.34 with the doses pooled, and cl 19.3, v 166, ka 0.57 with them separated, which is what that data fitted on its own gives.

Undeclared, model_report() says so on a study with more than one concentration, which is the case where it matters.

Inference promotes a second endpoint only on signal 4. Dose proportionality is the positive evidence that something is a drug: a concentration scales with the dose and a biomarker does not. Signal 2 is necessary and not sufficient, and it is satisfied by any endpoint with no pre-dose observation. The distinction matters because signal 4 is often not computable, and the required test treats that as passing — which is right for choosing one endpoint and wrong for promoting a second. Unknown dose proportionality does not promote a second endpoint into PK. Positive baselines independently identify warfarin’s prothrombin activity and onc_sim’s tumour size as PD.

A study whose two endpoints both pass the baseline signal and scale with dose is classified with two concentrations. If the remaining evidence cannot separate them, the caller must declare the endpoint roles.

endpoint_roles = list(pk = c("parent", "metabolite")) declares it where inference cannot. Declaring several is worth doing wherever it is true: a concentration demoted to a pharmacodynamic endpoint is fitted a time course with no dose term in it, so the generated values lose their dose ordering entirely — measured at 1.05 times across a three-fold dose range, against 2.7 in the source (SIM-080).

model_report(fit)$endpoints$signals
#>   endpoint compartment post_dose shape proportional
#> 1       cp          NA      TRUE  TRUE           NA
#> 2      pca          NA     FALSE FALSE           NA

What this cannot do. On a study with one continuous endpoint and one dose level the classification rests on signal 2 alone, and signal 2 passes for any endpoint with no pre-dose observation, including a white-cell count. The remedy is endpoint_roles, and the evidence is in the table above: a row reading “not computable” under proportional is the classification saying how little it had to go on.

Step 2: Detect the Design

Detection prunes the candidate set; the fit picks among what survives. Neither property below decides the model on its own.

Route

A declared route is not detected at all. pmx_roles(adm = "ADM", routes = c("1" = "iv", "2" = "extravascular")) names the administration column and says what its values mean, and the two are declared together or not at all: which id is which route is a convention of the dataset and cannot be read off the numbers, so a guess would put every dose in the wrong compartment without failing.

Where a study declares both routes, the route is a property of each dose record rather than of the study, and the candidates are the mixed models — the intravenous form for the doses given intravenously and the extravascular form for the rest, summed. That sum is exact rather than an approximation, because every model in the closed-form set is linear in dose, and it is what lets one patient receive both. Bioavailability scales the extravascular doses and not the intravenous ones, which is what bioavailability means and why it is identifiable in a study dosed both ways and not in one dosed a single way. It is fitted as a fixed effect with no between-subject term: a design with one dose each way identifies the contrast between the routes, not a per-subject distribution over it.

A mixed model is fitted with explicit compartments rather than the closed form, and bioavailability is the whole reason. Under linCmt() the f(depot) directive does not reach the doses the data sends to the depot: evaluated at known parameters it scaled the intravenous doses instead, so f was estimated as its own reciprocal and clearance and volume as cl/f and v/f — a fit whose own predictions tracked the data while reporting parameters that meant something else (SIM-077). Written as d/dt(depot) and d/dt(central), the compartments are real, f() binds to the one it names, and the predictions match the analytic profile on both routes. Only a mixed study pays the solver cost; a study dosed one way has no f to place and keeps the closed form.

Starting values for a mixed study are read one route at a time, for the same reason the fit has to be: pooled, the cohort’s median profile is an intravenous decline and an extravascular rise averaged together, and no non-compartmental quantity read off it describes either route. Disposition comes from the intravenous records, where f is not in the way; absorption from the extravascular ones, where there is a peak to find; and their two clearances give f itself a starting value, since an extravascular read returns cl/f.

Detection below runs only where nothing was declared. It answers with one route for the whole study, which is the right answer for a study dosed one way and no answer at all for a study dosed two ways — the rate test reads a single nonzero rate anywhere as “this is an infusion study”, so a mixed study came back as an infusion model with no absorption in it.

A rate role carrying a nonzero value is an infusion, and the infusion models are the candidates. Otherwise one property separates the two remaining routes, and it is not how many patients peak early.

A drug given by mouth cannot be in the blood at the moment it is swallowed. So a concentration observed at time zero after a dose is intravenous, and a median profile that rises before it falls is oral. A profile that declines from a first sample drawn some time after the dose is both of those at once, an intravenous bolus or an oral dose whose absorption finished before anybody looked, and no amount of reading the source separates them. There both routes are offered and AIC settles it.

Counting subjects instead is the mistake this replaces. In a study whose sampling starts well past the peak, every subject’s profile declines from its first point, and a per-subject vote reads an oral study as intravenous: it counts the sampling schedule rather than the drug. The share of subjects whose own profile rises is still reported, because it says how much of the cohort the median profile speaks for, but it decides nothing.

model_report(fit)$design$reason
#> [1] "the median profile rises to a peak at 9 before declining, and 31% of subjects do too"

Sampling Richness

Across the cohort, distinct nominal times after a dose holding an observation. Below six the fit still runs and warns: a one-compartment model is not identifiable from that sampling, so its parameters stay near their starting values and the simulated profiles are a plausible shape rather than this study’s. There is no simpler structural model to drop to — the simplest candidate is already the one-compartment model. synpmx_avatar() and synpmx_pca() carry sparse sampling without fitting a structure to it. No post-dose observation at all is an error, since there is then no curve to fit.

Per subject, distinct nominal times after a dose within the richest dose interval. Four or more, with two after the median peak, is labelled richer sampling. A sample before the peak is required for that label only where there is an ascending limb to sample, since an intravenous bolus peaks at the dose. This sampling heuristic does not establish parameter identifiability.

This count decides nothing. Both compartment counts remain available, and model_report() prints it beside the selection, because what it answers is how much within-subject time coverage the fitted curve rests on.

Choosing the Compartment Count

The default first fits two compartments for the detected route. A fit that passes the checks in Step 3 is accepted without fitting one compartment. If it fails, one compartment is tried under the same checks. If both fail, estimation stops and lists the reasons. A stable two-compartment fit may be accepted even when a one-compartment model could describe the study.

Where the route is ambiguous, each route gets its own fallback sequence. The accepted models are then compared by Akaike information criterion (AIC). pk = "1cmt_oral" or pk = "2cmt_oral" requests exactly that model, with no automatic substitution. A vector such as pk = c("1cmt_oral", "2cmt_oral") fits both and selects the acceptable model with the lowest AIC, without an extra preference for two compartments.

Step 3: Estimate the PK and PD Parameters

The population-fitting and acceptance operations below apply only to PK endpoints. A declared PD-only study proceeds directly to the PD time-course fits later in this step.

The population model is fitted to at most max_fit_subjects patients, 60 by default, drawn in proportion to the arms with the run’s seed. Fit time is linear in subjects and goes into the per-subject inner loop — on onc_sim, 40 patients fit in four minutes and 200 in twenty-five — while the parameters a synthetic study needs are settled long before the sixtieth patient. Every other model in Step 4 still reads the whole study, and the report states the count.

The candidates are the structural PK models Step 2 left standing. Each one is fitted as a population nonlinear mixed-effects model through nlmixr2, which estimates the fixed effects, the between-subject covariance matrix and the residual error, and where more than one candidate is offered they are compared by the rule above. The PD endpoints are estimated separately, and not as mixed-effects models.

The candidate set is drawn from the package’s closed-form set and never goes outside it. A candidate the fitter could estimate and the generator could not simulate would be a model that fits and then generates nothing, so the two lists are one list. That is also why no dose effect on clearance is fitted: it would mean a candidate with no closed form, and generation would need a solver.

A subject the endpoint was never measured in is not fitted. It adds no term to the likelihood, and a record carrying doses and no samples is at best ignored by the solver. Its dosing and its visits are still real and still reach the models of Step 4, which read the study rather than the table handed to the fitter. How many subjects that left out is said before the fit, and it warns where it takes the fitted cohort under min_subjects.

Each candidate carries log-normal between-subject variability on every structural parameter and a proportional residual error. An endpoint holding values at or below zero after censoring is handled gets an additive error instead, and the substitution is recorded on the fit.

The concentration is fitted with its censored observations censored. Where the study declares a censoring column, the rows below the assay limit reach nlmixr2 as censored – the limit in DV, the flag in CENS, and the other bound in LIMIT where the study reports one – so each contributes the probability of falling below the limit rather than a value nobody measured. The residual error determines the scatter of generated observations, so substituting invented concentrations in the likelihood would change the quantity that generation uses. The case1_pkpd example retains an optimizer false-convergence warning while passing the parameter and generation checks.

Everything else reads an imputed value, because nothing else has a likelihood to put censoring in: the visit model, the covariate model and the PD shapes are least-squares fits, and a stack of identical boundary substitutions would bend each of them toward the limit. Those rows are replaced by a uniform draw inside the censoring region – a draw rather than a fixed LLOQ/2, which would swap one artificial spike for another – and the boundary is put back when the synthetic data is emitted. The share of each endpoint sitting below the limit is reported with the fit either way.

Starting values come from a non-compartmental reading of the source rather than from a guess, because a population fit started far from the answer either converges slowly or reports the starting values back. The terminal slope is fitted by log-linear regression over the points after the peak, and everything leans on it. Clearance is the dose over the area extrapolated to infinity, not over the trapezoid alone, which on a study sampled to four half-lives understates the area by about a tenth and overstates clearance by the same. Volume is clearance over the terminal slope, which is the quantity the terminal phase identifies, since dose over the peak is not a volume of any kind for an oral dose. Absorption is solved from where the peak falls: for a one-compartment oral model the peak sits at log(ka/ke)/(ka - ke), one equation in one unknown once the slope has given ke. For a two-compartment model the steady-state volume comes from the mean residence time, with the absorption mean taken back off for an oral dose, and is split between the central and peripheral compartments.

The default tries two compartments, with one compartment as fallback. Allometric scaling, when requested, is folded into each fit.

synpmx_model_estimate(data, roles) # two compartments, then fallback if needed
synpmx_model_estimate(data, roles, pk = "2cmt_oral") # this model only
synpmx_model_estimate(data, roles, pk = c("1cmt_oral", "2cmt_oral")) # AIC

Acceptance Checks

A finite AIC alone does not establish convergence. For first-order conditional estimation with interaction (FOCEi), the optimizer must return a successful termination code and a finite AIC, except that false convergence alone is a warning. The report explains it as:

Warning: the optimizer stalled before convergence was confirmed. Check that the synthetic data reasonably reproduces the source data’s patterns and variability.

A fit with that warning can be accepted if its parameters and generated profiles pass the remaining checks; converged stays false and accepted is true. Iteration or evaluation limits and other failures still reject the candidate. Stochastic approximation expectation-maximization (SAEM) has no equivalent termination code: the check compares mean parameter values in two halves of the final fifth of its trajectory (at least ten rows). A change exceeding 50% of the larger of one and the first-half absolute mean rejects the fit. A trajectory with fewer than 20 usable rows or non-finite values in its checked tail also fails. This is a coarse drift screen, not proof that a stochastic fit has converged. A missing SAEM AIC only prevents comparison when multiple accepted candidates remain.

Fixed effects must be finite and positive, residual error finite and non-negative, and the between-subject covariance finite, symmetric and positive semidefinite within numerical tolerance. Standard errors and peripheral-parameter precision are not acceptance criteria. Estimates remaining within 1% of their starts warn; little movement alone does not prove failure.

The generation screen draws four new subjects for each fitted subject’s recorded dose and sampling design, using the generator’s population-effect, profile and residual-error functions. It uses expanded dose records, including infusion durations and mixed administration routes, and declared allometric scaling. It never uses the original subjects’ fitted random effects. An internal seed of 104729 makes the screen repeatable and preserves the caller’s random-number stream.

Missing observed designs, non-finite parameters or concentrations, nonpositive drawn parameters, and concentrations below -1e-8 fail. Tiny negative roundoff is clamped to zero. Each source reading supplies lower and upper bounds: an exact reading has identical bounds, a left-censored reading lies between zero and its limit, and an unbounded right-censored reading lies above its limit. A declared second interval boundary tightens those bounds. Inconsistent bounds fail. The medians of the lower and upper bounds enclose the source median. A generated median below 0.01 times a positive lower bound or above 100 times a finite positive upper bound fails. With uncensored data this is the ordinary median ratio check.

More than 10% of simulated profiles exceeding 100 times their source-design maximum also fails. The source maximum uses the upper observation bounds; designs with an unbounded right-censored reading are omitted from this check, and that omission is reported. These are engineering screens for gross failures, not validated goodness-of-fit cutoffs.

Less decisive discrepancies warn: a generated median more than twofold outside the source median bounds, an uncensored generated/source interquartile-range ratio outside 0.5 to 2, or a typical profile differing by more than fourfold at the median positive uncensored observation. Censoring disables the spread comparison, and the omission is reported. Numerical checks and comparisons supported by the remaining bounds still run. The pooled comparisons can miss time- or arm-specific errors; inspect the source-versus-synthetic plots for the profile’s rise, decline and spread.

model_candidates() lists attempted candidates in order, with separate converged and accepted columns and warning or rejection reasons in note. A fallback that was not needed has no row. model_report() states why the model was selected. The measured time includes fitting and acceptance checks.

model_candidates(fit)
#>       model converged accepted      aic seconds note
#> 1 2cmt_oral      TRUE     TRUE 922.2948  50.869

Covariates

The default "none" puts no covariate in the structural model. A synthetic study does not need the relationship: covariates are generated from the source’s own study-wide distributions either way, so a clearance that moves with weight buys the generator nothing.

covariate_effects = "auto" applies allometric scaling on clearance and volume where a weight-like covariate is declared, and fits nothing else. The exponents are the standard 0.75 and 1 rather than estimated ones, and the effect is asserted rather than tested: it is folded into each fit, not compared against a model without it. Testing it would double the cost of the whole call. The reference weight the scaling is centred on is the cohort’s median rounded to one significant figure, 70 rather than 71.7, because the median of an odd number of patients is one patient’s weight and the reference is stored with the fit.

The covariate is recognised by name, wt, weight, bw and the like, and must be numeric and positive. There is no way to recognise a body weight from its values alone.

The cost of this default is explicit. A covariate that influences the real profiles and is not in the model is generated independently of them, so the synthetic data carries no relationship between the two. model_report() reports the correlation between each declared covariate and the individual random effects, which is where an unmodelled relationship shows up.

Pharmacodynamics

Every remaining continuous endpoint is fitted against a constant, a linear and an exponential time course and selected on AIC, with between-subject variability on the baseline. pd = c(endpoint = "linear") names the shape and skips the search for that endpoint; a named shape that will not fit is an error rather than a silent fallback. There is no minimum number of observations: a shape is a candidate where it has one residual degree of freedom, which is what stops a line being fitted exactly through two points and then winning the comparison it was never tested by. Below any candidate the endpoint is a constant at its mean, which is what an endpoint measured once is, and the report says that rather than naming a search that did not run.

These shapes have no concentration term, so a PD endpoint driven by exposure is reproduced as a time course that happens to resemble the average subject’s response.

One shape for the whole cohort by default, and one per arm on request. The pooled fit predicts a single number for every arm, so a synthetic patient’s dose does not reach their response. pd_by_arm = TRUE fits the shape within each arm instead, which is the same per-arm summary the dosing and visit models already are and asserts no dose-response form – an arm is fitted on its own observations, on the same three-candidate ladder, or not at all. It is not the default because it is a different claim about the study: the pooled shape says these endpoints are a time course, the per-arm shape says each arm has its own. An arm holding too little of an endpoint to fit keeps the pooled shape, so no arm is left with nothing to generate.

The between-subject spread and the residual are read around each subject’s own curve, not around the population one. Taking the residual as the spread of every point about the typical curve puts all of the between-subject variation into it, and generation then emits that as independent noise per observation. On a study where subjects share a shape and sit at different levels, almost the whole of the structure comes back as scatter and no synthetic subject has a profile at all. .pd_profile() is linear in baseline for every shape, so a subject’s own baseline is a least-squares projection needing no second optimizer, and it is the same quantity generation draws: the between-subject term is the spread of those baselines on the log scale, and the residual is what is left around each subject’s own curve.

A shape that fails to fit stays in the table. The exponential is tried from four start sets rather than one, because the single median-based set fails on any response that falls and then recovers, which is the shape a turnover endpoint has. Where none of the four converges the candidate is reported with converged false and the solver’s message rather than dropped, so a flat line winning on AIC against two other flat lines is visible as that.

All three shapes are fitted by least squares rather than through nlmixr2, so all three together cost no measurable time, and the only compiled fits in the whole call are the PK candidates. model_report() breaks the wait down by fit: each concentration’s block carries what that fit took, one line per candidate where the search compiled more than one. A search that spent nine of its ten minutes in a second candidate is a run whose answer is to name pk.

Binary and ordinal endpoints are not fitted at all. They are drawn from one set of level frequencies over the whole study, as the shapes are one time course, or one set per arm under pd_by_arm = TRUE, so the response does not change over the visits. A level fewer than min_category_patients patients held is folded into the endpoint’s most common level first, as a rare categorical covariate level is excluded; an endpoint at which no level is held by that many patients is not drawn at all.

Step 4: Fit the Dosing Model and the Visit Model

Everything that is not the concentration-time curve comes from the models synpmx_pca_summarize() already builds, unchanged in what they represent. Dose reductions, interrupted cycles, discontinuation and missed visits are modelled here exactly as they are there, along with arm sizes, the censoring boundary and the schema. Baseline covariates use the study-wide distributions described below. The two generators call the same code rather than a copy of it; what differs is only the adapter that hands it the grid, and this generator writes its own over the nominal times the source holds.

The planned schedule is per arm, because it is the protocol. Everything read from how patients departed from it, and the visit model, is pooled over the arms, as the PD shapes are:

Model	What it holds	Drawn at generation as
Planned schedule, per arm	The nominal dose times enough of the arm reached, and the modal amount and infusion rate at each cycle among patients still on their starting dose. Where no amount is shared by min_arm_patients patients, as on a study dosed per kilogram, the arm’s mean to two significant figures	The cycle grid every subject starts from
Dose ladder	The levels patients dropped to, as ratios to their own starting dose to one significant figure, built from within-patient decreases: every level min_arm_patients patients of some arm share, one ladder for the study	The amount multiplier in force at a cycle
Reduction rate	Discrete-time hazard of stepping down a level, pooled over arms	Decided before the cycle is dosed
Interruption rate	Discrete-time hazard of skipping a cycle without ending treatment, pooled over arms	Decided at the cycle
Discontinuation rate	Discrete-time hazard of stopping treatment, pooled over arms	Decided after the cycle is dosed
Visit model	The visits each arm has, where at least min_arm_patients of its patients were observed, and one attendance rate per endpoint: the share of that endpoint’s scheduled visits at which a patient was observed	Attendance, drawn per visit at its endpoint’s rate

Pooling puts every patient behind each rate. Its cost is that an arm whose patients reduced, stopped or missed visits more often than the others’ takes the study’s rates, so a dose-dependent discontinuation is not reproduced, and that attendance no longer depends on when in the study a visit falls: a late visit few patients reached is generated as often as an early visit of the same endpoint, though each endpoint’s number of observations per patient is kept on average. An arm that never has a visit, such as a placebo arm without drug sampling, still has none.

A study where nobody reduces, skips or stops early has all three rates at zero and one level, and the model then reproduces the planned schedule exactly. No detector decides which kind of study this is, because the rates already say it.

Every arm, visit, rate and level rests on at least min_arm_patients patients on each side, 3 by default. A grid cell fewer patients were observed at is dropped, a dose-change rate whose events fewer patients had, across every arm, is zero, and an endpoint’s attendance rate is 1 where fewer patients missed any of its visits. Each arm also carries the values of its strata and keep columns, copied from the arm’s first patient; a keep value fewer than min_category_patients of the arm’s patients hold is written as missing, with a warning. Why both sides of a rate are floored is in the privacy article.

Baseline covariates are drawn independently from study-wide distributions. Each patient’s first row contributes one value, after small arms are excluded.

A continuous covariate is summarized by a mean and an SD: of its logarithm where every value is positive, of the value otherwise. Both are taken after the highest and lowest 5% of patients are set aside, at least one at each end, and the SD is divided by the SD a standard normal keeps under the same trimming, since the tails are what was removed. A cohort too small to leave three values after trimming is summarized whole.

For factor, character and logical covariates, levels held by fewer than min_category_patients patients (default 3) are excluded. The remaining counts are renormalized to probabilities, and every synthetic subject draws from that distribution. Set min_category_patients = 1 to retain all observed levels. Numeric category codes must be supplied as factors; numeric covariates keep their continuous distribution models. Binary and ordinal endpoints are folded by the same threshold at each visit (Step 3).

If no categorical level remains, generation keeps the column with missing values and estimation warns. Excluded factor labels are also removed from the stored schema. This changes the category proportions and can remove a covariate entirely from analysis. The threshold is applied during estimation, so stored fits must be re-estimated to use it.

Step 5: Store the Fit

Estimated by nlmixr2: the selected structural model, the candidate comparison table, fixed effects, the between-subject covariance matrix, the residual error, the covariate effects that survived, the PD shape and parameters per endpoint, and how long the fitting and the whole call took.

Not estimated: the dosing model, a planned schedule per arm with three rates pooled over arms, the visit model, arm sizes, the covariate distributions, the censoring boundary, the schema and the roles.

Only the fingerprint leaves the study. The fit also carries diagnostics for whoever ran it, such as the AIC of every candidate, the starting values and the correlations between covariates and random effects, and generation reads none of them. model_fingerprint() returns the part generation reads, and synpmx_model_generate() attaches that to its output rather than the fit. Printing the fit, or model_report(), gives the same account; vignette("pmxmodel-fingerprint") lists every number in a fingerprint.

Four operations protect the fingerprint as it is stored, each set out with its evidence in the privacy article:

1. Rounding. Every released PK and PD estimate is rounded to one significant figure, and the covariate summaries to two. The acceptance checks of Step 3 run on the fit before rounding, and the full-precision PK estimates stay on the fit, outside the fingerprint.
2. Single-patient influence. Estimation reads how far leaving out each patient would move each released estimate. A patient who moves one by 15 or more (in percent for a typical value, in points of between-subject SD on the log scale for a spread, and also by one standard error for a PD shape parameter) is left out of the PK fit, the PD fits and the covariate summaries, and estimation runs again. The patient stays in the dosing, visit and arm models. Rounds repeat until nobody is flagged, and stop before a tenth of the cohort has been left out; a reading still at 15 is a review and at 30 a failure. drop_influential = FALSE keeps every patient and only reports the reading.
3. No individual estimates. Each patient’s empirical Bayes estimates (EBEs) are read for the covariate correlations and the influence reading, then dropped. Generation draws random effects from the covariance matrix instead.
4. Rare levels. The schema stores no factor or discrete-endpoint level fewer than min_category_patients patients hold. An inferred binary or ordinal type left with fewer than two levels is stored as whole numbers.

model_privacy_checks() runs five checks on the result.

model_report(fit)
#> Summarized from the source, not estimated
#>   cohort             32 patients in 1 arm(s)
#>                      all (32)
#>   dose changes       none
#>   visit grid         2 endpoint(s) at 16 nominal time(s), 22 slot(s) in all
#>   visit attendance   share of scheduled slots attended, one rate per
#>                      endpoint: cp 55%, pca 91%
#>   covariates         wt lognormal, age lognormal, sex categorical, drawn
#>                      once for the whole study, independently of the
#>                      profiles
#>   columns emitted    id, time, ntime, dv, amt, evid, dvid, wt, age, sex
#> 
#> Values at the lower limit of what was observed
#>     cp                 0.2, no assay limit: half the lowest value several patients reached, rounded down
#>     pca                5, no assay limit: half the lowest value several patients reached, rounded down
#> 
#> Each non-PK continuous endpoint, fitted as constant, linear, or exponential
#>   pca                exponential
#>                        plateau          30
#>                        baseline         100
#>                        rate             0.1
#>                        between-subject  0.1 (SD on the log baseline)
#>                        residual         additive 10
#>                        chosen on AIC from constant, linear, exponential
#> 
#> PK endpoint for the PopPK model
#>   cp                 inferred from the following data characteristics:
#>                        absent before the first dose
#>                        rises to one peak and comes back down within one
#>                          dose interval
#>   route              oral: the median profile rises to a peak at 9 before
#>                      declining, and 31% of subjects do too
#>   sampling           median 6 distinct times after a dose, 6 after the
#>                      peak: sparse within-subject sampling
#> 
#> The PopPK model
#> 
#> Estimated by nlmixr2
#>   structural model   2cmt_oral 
#>   selected by        two-compartment model passed acceptance checks
#>   fitted on          all 32 patients with a concentration
#>   fixed effects      cl 0.1, v 7, q 0.1, v2 2, ka 0.4 
#>   between-subject    cl 0.265, v 0.2, ka 0.548, q 0.0775, v2 0.632 (as SD on the log scale)
#>   residual error     proportional 0.2 
#>   time to fit        50.9 s
#>   whole call         51.1 s, against 50.9 s in the fitter
#> 
#> Privacy
#>   one patient's pull pass: largest: cp: ka between-subject SD, 9.84 points;
#>                      `model_privacy_checks()` has every check

Step 6: Generate New Subjects

Per synthetic subject: assign an arm keeping the source arm shares, draw covariates from the study-wide covariate model, draw random effects from the covariance matrix, apply the covariate effects to the typical parameters, draw the dose schedule from the arm’s dosing model, draw the visits attended from the visit model, evaluate the profile at the attended times against the drawn schedule, add residual error, apply the censoring boundary and emit.

The order of the two draws is fixed. The schedule is drawn first and the profile is computed from it, so a synthetic subject who steps down a dose level has a lower exposure from that cycle on, and one who skips a cycle has the trough that implies. This is the one thing the generator does that neither synpmx_avatar() nor synpmx_pca() does. Drawing a profile and then a schedule would put a subject’s dose reduction in the dosing records and not in the concentrations.

The proportional multiplier is lognormal. Written as 1 + N(0, cv) it is non-positive with probability pnorm(-1 / cv), which is 2.3% of draws at a coefficient of variation of 0.5 and 7.7% at 0.7, and each of those is a concentration clamped to zero. A coefficient of variation that high is a statement that the structural model does not describe the data, and it belongs in the output as a wide band rather than as a scatter of zeros on the floor of a log axis. exp(N(0, sqrt(log(1 + cv^2)))) has the same coefficient of variation and a median of one, so the spread the fit estimated is preserved and no draw reaches zero.

A value is floored near the smallest the study reported. The residual is not the only thing that can put a synthetic value below anything the assay could return: a one-compartment profile evaluated late in a long dose interval underflows on its own. A study that declares a censoring column says where its assay stopped, and that boundary is put back. A study that declares none still had an assay, and its lowest values are the only evidence of where the limit sat. The floor is half the lowest value that at least min_arm_patients patients reached, each patient’s own smallest value taken first, rounded down to 1, 2 or 5 times a power of ten, the way assay limits are written. It is emitted wherever a draw falls below it, so it is not half the single smallest value, which is one patient’s measurement. An endpoint that reports a zero is given no floor, since a zero is a value a floor would contradict. model_report() carries the floor per endpoint.

A floor that catches a large share of the output is a warning, not a repair. Raising a value is still changing it, and a floor doing that to much of a dataset is hiding a fitted model that does not describe the low end of the data. What the floor caught is recorded on the generated dataset as the pmx_floored attribute, and above one observation in twenty the generator says so.

Gates

Only a study the generator cannot fit at all is refused. A study it can fit badly is fitted, and told about.

Gate	Threshold	Result
Cohort size	20 subjects	Warns and fits. A covariance matrix fitted to a handful of subjects describes those subjects, and nothing downstream will say so: the scorecard asks whether the output copies anybody or changed the study’s shape, and a small-cohort fit does neither. synpmx_pca()’s floor is 10 and this one is higher because a parameter estimate concentrates on its cohort faster than a score does.
Fit cap	max_fit_subjects = 60	The PK model is fitted to a subset drawn in proportion to the arms; the dosing, visit and covariate models read every subject. Reported in model_report(). Cannot be set below min_subjects.
PK cohort time coverage	6 distinct nominal times after a dose	Warns and fits. Below it a one-compartment model is not identifiable and the parameters sit close to their starting values, which is the caller’s call to accept.
No observation after a dose	PK endpoints only	Errors, naming which of the role columns the empty count came from: nothing selected as an observation, nothing selected as a dose, or the two never meeting in one subject.
Arm size	3 patients	Warns and drops those patients before anything is fitted, so the arm is absent from the model and from the data generated from it. Inherited from the dosing and visit models, which are summaries of an arm: an arm of one or two has no rates to pool.
nominal_time undeclared	—	Errors. The grid is a statement about the protocol only the caller can make.
Rare categorical level	min_category_patients = 3	Exclude levels held by fewer distinct patients; renormalize remaining counts. If none remain, warn and generate missing values. Set to 1 to retain all observed levels.
Rare discrete-endpoint level	min_category_patients = 3 over the whole study, or per arm under pd_by_arm	Fold the level into the endpoint’s most common level. An endpoint with no level held by that many patients is not drawn.
Single-patient influence	a move of 15 or more, in % or points; for a PD shape parameter, also one standard error or more	Leaves the patient out of the PK, PD and covariate estimates and estimates again, at most a tenth of the cohort, naming them on the console. A move still at 15 is a review and at 30 a failure that warns; drop_influential = FALSE skips the leaving out. Read by model_privacy_checks() and scorecard row E3; nothing per patient is stored.
Thin dose-change rate	fewer than min_arm_patients patients with the event, over every arm	The rate, pooled over arms, is zero.
Attendance rate	fewer than min_arm_patients patients who missed any of an endpoint’s visits; fewer than min_arm_patients of an arm’s patients at a visit	That endpoint’s rate is 1. An arm with too few of its patients at a visit does not have it.
Rare kept value	a keep value of any type fewer than min_category_patients patients in the arm hold	Written as missing, with a warning. strata columns are carried unchanged.
Rare level in the schema	a factor column’s level, or a binary or ordinal endpoint’s level, fewer than min_category_patients patients hold	Not stored. A column copied onto arms keeps the values its arms carry. An inferred binary or ordinal type left with fewer than two levels is stored as whole numbers.
No grid cell shared	min_arm_patients	The cell is dropped. A nominal time one patient attended is that patient.
Administration column	adm and routes together	Errors on either alone. What an administration id means is a convention of the dataset, and reading it wrong routes every dose to the wrong compartment silently.
PD shape candidacy	1 residual degree of freedom	The shape is dropped from the comparison, not the endpoint from the study. Below every candidate the endpoint is generated as a constant at its mean. A shape fitted exactly through its own points has AIC -Inf and would win any comparison it entered.
No PK endpoint identified	—	All continuous endpoints fail the baseline signal, or there are only discrete endpoints: infer PD-only and skip PK fitting. Other unresolved classifications error and name endpoint_roles.
Optimizer convergence	Successful termination or false convergence alone; finite AIC for FOCEi	False convergence warns and leaves acceptance to the remaining checks; other failures reject the candidate. SAEM uses the trajectory screen described in Step 3.
Parameter validity	Positive finite fixed effects; non-negative finite residual; usable covariance	Rejects the candidate.
Generation validity	Finite parameters and concentrations; numerical negatives below -1e-8 rejected	Rejects the candidate, using new population draws on the recorded designs.
Gross generation discrepancy	Median more than 100-fold outside source median bounds, or over 10% of profiles above 100 times their source-design maximum	Rejects the candidate. Censoring supplies bounds; unbounded source maxima omit the affected designs from the upper-tail check. Omissions are reported.
Moderate generation discrepancy	Median more than twofold outside source median bounds; uncensored interquartile-range ratio outside 0.5–2; median typical-profile discrepancy over fourfold on positive uncensored observations	Warns; does not trigger fallback. Censoring omits the spread comparison.
No candidate accepted	—	Errors with each attempted model’s rejection reason.
The between-subject terms did not move	1% of the eta init, fixed effects moved	A line under the parameter block in model_report(), not the banner: the fixed effects are estimates, but the spread synthetic subjects get is the starting value rather than this study’s.
Two endpoints tie as concentrations	neither tie-break separates them	Errors, naming both and endpoint_roles. Naming several is accepted, and each is fitted its own model.
The fit did not move	1% of every starting value	Warns, and prints !! THE FIT DID NOT MOVE !! above the parameter block in model_report(). A converged-looking fit is not necessarily an estimated one: where the optimizer takes no effective step, nlmixr2 returns an objective, an AIC and a full table of starting values, and nothing else downstream contradicts them — the generator simulates from them and the scorecard passes, because it asks whether the output copies anybody or changed the study’s shape and a fit that never moved does neither. Each parameter’s start, estimate and percent change are listed, because “did not move” is a claim the reader has to be able to check.
A fitted model as a public input	—	pmx_prior(), synpmx_prior() and synpmx_calibrated() refuse a pmx_fitted_model, because its parameters were estimated from the confidential study.