Why Make the Fitting Script Reproducible?
A joint-model result depends on more than the model formula. It also depends on which records entered the analysis, how variables were constructed, which software versions performed the fit, which numerical settings were used, and, for simulation-based procedures, how random numbers were generated. If those choices are scattered across an interactive session, a saved model object alone may not reveal how the result was obtained.
A reproducible fitting script makes the computational path inspectable: a reader can identify the inputs, follow each transformation, see the model and estimation settings, and rerun the fitting code in a sufficiently similar environment. Reproducibility does not mean that every computer will produce bit-for-bit identical output. Numerical libraries, software versions, parallel computation, and stochastic algorithms can affect the last digits or the random-number sequence. The goal is to record and control relevant choices, and to make any remaining variation visible.
The script should be treated as part of the analytical record, not as a convenient way to launch software. A saved object can preserve a fit for later inspection, but it cannot by itself establish which input files or processing decisions produced that fit.
What the Script Needs to Record
Keep the script in clearly separated stages. The exact function call depends on the joint-model software selected; the important practice is to make every consequential argument explicit and to save the call and its context.
Record the data-file paths, source or release, date or version where available, and any access or extraction details needed to locate the same data. Do not rely on an undocumented file already loaded in an interactive session.
Write out filtering, recoding, unit conversion, time-origin construction, and joining operations. Check the resulting identifiers, event coding, visit times, and missingness before fitting.
Record the longitudinal and event formulas, association structure, baseline-hazard choice, and any other modeling options. Do not rely on defaults when a default changes the scientific or computational specification.
Record optimizer, iteration limits, tolerances, quadrature settings, or MCMC settings that apply to the chosen method. Preserve the exact values passed to the fitting function.
Capture package and language versions, random-number settings when relevant, the fitted object, and a metadata record connecting the object to the code and input data.
For an R analysis, a script can record the environment with commands such as sessionInfo() and package versions with packageVersion("packageName"). A project lockfile can preserve package dependency versions for later installation. The lockfile and session record serve related but distinct purposes: the lockfile helps recreate dependencies, while the session record documents the environment used for a particular run.
When random-number generation is involved, set the seed before the operation that consumes random numbers and record the random-number generator settings as appropriate. For example, a Bayesian sampler may use a seed, whereas a deterministic maximum-likelihood fit may not use randomness at all. Do not claim that setting a seed makes an entire analysis reproducible if the algorithm uses parallel streams, nondeterministic hardware operations, or a different software implementation.
Save outputs to a known location and avoid silently overwriting the only copy of an earlier fit. A useful saved record includes the fitted object, the script or its version-control identifier, a data identifier, environment information, and a concise record of the fitting arguments. If the input data cannot legally be redistributed, the record should still explain how an authorized analyst can identify and obtain the required input.
A Minimal Script Structure
The following R-style outline shows how to make the stages visible. The final fitting call is deliberately described as package-specific: joint-model packages use different function names and argument conventions, so the actual call must match the selected package’s documentation and the model already intended for the analysis.
Inputs: set a project directory; read the longitudinal and event data from named files; retain file names and, where feasible, checksums or release identifiers as metadata.
Processing: apply documented inclusion rules; construct variables using explicit code; join the longitudinal and event records by the subject identifier; and check that the join did not unexpectedly duplicate or remove subjects.
Specification: store the longitudinal formula, event formula, association form, and baseline-hazard choice as readable text or explicit arguments. For example, a project might record “longitudinal outcome: biomarker as a function of time and treatment with subject-specific random intercept and time slope” and “event outcome: time to event with treatment as a covariate; current-value association.” Those descriptions must agree with the actual package call.
Fit: call the selected package’s documented fitting function with the prepared data and explicit estimation controls. Keep the call in the script rather than copying it into a console history that may not be retained.
Save: save the model object and a metadata object containing the script version, input identifiers, software versions, formulas, controls, and seed or generator settings if used. Confirm that the saved file can be read back and that the expected object is present.
Writing out this structure is more reliable than storing a single long command. It also gives each stage a natural place for validation checks and makes it easier to determine whether a mismatch arose during data preparation, model specification, or fitting.
Worked Example: Recording Data Processing
Suppose a hypothetical cohort has 12 subjects. The longitudinal table contains 36 visits, and the event table contains one row per subject. One subject has no event time recorded and is excluded under a prespecified rule requiring an observed follow-up time. Of the remaining 11 subjects, one has no longitudinal measurements and is also excluded because the planned joint fit requires biomarker observations.
Processing record. The script should state the two rules and apply them in order. After the event-time rule, \(12-1=11\) subjects remain. After the no-longitudinal-measurement rule, \(11-1=10\) subjects remain. If the 36 visits include 3 visits belonging to the subject excluded for missing event time and 2 visits belonging to the subject with no usable biomarker measurements, the final visit count is \(36-3-2=31\).
Checks. After processing, the script should confirm that the event table has 10 unique subject identifiers, that all 31 longitudinal rows refer to those 10 identifiers, and that each retained subject has the required event-time information. The arithmetic alone does not prove the data are correct; these checks test whether the implemented filters produced the intended records.
Metadata to save. Record the original row counts, the two exclusion rules, the retained subject and visit counts, and the input file identifiers. If a later rerun produces 9 subjects or 30 visits, the discrepancy is then detectable before interpreting a model output.
Making Model and Estimation Choices Auditable
Formula strings are helpful, but they are not enough if the fitting function also has arguments that alter the estimand or numerical procedure. Record both the scientific specification and the actual call. For instance, distinguish a current-value association from a slope association, and document whether the hazard uses a parametric, spline-based, or other supported baseline specification. These are substantive choices, not merely software settings.
Also save estimation controls in a form that can be compared with the fitting call. A value listed in a note but not passed to the function is not a controlled setting. Conversely, a default that is not recorded may change between software releases. When the package uses defaults that are unavoidable or not exposed, record the package version and identify that the setting was left at its default.
Worked Example: Checking That Controls Match the Fit
Consider a hypothetical maximum-likelihood fit configured to allow at most 800 optimizer iterations, with tolerance \(10^{-6}\), and 15 quadrature points. The script metadata says “800 iterations, tolerance \(10^{-6}\), 15 points,” but the actual call passes 500 iterations, tolerance \(10^{-5}\), and 15 points.
Compare each setting. For iterations, the recorded value is 800 and the supplied value is 500, so the discrepancy is \(800-500=300\) iterations. For tolerance, \(10^{-5}/10^{-6}=10\), so the supplied tolerance is ten times larger (less stringent) than the recorded value. The quadrature setting agrees: \(15=15\).
Conclusion. The fit was not run under the documented controls. The correct repair is to reconcile the script and metadata, rerun if the intended settings were not used, and save the corrected call with the resulting object. It would be misleading to describe the existing object as having used the documented settings.
This kind of comparison can be automated by creating one named control object and passing that same object to the fitting function and to the metadata record, where the software interface allows it. If the package accepts settings only as separate arguments, create the metadata from those same argument values rather than retyping them independently.
Random Seeds, Versions, and Saved Objects
A random seed initializes a pseudorandom-number sequence; it does not specify every condition needed to recreate that sequence. The same seed can yield different results if the algorithm, software version, random-number generator, execution order, or parallel configuration changes. Record the seed near the code that uses it, and record relevant generator or parallel settings when the procedure depends on them. For MCMC fits, also retain the chain configuration and diagnostics needed to interpret whether sampling was adequate; a seed is not evidence of convergence.
Version records matter because software updates can change numerical routines, defaults, bug fixes, or object structures. A strong record includes the statistical language version, the joint-model package version, and versions of important dependencies. A lockfile can help reconstruct the dependency environment, but it does not replace the script, input-data identification, or saved fitting arguments.
Save the fitted object in a format that preserves the object for the software used to create it, and test the save-and-read step. Include a metadata record that identifies the fit unambiguously. For example, the object name should not merely be “final”; it should be associated with a date or run identifier, the candidate specification, and the script revision. If a fit is later compared with another candidate using likelihood-based criteria, preserving the exact model calls and environment helps establish that the comparison is interpretable.
Diagnostics, Alternatives, and Limitations
Reproducibility checks begin before fitting. Look for unexpected counts after filtering and joins, duplicate subject identifiers in a table meant to contain one row per subject, missing required fields, impossible time values, and discrepancies between the recorded formulas and the actual call. After fitting, record warnings, convergence information, and relevant diagnostics alongside the object; do not suppress warnings merely to make a script run quietly. Numerical warnings and failed fits require substantive review, not just a more convenient saved file.
Version control is a useful complement to a script: it records changes to code over time and allows a run to be connected to a particular revision. A project lockfile complements it by recording dependency versions. Neither replaces a well-documented data-processing pipeline, and neither guarantees that unavailable or changing source data can be recovered. A workflow platform or container can further control the computational environment, but it adds maintenance and does not correct an incorrect model specification.
Interactive analysis can be valuable for exploration, but it is a weaker final record when the steps depend on hidden objects, manual edits, or an untracked execution order. A script is more inspectable, while a saved model object is more convenient for later summaries. The two should be used together. Even then, exact numerical identity can be limited by floating-point arithmetic, hardware, compiler, or parallel execution. Report that limitation rather than implying stronger reproducibility than the environment supports.
Worked Example: Auditing a Saved Fit
A hypothetical analyst saves a fit as “joint_model.rds.” The accompanying note reports a current-value association and 10 quadrature points. The script instead specifies a slope association and 12 points. The saved data-processing record says 10 subjects were analyzed, while a check printed at fitting time reports 9.
Audit the discrepancies. The association form differs: current value is not the same specification as slope. The quadrature count differs by \(12-10=2\) points. The subject count differs by \(10-9=1\) subject. Each mismatch is evidence that the note, script, or executed run does not fully describe the same analysis.
Resolution. Do not choose whichever version seems most plausible after seeing the result. Inspect the executed script and logs, determine which data and arguments produced the saved object, correct the metadata, and rerun if the intended specification cannot be established. Save the corrected object with a run identifier and preserve the original artifact until the discrepancy has been resolved.
Check Your Understanding
Answer each question, then compare with the solutions below.
- Why does setting a random seed not by itself guarantee identical joint-model results on every computer?
- A data-processing script starts with 18 subjects, excludes 2 under a stated rule, then excludes 1 more for having no longitudinal measurements. How many subjects should the final event table contain?
- The metadata record says 20 quadrature points, but the fitting call uses 16. What should be done before reporting the fit as using 20 points?
- Name two records, in addition to the saved model object, that help another analyst rerun a fit.
- What is one limitation of a package lockfile?
Solutions. (1) The algorithm, software version, generator, execution order, parallel settings, or hardware can differ, so the same seed may not reproduce the same sequence or result. (2) \(18-2-1=15\) subjects. (3) Reconcile the record with the executed call; if 16 points were used, correct the metadata, and if 20 were intended, rerun with 20 and save the resulting object and call. (4) Examples include input-data identifiers, the processing script, the model formulas and controls, software versions, and random-number settings. (5) A lockfile can help recreate dependencies, but it does not identify or preserve the input data, document processing, or guarantee identical numerical results across hardware and software environments.