Simulation-Based Parameter Recovery

Why Simulate Parameter Recovery?

A joint-model fit can converge and still estimate a parameter poorly. When information is weak, estimates may be biased, intervals may miss their target too often, or an optimizer may fail to find a usable solution. A parameter-recovery simulation makes these concerns concrete: specify a model and known parameter values, generate many hypothetical data sets from them, refit the model, and compare the estimates with the values used to generate the data.

The purpose is not to show that a model is correct for an actual study. Rather, it is to assess how an estimator and fitting procedure behave under specified conditions. Recovery depends on the whole setup: the assumed data-generating model, sample size, visit schedule, event rate, censoring, parameter values, and numerical fitting procedure. Good recovery in one scenario does not establish good recovery in every scenario.

Specify the Model Before Simulating

Consider a shared-random-effects joint model with one continuous longitudinal marker and a time-to-event outcome. For subject \(i\), let \(Y_{ij}\) be the marker measured at time \(t_{ij}\), and let \(T_i\) be the event time. One possible model is

The latent marker trajectory is \(m_i(t)=\beta_0+\beta_1t+b_{0i}+b_{1i}t\). In this example, the event hazard conditional on the random effects and a baseline covariate \(X_i\) is

Here \(h_0(t)\) is the baseline hazard, \(\gamma\) is the covariate log-hazard coefficient, and \(\alpha\) is the association parameter linking the latent marker to the hazard. A one-unit increase in \(m_i(t)\) corresponds to a hazard ratio of \(\exp(\alpha)\), conditional on the covariate and the model. This is a conditional model statement, not a causal conclusion.

A simulation should state its assumptions explicitly. In this example, subjects are independent; random effects follow the specified bivariate normal distribution; residual errors are independent normal variables with variance \(\sigma^2\), independent of the random effects; and the conditional event hazard follows the stated model. Visits may be scheduled at fixed times, with marker measurements stopping at the event or censoring time. For a simple design, censoring can be administrative at a fixed study end \(\tau\). The observed time is then \(U_i=\min(T_i,\tau)\), with event indicator \(\delta_i=I(T_i\leq\tau)\). If censoring is random instead, its generation and assumptions should also be specified.

The fitted model should match the data-generating model when the aim is to assess recovery under correct specification. If the simulation instead deliberately generates data under a different model, that is a misspecification study and should be described as such. Before interpreting a poor result as weak identification, check that the generating and fitting code use the same parameter definitions, time scale, covariate coding, baseline hazard, and censoring convention.

Designing the Repeated-Sampling Experiment

Choose the target parameters and values before generating data. Select study features that matter for information: number of subjects, visit times, marker measurement error, event frequency, censoring, and the random-effects covariance. For a shared-parameter question, vary \(\alpha\) across scenarios while keeping the other features controlled. If the random-slope variance is near zero, for example, the data may carry little information about the slope-related part of the trajectory and its contribution to the hazard.

For each replicate, generate subjects, marker observations, event and censoring times according to the specified design, then fit the prespecified model. Record the estimate, its standard error or interval, convergence status, and any predeclared diagnostics. A convergence flag alone is not enough to define a usable fit. The rule might also require a finite estimate, a finite uncertainty estimate, and a gradient norm below a stated threshold. Apply the same rule in every scenario, and report it.

Recovery Measures and Their Interpretation

Suppose \(R\) replicates are generated, the true value of a scalar parameter is \(\theta\), and \(\hat{\theta}_r\) is its estimate in replicate \(r\). If all \(R\) fits are usable, the empirical bias is the average estimation error:

A positive value means estimates tend to exceed the generating value in this simulation; a negative value means they tend to fall below it. Bias near zero does not guarantee useful estimation: large variability can still make estimates unreliable. Report the empirical standard deviation of the estimates as well, and distinguish it from the average model-based standard error.

For nominal \(100(1-\alpha)\%\) intervals \([L_r,U_r]\), empirical coverage is the proportion of intervals that include the true value:

Coverage near the nominal level is desirable over repeated samples under the specified scenario. It does not mean that a particular reported interval has that probability of containing a fixed parameter. Also report interval width: a procedure can obtain high coverage simply by producing very wide intervals.

Let \(F\) of the \(R\) fitting attempts meet the prespecified failure rule. The empirical optimizer failure rate is \(F/R\). State what counts as failure and whether failed fits are excluded from the bias and coverage calculations. Excluding them can make the reported summaries look better than the procedure as a whole. Report failure rates separately, and, when reporting conditional summaries among successful fits, label them clearly.

Worked Examples

A Short Justification for Monte Carlo Precision

Simulation summaries are themselves estimates because only finitely many replicates are run. For coverage, define \(C_r=I(L_r\leq\theta\leq U_r)\), which equals 1 for coverage and 0 otherwise. If replicates are generated independently under the same scenario and the probability of coverage is \(q\), then \(C_r\) has variance \(q(1-q)\). The sample average of \(R\) such indicators has variance \(q(1-q)/R\), giving the estimated Monte Carlo standard error

This formula explains why a small simulation can give a noisy coverage estimate even when the fitting procedure is unchanged. For estimated bias, if the replicate estimates have empirical standard deviation \(s_{\hat{\theta}}\), the Monte Carlo standard error of their mean error is \(s_{\hat{\theta}}/\sqrt{R}\), provided the same set of replicates is used. Increasing \(R\) reduces Monte Carlo uncertainty; it does not fix a flawed model or an inadequate study design.

Reading Patterns, Comparisons, and Failure Modes

A useful recovery study compares a small number of planned scenarios rather than reporting one isolated result. For example, compare a larger and smaller sample size, or a stronger and weaker association parameter, while holding other design features fixed. If weak association scenarios show more variable \(\hat{\alpha}\), poor coverage, or more optimizer failures, that pattern is evidence that this model-and-design combination provides less reliable estimation of \(\alpha\). It is not, by itself, proof that no estimator can identify the parameter.

Parameter-recovery simulation complements, rather than replaces, diagnostics from fitting the observed data. Information-matrix eigenvalues and Hessian-based uncertainty estimates describe local features of a fitted likelihood or posterior approximation. Recovery evaluates repeated behavior under the chosen simulation scenarios. Disagreement between these views is a reason to investigate, not a reason to select whichever summary is more reassuring. For Bayesian fitting, define in advance whether intervals are credible intervals and assess the relevant posterior summaries across replicated data sets; do not treat optimizer failure criteria as a substitute for MCMC diagnostics.

Several practical checks help distinguish genuine weak recovery from a simulation or implementation error:

• Verify generation. Check that simulated event and marker data have the intended scales, event frequency, censoring, and random-effects variation.
• Verify parameter mapping. Confirm that a covariance parameter, standard deviation, and association coefficient mean the same thing in the generator and fitted model.
• Use stable fitting settings. Keep scaling, integration accuracy, starting-value strategy, and stopping rules consistent with the intended analysis.
• Inspect more than averages. Plot or tabulate estimates, interval widths, boundary estimates, and failures; a small mean bias can hide extreme or multimodal estimates.
• Increase replication when needed. If Monte Carlo standard errors are large relative to the differences between scenarios, the simulation does not yet resolve those differences.

The main alternatives answer different questions. An information-matrix or Hessian diagnostic examines local information at a fitted solution; it is faster but can miss finite-sample bias or failure to reach that solution. A recovery simulation directly examines repeated-sample behavior under specified conditions, but its conclusions are only as relevant as its scenarios and its data generator. A larger sample size may improve some forms of recovery, but it cannot automatically resolve model misspecification, a poorly chosen time scale, or a parameter that the specified design barely informs.