arXiv CSopen access

Overcoming Selection Bias in Statistical Studies With Amortized Bayesian Inference

2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

machine learning, deep learning, neural networks

Overcoming Selection Bias in Statistical Studies With Amortized Bayesian Inference Jonas Arruda1,2 , Sophie Chervet3,4 , Paula Staudt5,6 , Andreas Wieser7,8,9,10 , Michael Hoelscher7,8,9,11 , Isabelle Sermet-Gaudelus12,13,14 , Nadine Binder6,15 , Lulla Opatowski3,4 , and Jan Hasenauer∗1,2

arXiv:2604.18319v1 [stat.ML] 20 Apr 2026

Bonn Center for Mathematical Life Sciences, University of Bonn, Bonn, Germany 2 Life & Medical Sciences Institute, University of Bonn, Bonn, Germany 3 Epidemiology and Modeling of Antibiotic Evasion Unit, Institut Pasteur, Paris, France 4 Université de Versailles Saint-Quentin-en-Yvelines, Université Paris Saclay, Inserm U1018, Team Infectious Diseases, Interactions and Antimicrobial Resistance, Paris, France 5 Institute of Medical Biometry and Statistics, Faculty of Medicine and Medical Center, University of Freiburg, Freiburg, Germany 6 Freiburg Center for Data Analysis, Modeling and AI, University of Freiburg, Freiburg, Germany 7 Institute of Infectious Diseases and Tropical Medicine, LMU University Hospital, Munich, Germany 8 German Center for Infection Research, Partner Site Munich, Munich, Germany 9 Fraunhofer Institute ITMP, Immunology, Infection and Pandemic Research, Munich, Germany 10 Max von Pettenkofer Institute, LMU Munich, Munich, Germany 11 Unit Global Health, Helmholtz Zentrum München, German Research Center for Environmental Health (HMGU), Neuherberg, Germany 12 Centre de Référence Maladies Rares, Mucoviscidose et Maladies Apparentées, Site Constitutif Pédiatrique, Hôpital Necker Enfants Malades, Paris, France 13 Université de Paris, CNRS, INSERM, Institut Necker-Enfants Malades, Paris, France 14 European Rare Disease Network–Lung, Frankfurt, Germany 15 Institute of General Practice/Family Medicine, Faculty of Medicine and Medical Center, University of Freiburg, Freiburg, Germany

April 21, 2026

Abstract Selection bias arises when the probability that an observation enters a dataset depends on variables related to the quantities of interest, leading to systematic distortions in estimation and uncertainty quantification. For example, in epidemiological or survey settings, individuals with certain outcomes may be more likely to be included, resulting in biased prevalence estimates with potentially substantial downstream impact. Classical corrections, such as inverse-probability weighting or explicit likelihoodbased models of the selection process, rely on tractable likelihoods, which limits their applicability in complex stochastic models with latent dynamics or high-dimensional structure. Simulation-based inference enables Bayesian analysis without tractable likelihoods but typically assumes missingness at random and thus fails when selection depends on unobserved outcomes or covariates. Here, we develop a bias-aware simulation-based inference framework that explicitly incorporates selection into neural posterior estimation. By embedding the selection mechanism directly into the generative simulator, the approach enables amortized Bayesian inference without requiring tractable likelihoods. This recasting of selection bias as part of the simulation process allows us to both obtain debiased estimates and explicitly test for the presence of bias. The framework integrates diagnostics to detect discrepancies between simulated and observed data and to assess posterior calibration. The method recovers well-calibrated posterior distributions across three statistical applications with diverse selection mechanisms, including settings in which likelihood-based approaches yield biased estimates. These results recast the correction of selection bias as a simulation problem and establish simulation-based inference as a practical and testable strategy for parameter estimation under selection bias. ∗ Corresponding author: [email protected]

Introduction Population-level inference is central across many disciplines, including epidemiology [Rothman et al., 2008], pharmacology [Strom, 2019], medicine [Forrest et al., 2025], and economics [Heckman, 1979]. However, in many such settings, the ability to generalize to the overall population cannot be guaranteed due to structural features of study design or data collection [Rudolph et al., 2023]. For instance, if women are less likely to participate in a seroprevalence study but are more likely to be infected, the sample will overrepresent men, and the estimated prevalence will be biased downward. Such lack of representativeness is commonly referred to as selection bias and arises when the probability that an observation is included in a dataset depends on characteristics related to the quantities of interest, such as observed covariates, unmeasured confounders, outcome-dependent sampling, informative censoring, or latent variables that jointly influence inclusion in the study and outcomes [Cochran, 1977, Kleinbaum et al., 1981, Rothman et al., 2008, Smith, 2020]. When selection bias is ignored, analyses can yield systematically biased inference results [Little and Rubin, 2019]. At the same time, multiple sources of bias are rarely modeled jointly, and selection bias is often neglected relative to other sources such as misclassification or uncontrolled confounding [Petersen et al., 2021]. With modern large-scale studies, where random errors diminish as sample size increases, systematic distortions arising from selection and study design can become the dominant source of estimation errors [Kaplan et al., 2014]. Most approaches to mitigating selection bias rely on explicitly modeling the data-generating and sampling processes under assumptions that yield tractable likelihoods or correction formulas. Methods such as inverse probability weighting [Horvitz and Thompson, 1952], the Heckman correction [Heckman, 1979], propensity score techniques [Austin, 2011], and instrumental variables approaches [Streeter et al., 2017] all require parametric structures or tractable representations of the selection mechanism to adjust estimates. Bayesian formulations extend these paradigms by jointly modeling outcome and selection processes within a probabilistic framework [Copas and Li, 1997, Scharfstein et al., 2003, Linero and Daniels, 2018, Kawabata et al., 2024]. However, they likewise depend on tractable likelihoods for both the outcome and the selection process. These requirements impose restrictive parametric assumptions and limit applicability when the underlying processes involve latent stochastic dynamics, nonlinear interactions, or high-dimensional hidden states that preclude tractable computation. As statistical models become more complex, the necessity of tractable likelihoods becomes a central bottleneck, rendering bias correction infeasible. Simulation-based inference (SBI) has emerged as a general framework for Bayesian inference in models with intractable likelihoods by approximating posterior distributions only using simulations generated by a statistical model. A broad class of approaches relies on neural density estimation to approximate the posterior distribution, including neural posterior estimation, neural likelihood estimation, and neural ratio estimation, and has been shown to be effective across various different applications [Lueckmann et al., 2021]. These methods can be applied either in a dataset-specific manner or in an amortized setting. In amortized Bayesian inference, a generative neural network called a neural posterior estimator (NPE) is trained once on simulated data and can then be reused across arbitrarily many datasets at test time without retraining. Because inference relies only on the ability to simulate from the model, SBI methods have been developed to work with complex, high-dimensional, and implicit generative processes [Arruda et al., 2025]. Recent extensions have improved the robustness of SBI to practical challenges such as outliers [Schälte et al., 2021, Bharti et al., 2026, Khoo et al., 2026], domain shifts [Elsemüller et al., 2025], and certain forms of missing data [Wang et al., 2024, Gloeckler et al., 2024]. Despite these advances, SBI methods generally assume that the observed data constitute unbiased realizations of the modeled generative process. In practice, however, the data available for inference often result from structured selection mechanisms. When such mechanisms are not explicitly represented in the simulator, applying SBI to selectively observed data can lead to systematically biased posterior estimates and miscalibrated uncertainty. Detecting such violations is itself challenging, as the true parameters are typically unknown and standard posterior predictive checks and coverage tests cannot reveal biases induced by the selection mechanism. To date, no amortized framework explicitly accounts for outcome-dependent or dynamically coupled selection mechanisms while simultaneously providing tests to detect bias and assess posterior calibration, leaving bias-aware likelihood-free inference under structured selection an open methodological challenge. To address this limitation, we develop a bias-aware amortized Bayesian inference framework that embeds the selection mechanism directly within the generative model. Selection, censoring, and missingness are

S Selection

Bias depends on

e.g., study with a nonrepresentative sample of the population

C Covariates

θ Y S forward process inverse process

& e.g., selection depends on individuals still being observable, i.e., not dead

θ Model parameters &

e.g., only infected indiviudals are included in a study

Y Outcomes

selection process

Figure 1: Different sources of selection bias in statistical studies. The forward process generates outcomes from model parameters and covariates, with selection acting as an additional component that determines which individuals are observed. Selection may depend on covariates alone (e.g., a non-representative sample), on latent states such as survival and hence parameters (e.g., individuals must be alive to remain observable), or on outcomes (e.g., only infected individuals are enrolled). The inverse process must account for the selection in the forward process for bias-aware inference. formalized as structured observation operators acting on latent population-level processes, such that the population dynamics and the bias-inducing selection process are jointly simulated and NPEs are trained to infer population parameters from biased data. This likelihood-free construction requires only forward simulations and therefore avoids the need for tractable likelihoods. Because the framework can be amortized across similar selection mechanisms, it enables inference under heterogeneous and complex study designs within a unified model. Simulating both the data-generating and selection mechanism further allows systematic assessment of posterior calibration using simulation-based calibration (SBC) [Talts et al., 2020] and classifier two-sample tests (C2ST) [Lopez-Paz and Oquab, 2017], which enable assessment of bias induced by a selection mechanism and detection of bias in real datasets. We evaluate the approach in three representative epidemiological settings reflecting distinct forms of structured selection: prevalence estimation under biased sampling with partially missing outcomes, incidence estimation under informative censoring due to death with coupled event and survival processes, and transmission modeling in which covariates and infection outcomes jointly influence study inclusion. Across these scenarios, the proposed approach yields calibrated inference even in regimes where classical corrections become unreliable due to model complexity or intractable likelihoods.

Results A Framework for Bias-Aware Amortized Bayesian Inference Under Selection To provide bias-aware inference under selection, we propose an amortized Bayesian inference framework that explicitly incorporates the selection mechanism into the generative model. The framework builds on SBI using neural posterior estimation and addresses the challenge that observational datasets may arise from selection mechanisms that filter the underlying population. We therefore formulate the inferential target as the conditional distribution p(θ | Y, S = 1; C), where outcomes Y arise from a model with parameters θ and covariates C, and inclusion is governed by a binary selection process S. The process S is an outcome-, covariate-, or parameter-dependent stochastic selection process (Figure 1). In this setting, the observed dataset constitutes a filtered realization of the underlying generative process. A key innovation of our approach is to embed the selection mechanism directly within the forward simulation used for inference. Parameters and covariates are sampled from the prior, outcomes are generated from a mechanistic population model, and the selection operator is subsequently applied to 3

Age

Proportion Proportion

0.2

Infection model

Infection parameters θ

0.0

Logistic model Missing completely at random

Covariates C

Age, sex, household size, birth country

19 20 34 35 49 50 64 65 79

0.5

Inverse Infections Y

Sample S

R3 0

R4 R5 R1 R2 Simulated round (with missingness)

Unadjusted

Female

Male

0.0

Household size KoCo19 R1 y missing Munich

R2 R3

Other Germany

Birth country R4 R5

0.0

Estimated Prevalence (%)

Absolute error (%)

0.5

0.0

Missing at random

Sex

0.5

5 0

Simulated round Weighted

KoCo19 Round Bias-aware NPE

Figure 2: Estimation of prevalence from biased population sample on simulated data and the KoCo19 Study. (A) Visualization of the infection model and selection bias based on covariates. (B) Covariate distribution in the KoCo19 cohort for all 5 rounds compared to Munich, with indication for whom the infection status y is missing. (C) Absolute error of estimated prevalence across 1000 simulated datasets (including selection and missingness), each with five rounds of the KoCo19 Study using unadjusted counts of infections, inverse-probability weighting, and the mode of the posterior estimated with bias-aware NPE. (D) Estimated prevalence across five rounds of the KoCo19 Study using unadjusted counts of infections, inverse probability weighting, and bias-aware NPE. For the unadjusted and weighted approach, we bootstrapped the estimate, and for NPE, we show the full posterior. obtain the observed sample. Training the NPE on retained observations yields simulations from the joint distribution p(θ, Y, S = 1, C), so that the learned posterior approximation directly targets the correct conditional distribution under the assumed selection mechanism. Because this construction relies only on forward simulations, it avoids the need for tractable likelihoods for either the outcome model or the selection process (see Methods). While the NPE is trained on simulated data, successful bias removal on real data needs to be validated. The amortized nature of the framework enables repeated fast inference and hence systematic validation of this aspect by repurposing established simulation-based diagnostics for the assessment of selection bias. First, SBC is used to assess whether posterior ranks of ground-truth parameters are uniformly distributed across repeated simulations from the joint model, thereby verifying calibration with respect to the selection mechanism. Second, a C2ST is integrated to evaluate whether samples drawn from the learned posterior together with the observed data are statistically indistinguishable from parameter-data pairs generated directly from the joint distribution during training. Together, these diagnostics formalize the criteria under which inference is considered calibrated and bias-aware. The trained bias-aware NPE defines a global mapping from observed datasets to parameter distributions, and, due to inference being amortized, it can be evaluated across heterogeneous selection regimes without retraining. This enables systematic analysis of bias, calibration, and robustness under distinct forms of selection. In the following sections, we train bias-aware NPEs on three representative classes of selection mechanisms and examine the framework’s inferential behavior in each setting.

Amortized Inference Reduces Selection Bias in an Incomplete Data Setting To assess the proposed framework, we first consider disease prevalence estimation under sampling bias and covariate and outcome missingness (Figure 2A). For this class of problems, several statistical methods are available for evaluating bias correction. Here, we consider data from the KoCo19 Study, which was conducted to estimate the seroprevalence of SARS-CoV-2 antibodies in the general population of Munich over time [Radon et al., 2020]. The study consists of 5 rounds with increasing missingness of outcomes (Figure 2B). We first evaluated the proposed framework in a controlled yet realistic simulation setting and subsequently applied it to the observed data. Ground-truth prevalence was obtained from a synthetic population matching Munich’s demographic structure, with infection outcomes simulated from a logistic model (see Methods). Non-representative samples were then generated by resampling according to the original study sample weights and by imposing the observed round-specific missingness patterns (Figure 2B). To assess the framework’s performance, we trained a NPE on such biased subsamples paired with their corresponding true prevalence and estimated prevalence with an unadjusted estimator and inverse probability weighting (see Methods). This bias-aware NPE learned to correct both for non-representative sampling and outcome missingness (Figure 2C). On 1000 simulated datasets, bias-aware NPE yielded more accurate prevalence estimates than both the unadjusted estimator and inverse probability weighting, while retaining similarly fast inference due to amortization. The latter two approaches relied here on complete case analysis and therefore exhibited bias and increased variance in the presence of missing outcomes and a reduced number of samples in the data (Figure 2B). When outcome missingness was removed from the simulation, the NPE recovered prevalence estimates that were indistinguishable from those obtained via inverse probability weighting, confirming that additional gains arise from explicitly accounting for missingness (Figure S1A). The application of the bias-aware NPE, the unadjusted estimator, and inverse probability weighting to the real study data revealed that for the first 2 rounds, the prevalence estimates were broadly consistent across the methods (Figure 2D). Here, the same amortized NPE trained exclusively on synthetic data was applied to all 5 study rounds without retraining. In contrast, for rounds 3 and 5, the bias-aware NPE produced estimates that differed systematically from inverse probability weighting, as reflected by the different shapes of the uncertainty distributions. In round 3, inverse probability weighting yielded higher median prevalence estimates (6.70% vs. 5.96%), whereas in round 5 it produced lower median estimates compared to the bias-aware NPE (11.19% vs. 11.73%). This discrepancy arises from the increasing degree of outcome missingness across rounds and suggests that inverse probability weighting alone is insufficient when missingness is informative (Figure 2C). The assessment of the reliability of the posterior using C2ST revealed that posterior samples generated by the bias-aware NPE were statistically indistinguishable from samples drawn from the joint simulator. Moreover, on simulated validation datasets, the classifier operated at chance level, which confirms convergence of the neural network and calibration with respect to the simulated selection mechanism. When applied to real study rounds, the classifier likewise performed at chance level (Figure S1B), indicating that posterior samples were consistent with the assumed generative model and selection mechanism. These results support that the inferred posteriors remain calibrated despite the presence of sampling bias and outcome missingness.

Amortized Inference Enables Calibrated and Testable Time-to-Event Analysis Under Missing Disease Information To further assess the proposed framework, we consider time-to-event analysis under missing disease information due to death. In this setting, the event of interest cannot be observed once death occurs, which can introduce substantial bias in incidence estimates [Leffondré et al., 2013, Binder and Schumacher, 2014]. Such mechanisms can be seen as selection of observed data and can arise in any competing risk model and are particularly relevant for neurodegenerative diseases, such as dementia, where mortality competes with disease onset. We therefore consider data from the long-term population-based Framingham Heart Study, which has been used to study temporal trends in dementia incidence [Satizabal et al., 2016]. We base the analysis on an established illness-death model (IDM) with states corresponding to healthy, dementia, and death, where death may obscure the transition from healthy to dementia (Figure 3A,B). We first evaluated the proposed framework in a controlled simulation setting and subsequently applied it to the

10 1

Hazard parameters θ 1 10

Healthy

100Death S

Dementia/death cumulative hazard per 100 persons

10 1

states

Initial states

Censored, observed only at visit

Intermediate

Age, sex

Dementia unknown 2

Hazard NRMSE

102 101 100

10 1

h01

Inverse

Dementia Y

Covariates C

Final states

Death cumulative hazard per 100 persons

Multi-state model

Death Dementia

Dementia/death cumulative Death cumulative Dementia cumulative hazard per 100 persons hazard per 100 persons hazard per 100 persons

Dementia cumulativ per 100 perso

Time-to-event model

100

Epoch 1

Epoch 2

Epoch 3

Epoch 4

100

10 1 101 100

10 1

101

100

10 1

1h 2 02

NPE (full data)

102 101 100

101 1 2 3 4 5 1 2 3 4 5 1 12 32 4 35 41 25 3 4 5 1 2 3 4 5 Follow up years since entry in epoch Follow up years since entry in epoch

5 h 1 12

NPE (observed data)

Bias-aware NPE

IDM

Figure 3: Correcting for bias due to missing disease information because of death on simulations and the Framingham Heart Study. (A) Visualization of the time-to-event model and bias due to missing dementia because of death. (B) Observable and unobservable transitions between states from over 2200 patients in the Framingham Heart Study summed over all 4 epochs. Parts of the transitions to death are interval-censored since last visit with unknown dementia status. (C) Normalized root mean squared error (NRMSE) of the recovered transition hazard (hkl : 0 → 1 dementia onset, 0 → 2 death without prior dementia, and 1 → 2 death after dementia) adjusted by covariates on simulated data with and without censoring. For each posterior sample the mean over time is computed and the error normalized by the mean of the corresponding ground truth hazard. Posterior samples are aggregated with the median, and the recovery error across datasets is shown. (D) Cumulative hazards adjusted by covariates estimated using real data. We compared our bias-aware NPE approach (median cumulative hazard and 95% credible intervals) against the NPE trained on full data and a penalized likelihood approach. For the latter, the maximum likelihood estimate is shown together with intervals obtained from the lower and upper intensity curves (interpreted as 95% bounds of the hazard function) provided by the underlying IDM. observed data. For the simulation study, we generated data by simulating trajectories from the IDM using covariates observed in the Framingham cohort (full data) and subsequently induced missing dementia information due to death (observed data) (see Methods). The simulated data were analyzed using two different approaches: (1) a standard NPE trained on the full data and (2) the proposed bias-aware NPE trained on the data as observed. We analyzed recovery accuracy using the normalized root mean squared error of the covariate-adjusted transition hazards. The assessment of the two methods using simulated data revealed clear performance differences. For simulated data with induced missing dementia information, only the proposed bias-aware NPE achieved good performance (Figure 3C). Indeed, the bias-aware NPE recovered all transition hazards with accuracy comparable to the standard NPE trained on full data. Consistent with this result, both models achieved similar C2ST statistics on their respective validation sets when the training regime matched the datagenerating process (0.56 for observed data and 0.57 for full data). In contrast, applying the standard NPE to observed datasets substantially degraded hazard recovery. The largest errors occurred for the healthy-to-death h02 transition, but bias was also evident in the remaining transitions (Figure 3C). The dementia-to-death transition h12 is in general hard to recover, as indicated by a low contraction of the

θ Posterior

Trained with S

Classifier

Y Outcomes

4.0

×103

0.0

Trained without S

Epoch 1

Mean C2ST=0.61 p-value=0.50 4.0

Epoch 1

×103

Epoch 2

Mean C2ST=0.64 p-value=0.50 4.0

Epoch 2

×103

Epoch 3

Mean C2ST=0.65 p-value=0.50 5.0

Epoch 3

×103

Epoch 4

Mean C2ST=0.63 p-value=0.60

Epoch 4

×10 ×10 2.5×10 2.0 2.0 Mean C2ST=0.64 Mean C2ST=0.61 Mean C2ST=0.65 Mean C2ST=0.63 p-value=0.60 p-value=0.50 4.0 p-value=0.50 4.0 p-value=0.50 5.0 4.0 0.0 0.0 0.0 0.0

2.0

C 2.0

θ Prior

×103

1.0 0.5

×105

×104 ×104 2.5 ×104 2.0 2.0 Mean C2ST=0.92 Mean C2ST=0.96 Mean C2ST=0.96 Mean C2ST=0.98 p-value=0.00 p-value=0.000.0 p-value=0.000.01.0 p-value=0.00 0.0 1.0 2.0

0.5

0.0 0.0 0.0 0.0 0.0000 0.0005 0.0010 0.0000 0.0005 0.0010 0.0000 0.0005 0.0010 0.0000 0.0005 0.0010

a01

1.0

score C2ST score C2ST scoreC2STper bin) (mean per bin) (mean per (mean bin)

B Cl

Density of a01 Density of a01 Density of a01

0.8 1.0 0.6 0.8

1.0 0.6 0.8 0.6

a01

Figure 4: Validating bias correction on the Framingham Heart Study. (A) A classifier is trained on parameter-data pairs to predict if the pair belongs to the posterior or the joint. (B) Marginal posteriors of the bias-aware NPE for the baseline transition scale parameter a01 (healthy to dementia) for all 4 epochs. Posterior samples were tested using a classifier-based diagnostic, where a C2ST score of 0.5 means that samples cannot be distinguished from samples of the joint distribution. (C) Marginal posteriors of a NPE trained only on full data. Colors indicate classifier-based diagnostics, also trained only on data that was not censored. hazard scale parameter (Table S1) due to the low number of observed transitions. These results show that explicitly incorporating the observation process into the simulator is necessary for accurate parameter recovery under death-induced selection. The real data were analyzed using a further approach: a tailored spline-based full likelihood approach derived from the same multi-state formulation, which was previously used to analyze the data [Binder et al., 2019]. Applying all methods to the real Framingham data yielded estimates for the ageand sex-adjusted cumulative hazards for dementia, direct death, and death after dementia (Figure 3D) as well as the impact of covariates (Figure S2). For the age- and sex-adjusted cumulative death hazards, we observed good agreement between the spline-based likelihood approach and our proposed bias-aware NPE, but not with the NPE trained on full data, in line with our previous findings on simulations. In the first epoch, the cumulative dementia hazard estimated with the bias-aware NPE was slightly higher than with NPE trained on full data or the likelihood approach. For this transition, the NPE trained on full data showed more uncertainty than the other approaches. All approaches showed high uncertainty for cumulative dementia-to-death transition hazards (Figure 3D) due to the low number of known dementia-to-death transitions, reflecting limited identifiability of parameters governing death rates after dementia onset. For the impact of covariates, we found overall good agreement between methods with high uncertainty for the estimates related to sex of the NPE trained on full data (Figure S2). We additionally analyzed the data using a naive Cox proportional hazards model, which censors individuals who die without a recorded dementia diagnosis at their last dementia-free visit. This ignores the possibility of dementia onset between the last dementia-free observation and death, leading to systematic bias [Binder et al., 2019]. However, the resulting bias differs from that of the NPE trained on full data (Figure S2). The Cox model tends to overestimate cumulative death hazards while underestimating dementia hazards, especially in later epochs, as deaths are fully observed but dementia events are partially unobserved. In contrast, the NPE trained on full data underestimates death hazards compared to all other methods due to a mismatch between training and observed data. This shows that simulation-based methods can fail in unpredictable ways when simulated data and real data do not match. To provide an additional verification of the reliability of the estimates obtained using the bias-aware NPE, we assessed the calibration of the inferred baseline transition scale parameter from healthy to dementia, a01 , across study epochs using C2STs (Figure 4A). Posterior samples from the bias-aware NPE generated from the real data were not statistically distinguishable from samples drawn from the joint simulator, with C2ST scores around 0.65 across epochs (Figure 4B). In contrast, posterior samples from the NPE trained exclusively on full data were clearly distinguishable from the joint distribution when evaluated on the real censored data, with C2ST scores close to 1.0 and statistically significant discrepancies (Figure 4C). This mismatch indicates that ignoring censoring during training leads to an incorrect joint distribution and invalid posterior inference, while C2ST successfully detects the resulting bias even though neither the inference network nor the classifier was trained on censored data.

Transmission model

A Infection parameters θ Continuous-time household transmission model

Covariates C Age, vaccination

Inverse

Infections Y Missing not at random

Selection S School

1.0

Fraction of age group Omicron Alpha

Parameters

0.5 0.0 1.0 0.5 0.0

Infant

Random Selection Real Data

Child

transm pro acq pro Child sus Infant sus asym Adult inf asym Child inf asym Infant inf sym Child inf sym Infant inf

Adult

Child selection Adult selection

Alpha parameter value MCMC

Omicron parameter value Bias-aware NPE

Figure 5: Correcting for bias under multiple selection procedures on simulated data and the PedCovid Study. (A) Visualization of selection bias based on outcome and covariates. (B) Fraction of household members who were the first to test positive in each age group, shown for real data and for different simulated selection mechanisms (infections occurring on the same day are assigned to the younger age group). A total of 300 datasets with identical parameters were simulated for each mechanism. (C) ECDF of MCMC and NPE plots for the 3 different selection procedures. (D) MCMC and NPE results for real data with child selection procedure for Alpha (left) and Omicron (right) variants.

Enabling Unbiased Inference in Complex Stochastic Models and Study Designs with Intractable Likelihoods Having demonstrated reliable performance in applications with established estimators, we next assess whether the proposed framework can address problems for which a statistically coherent treatment remains challenging. We therefore consider inference under selection bias in complex stochastic simulation models where structured study inclusion distorts the observed data and the complexity of the underlying processes makes explicit likelihood-based correction infeasible (Figure 5A). As a representative example, we consider the PedCovid Study, a prospective longitudinal household study of SARS-CoV-2 transmission in which households were enrolled if there was a child that tested positive prior to the study inclusion date [Delaunay-Moisan et al., 2022]. This outcome- and covariate-dependent inclusion criterion induces strong selection bias when inferring transmission parameters. We evaluate the proposed framework using both simulated datasets and the observed PedCovid data. For the simulation study, household-level transmission data were generated using a mechanistic transmission model (see Methods). Observed households were replicated multiple times to construct synthetic populations, and infection dynamics were simulated conditional on parameters drawn from a prior distribution. To investigate the effect of study inclusion on parameter inference, datasets were generated under multiple selection schemes, including random household inclusion, child-dependent inclusion, and adult-dependent inclusion. Under random inclusion, households were eligible if the selected member tested positive before the inclusion date. Under child-dependent inclusion, the selected member additionally had to be younger than 18 years, whereas under adult-dependent inclusion, the selected member had to be 18 years or older. These inclusion criteria alter the composition of the observed households and thereby affect the distribution of infection characteristics, such as the age of the first infected household member (Figure 5B), a key determinant of transmission dynamics given known differences in immune responses across age groups [Manfroi et al., 2024]. Inference was performed using a bias-aware NPE trained on pairs of biased datasets and ground-truth 8

parameters, with an explicit indicator of the selection mechanism included as part of the training data. For comparison, we additionally performed likelihood-based inference using Markov chain Monte Carlo (MCMC) with the No-U-Turn Sampler implemented in Stan [Carpenter et al., 2017]. This likelihoodbased formulation assumes random household inclusion, which is computationally tractable but does not account for outcome-dependent selection. MCMC inference required over 12 hours of wall-clock time using 96 CPUs in parallel, whereas inference with the trained bias-aware NPE was performed in a matter of seconds on a single GPU. Evaluation on simulated data showed that under random household inclusion, both MCMC and the bias-aware NPE recovered posterior distributions consistent with the ground truth. Under outcomeor covariate-dependent inclusion, the bias-aware NPE consistently recovered the true parameter values (Figure S3A). However, MCMC produced systematically biased posterior estimates as reflected in the simulation-based calibration results (Figure 5C): the empirical cumulative distribution functions (ECDFs) for MCMC deviated markedly from the uniform distribution under biased sampling, while those obtained from the amortized approach remained well calibrated. By conditioning explicitly on the study inclusion mechanism, the amortized approach yields parameter estimates that remain coherent across virus variants and selection mechanisms. The trained bias-aware NPE was subsequently applied to the real PedCovid data. Across both the Alpha and Omicron variant cohorts, posterior estimates obtained from the bias-aware NPE and from the likelihood-based MCMC approach differed systematically, with the largest discrepancies observed in parameters governing age-related differences in infection risk (µinf ) and the household-size-dependent transmission parameter δ (Figure 5D). In the Alpha variant, infection-related modifiers, particularly those associated with asymptomatic individuals, were substantially larger under the bias-aware approach, whereas the MCMC estimates shrank these effects toward values near or below one (Figure 5D). Protection parameters also shifted qualitatively for both variants, with bias-aware inference generally implying weaker protection against transmission and smaller effects on acquisition relative to MCMC for the Alpha variant and stronger protection and acquisition effects for Omicron. These patterns suggest that ignoring the selection mechanism induces structural bias in several transmission parameters while leaving the baseline transmission rate largely unaffected. Simulation-based calibration further supports this interpretation. Under MCMC, symptomatic infectionrelated modifiers were systematically underestimated, whereas susceptibility-related modifiers were on average overestimated across simulated datasets for the child selection procedure (Figure 5C). On the Infant Child real data, µsym and µsym appear underestimated relative to the bias-aware NPE posterior inf inf (Figure 5D), which aligns with the SBC direction. Finally, classifier-based diagnostics confirm the consistency of the bias-aware posterior with the assumed generative model. For the real PedCovid data, the classifier accuracy was 0.719 (p = 0.50) for the Alpha variant and 0.766 (p = 0.20) for the Omicron variant, indicating no statistically significant discrepancy between posterior samples obtained with the bias-aware NPE and samples drawn from the joint simulator. Together, these results demonstrate that incorporating the selection mechanism directly into the simulator enables calibrated posterior inference even in complex stochastic models where explicit likelihood-based correction of selection bias is infeasible.

Discussion In this study, we developed a general amortized Bayesian inference framework that accounts for selection bias in parameter estimation tasks. By embedding selection, censoring, and missingness within forward simulations and training neural networks on the resulting joint distribution, the proposed approach yields calibrated and testable posterior inferences, even under complex selection mechanisms. Across three simulated and real epidemiological applications, the method matches analytic solutions where such solutions are available and remains applicable in regimes where the formulation of likelihoods accounting for the selection process is difficult. A key strength of the framework is that it requires only a simulator of the population process together with the selection mechanism. In many application areas, such simulation models are routinely developed to study bias through Monte Carlo experiments, sensitivity analyses, or study design investigations [Burton et al., 2006, Rothman et al., 2008, Kawabata et al., 2024]. Our approach leverages this established modeling practice for inference rather than using simulations solely for validation or sensitivity analysis. Because the method operates entirely through forward simulation, it can accommodate arbitrary

selection operators, including nonlinear, stochastic, and latent-state-dependent mechanisms, provided that they can be simulated. This flexibility allows the framework to be applied in settings where analytical likelihood corrections are unavailable or difficult to derive. More broadly, the growing role of simulation-based approaches in statistical workflows [Bürkner et al., 2025] highlights the potential of such methods to address complex data-generating processes that cannot be easily expressed in closed form. An additional advantage of the framework is amortization. Once trained, the inference network can be applied to multiple datasets without retraining, enabling rapid posterior estimation across simulated scenarios and real-world studies in contrast to neural likelihood-based approaches [Boyd et al., 2024]. This capability makes the approach particularly useful for systematically exploring the consequences of selection bias under a range of plausible assumptions. For example, once a simulator has been specified and an inference network trained, researchers can efficiently evaluate how alternative data-collection procedures or selection mechanisms affect inference without re-deriving likelihoods or repeating computationally expensive analyses. In this way, amortized simulation-based inference provides a practical tool not only for bias-aware parameter estimation but also for sensitivity analysis [Elsemüller et al., 2024] and experimental design [Bracher et al., 2025]. A potential concern is misspecification of the selection mechanism. While the proposed framework requires specifying a structural form for the selection process, such misspecification does not remain silent: discrepancies between the training data and the observed data can be detected through the C2ST diagnostics. Importantly, the selection mechanism need not be fully known. Unknown components, such as parameters governing prevalence or inclusion probabilities, can be incorporated into the generative model as additional latent variables and inferred jointly within the same amortized framework. In this way, uncertainty about the selection process itself can be propagated through the posterior, provided that a plausible structural form can be specified. This capability distinguishes the proposed approach from most existing methods addressing bias or missingness. Previous work has primarily focused on settings with partially observed covariates or missing data under missing-at-random assumptions [Lueckmann et al., 2017, Wang et al., 2024, Gloeckler et al., 2024, Verma et al., 2025, Simkus and Gutmann, 2025], where the observation process is either assumed known or modeled through auxiliary imputation mechanisms. Other approaches target specific bias structures, such as confounding addressed through instrumental-variable assumptions [Braun et al., 2025], or aim to improve robustness to model misspecification [Bharti et al., 2026, Khoo et al., 2026]. These methods typically rely on additional identification conditions or prespecified observation models and therefore do not allow uncertainty in the selection mechanism itself to be inferred jointly with the parameters of interest. By contrast, our framework treats the selection process as an explicit component of the generative model, allowing unknown aspects of the selection mechanism to be parameterized and marginalized over during inference. Consequently, the approach does not eliminate the need for substantive modeling assumptions; rather, it reframes selection bias from a problem of post hoc statistical correction to one of explicit generative modeling, where uncertainty about both the outcome and selection processes is propagated coherently through the posterior. As with other neural posterior estimation methods, performance depends on the choice of network architecture, the available simulation budget, and the fidelity of the underlying simulator. Insufficient coverage of the parameter space or simulator misspecification may lead to biased or poorly calibrated posteriors [Schmitt et al., 2023]. In addition, the classifier two-sample tests used for diagnostic evaluation operate on the low-dimensional embedding learned by the summary network rather than on the raw data. While such dimensionality reduction is necessary for computational tractability in high-dimensional settings, discrepancies that are not captured by the learned summary statistics may remain undetected. Consequently, both posterior quality and diagnostic reliability depend on the expressiveness of the neural networks and on how well the simulated training distribution reflects the relevant range of real-world data-generating processes. In this context, well-behaved simulation-based calibration results provide empirical evidence that the chosen network architectures are sufficiently expressive for the problem at hand, given the available simulation budget. In summary, we establish a simulation-based framework for Bayesian inference under selection bias that embeds the selection mechanism directly within the generative model. By combining amortized neural posterior estimation with explicit modeling of the data-collection process and formal diagnostic evaluation, the approach yields calibrated inference in settings where likelihood-based corrections are unavailable or impractical. We provide a reusable implementation directly embedded in state-of-the-art tools with well-documented examples. As simulation-based methods become increasingly central in

applications involving complex and structured data acquisition, this framework provides a principled basis for integrating selection mechanisms into routine probabilistic analysis.

Methods Problem Description We consider posterior inference p(θ | Y, S = 1; C) given the selected observations Y = {yi }N i=1 with S=1 and the corresponding covariates C = {ci }N . Selection is considered a missing data problem i=1 [Howe et al., 2015]. We write Yfull = (Y, Ymis ), where Y denotes the observed values and Ymis denotes the missing values (due to selection). Then the likelihood has the form ZZ p(Y, S | θ; C) = p(Y, Ymis | θ; C, Cmis ) p(S | Y, Ymis , θ; C, Cmis ) dCmis dYmis , (1) where p(S | Y, Ymis , θ; C, Cmis ) reflects the probability of a sample being selected during the selection process. In addition, the selection probability might depend on the latent variables, which can be computed from model parameters and covariates. To get the unbiased posterior distribution p(θ | Y, S; C) ∝ p(θ)p(Y, S | θ; C),

(2)

all unobserved variables must be marginalized. For a detailed introduction to selection bias, we refer to [Little and Rubin, 2019]. We differentiate between three scenarios according to the respective missingness type. First, if the selection probability p(S | Y, Ymis , θ; C, Cmis ) depends on neither Y or Ymis nor the covariates C or Cmis , then this is called missing completely at random. This is the simplest case, because missingness can be ignored. Second, if the selection process depends only on Y or C and not on Ymis and Cmis , the integrals in (1) simplify to p(Y, S | θ; C) = p(Y | θ; C)p(S | Y, θ; C). This scenario is called missing at random. Furthermore, if θ = (ϕ, ψ) and p(Y, S | θ; C) = p(Y | ϕ; C)p(S | Y, ψ; C), and the prior factorizes as p(θ) = p(ϕ) p(ψ), then the missing-data mechanism is ignorable for Bayesian inference [Little and Rubin, 2019]. However, inferred population-level quantities are only valid for the observed covariates. If the target of inference is a population-level parameter, such as the prevalence of infection in the full population rather than in the selected sample Y, then the selection depending on C can no longer be ignored. In this case, the observed covariate distribution p(C) differs from the population distribution p(Cfull ) and must be explicitly corrected. Third, when the selection mechanism depends on unobserved values such as Ymis and Cmis the data are not missing at random. In this case, the selection process cannot be ignored for inference, and all unobserved quantities involved in the selection mechanism must be explicitly modeled and marginalized to obtain unbiased posterior inference. Failure to account for the dependence of selection on unobserved variables generally leads to biased parameter estimates and invalid uncertainty quantification [Little and Rubin, 2019].

Amortized Bayesian Inference Under Selection Bias We train a neural posterior estimator (NPE) to learn a global mapping from datasets to posterior distributions across repeated simulations from a generative model. Let θ ∼ p(θ) denote parameters drawn from a prior distribution, and Y ∼ p(Y | θ; C) denote data generated from the forward model. We sample C from the data, resulting in an implicit prior. We simulate parameter-dataset pairs (θ, Y; C) from the joint distribution p(θ, Y, C) and train a NPE q(θ | Y; C) to approximate the true posterior p(θ | Y; C). Here, we instantiate the NPE using two neural networks: a summary network and a generative neural network as an inference network, for example, a flow matching or consistency model (see [Arruda et al., 2025] for a detailed introduction to generative neural networks and the corresponding objective functions). Both networks can be trained jointly, where the summary network provides a lower-dimensional summary of the simulation to the inference network. Under sufficient model capacity and training data, this procedure yields a consistent approximation of the posterior [Radev et al., 2020]. 11

In the presence of selection bias, we augment the generative model with an explicit selection mechanism. Parameter-dataset pairs are generated as before from the structural population model, yielding a latent population-level dataset. A selection mechanism is then applied to this dataset, resulting in a random variable S, representing sampling or missingness processes that determine which units are observed. In this way, selection acts on the fully simulated population data prior to observation. Because the framework relies only on forward simulation, it accommodates arbitrary and potentially complex selection mechanisms provided that they can be simulated. The resulting estimator approximates p(θ | Y, S = 1; C) without requiring the evaluation of a tractable likelihood or explicit estimation of selection probabilities. Moreover, we train the NPE across multiple selection mechanisms by parameterizing or indexing the selection process and conditioning the network on the corresponding selection indicator or mechanism-specific variables. The specific NPE architecture used in each experiment is described in the following sections.

Simulation-Based Diagnostics for Bias Detection To assess the calibration of the approximate posterior distributions, we employed simulation-based calibration (SBC) [Talts et al., 2020]. SBC evaluates whether posterior uncertainty is correctly quantified by repeatedly simulating parameters and datasets from the prior and generative models, performing inference, and computing the rank of each true parameter value within its corresponding posterior sample. If the posterior is calibrated, these ranks are uniformly distributed across repeated simulations. In practice, we generated M independent parameter draws θ (m) ∼ p(θ), simulated the corresponding (m)

datasets Y(m) , and obtained posterior samples θ̃ k from the trained NPE (or any other inference algorithm). For each parameter component, we computed the rank of θ (m) among the K posterior samples. Deviations from uniformity in the empirical ranks indicate systematic bias, overconfidence, or underdispersion in the posterior approximation. This deviation can be graphically assessed by computing the empirical cumulative distribution function (ECDF) and comparing it with the known uniform CDF [Talts et al., 2020]. In selection settings, SBC is performed under the full generative model, including the selection mechanism, allowing controlled quantification of bias across different selection regimes. As an additional diagnostic, we used a classifier two-sample test (C2ST) [Lopez-Paz and Oquab, 2017] to assess the fidelity of the learned posterior approximation. C2ST evaluates whether samples drawn from the learned posterior, combined with simulations, are statistically indistinguishable from samples drawn directly from the true joint distribution [Yao and Domke, 2023, Linhart et al., 2023]. Specifically, we constructed two sets of samples: (i) samples obtained by drawing a single θ ∼ q(θ | Y) for each dataset and pairing it with the corresponding Y (ii) joint samples (θ, Y) generated from the prior and forward models, including the selection mechanism. A binary classifier was then trained to distinguish between these two sample sets. If the posterior approximation is accurate, the classifier should achieve chance-level performance on the held-out data. In practice, classifiers were instantiated as a neural network and trained using 5-fold cross-validation on validation data not used for training the NPE. Classification accuracy substantially above chance indicates discrepancies between the learned posteriordata pair and the true joint distribution, suggesting insufficient training of the NPE or remaining bias. To assess the statistical significance for a given observed dataset, we additionally trained B=10 classifiers on label-permuted data, following established C2ST procedures [Linhart et al., 2023]. We define the test statistic as 2 K 1 X 1 T = ak − , (3) K 2 k=1

where ak denotes the classification accuracy for the parameter-data pair k. The empirical distribution of test statistics Tb under permutations b = 1, . . . , B provides a reference against which the observed test statistic Tobs is compared, yielding a p-value:   B 1 X p= I(Tb ≥ Tobs ) . (4) B b=1

Beyond validation on simulated data, C2ST can be applied to real datasets by applying the classifier to observed data paired with posterior samples. This enables the detection of bias or violations of 12

modeling assumptions even when the true parameters are unknown. To simplify the training of the classifiers, we first projected the data Y onto a lower-dimensional embedding using the summary network of the NPE. The classifiers then consisted of two hidden multi-layer perceptrons with widths such that the width was larger than 10 times the input dimension, as implemented in the BayesFlow library [Kühmichel et al., 2026].

Estimating Prevalence Under Selection Bias We consider estimation of the population-level prevalence ρ = E[Yfull ] under the selection p(S | Cfull ) motivated by the KoCo19 Study [Radon et al., 2020]. The simulation framework incorporates covariatedependent infection risk, non-representative sampling, outcome and covariate missingness, and imperfect diagnostic testing. Logistic Infection Model and Population Construction Individual-level infection status was modeled using logistic regression with odds ratio parameterization. The linear predictors included sex, age group, country of birth, and household size, with reference categories corresponding to male sex, age 20–34 years, birth in Germany, and single-person households. The model parameters consisted of an intercept on the log-odds scale and log-odds ratios for each non-reference category: ỹi ∼ Bernoulli(pi ), logit(pi ) = β0 + β ⊤ ci . (5) For simulation-based experiments, parameters were drawn independently from Gaussian prior distributions on β0 with mean −3 and standard deviation of 1.0 and on β with mean 0 and standard deviation of 0.5. The prior was chosen to be mostly uninformative and to obtain simulated prevalence below 50% of the population. The observed test results were then generated from ỹi using a misclassification model parameterized by externally estimated test sensitivity (Se=0.886) and specificity (Sp=0.997) [Olbrich et al., 2021], yielding apparent infection status yi . To construct latent population-level data, observed KoCo19 participants were oversampled using iterative proportional fitting, that is, p(S | C) was estimated to match known marginal population distributions of Munich from [Radon et al., 2020] with respect to age, sex, household size, and country of birth. For simulation experiments, we generated a synthetic population equal to 10% of the Munich population; for the real data analysis, we generated the full 1.5 million inhabitants. Missing covariates were imputed during this step by sampling from the corresponding population-level target distributions, as they were assumed to be missing completely at random. The resulting oversampled dataset approximated a synthetic population consistent with known census margins, where true infection outcomes were then simulated using the above logistic model and prevalence calculated as ρ = E[Yfull ]. Outcome missingness was imposed by reproducing the missingness pattern observed in the original KoCo19 data. Hence, individuals with missing outcomes in the original study also had missing outcomes after oversampling, thereby preserving the empirical missing at random structure. To emulate the original study design, the synthetic population was subsequently downsampled to the original cohort size using inverse-probability subsampling based on oversampling weights, yielding a biased sample that matched the observed KoCo19 cohort. This procedure preserved both the selection bias and missingness structure of the real data. Prevalence Estimation Two baseline prevalence estimators were considered. First, an unadjusted estimator based on complete cases was computed from the apparent test outcomes of the (simulated or real) KoCo19 cohort and corrected for test misclassification using the Rogan-Gladen estimator [Rogan and Gladen, 1978]: ρ̂RG =

ρ̂obs + Sp − 1 . Se + Sp − 1

(6)

Second, inverse probability weighting was applied to account for unequal sampling probabilities induced by the subsampling step, followed by misclassification correction. For bootstrap-based uncertainty quantification [Efron and Tibshirani, 1994], the entire pipeline was repeated with 100 bootstrap resamples 13

of the original KoCo19 data. In the analyses of the real data, simulated outcomes were replaced by observed test results, while the same inference procedures were applied unchanged. Neural Posterior Estimation To perform amortized Bayesian inference, we used a NPE with a simulation-based workflow tailored to the KoCo19 cohort. We used the simulation pipeline described above to generate training data consisting of J = 10,000 pairs {(Yj , ρj )}Jj=1 of subsampled cohorts Yj with a known prevalence ρj and 1000 additional validation pairs. Here, we aim to directly estimate the population-level prevalence and not the underlying parameters of the logistic model, which we marginalize implicitly. The observation Yj includes the observed test results yj,i and covariates cj,i for all individuals i ∈ {1, . . . , 5577}. As all covariates ci were categorical, we converted them to non-negative integers. Missing values were encoded using a fixed value of −1, following the approach of [Wang et al., 2024], allowing the network to learn to distinguish between observed and unobserved inputs. The workflow then proceeded by training a NPE, consisting of an inference network and a summary network, which were trained jointly using the training data. The input to the summary neural network consisted of Yj and the index of the epoch ej ∈ {1, 2, 3, 4, 5} used as the basis for oversampling, which determined the distribution of the covariates and the amount of missingness. The summary network was instantiated as a deep set encoder that can learn permutation-invariant representations of set-based data [Zaheer et al., 2017], with a summary dimension of 4. This choice ensured that the summary network was agnostic to the order of individuals in the cohort. For the inference network, we used a small conditional flow matching model [Wildberger et al., 2023], as it also allows for posterior density estimation, with 2 conditional multi-layer perceptrons of width 128 as the backbone and a dropout rate of 0.1. For an in-depth introduction to SBI with generative models, refer to [Arruda et al., 2025]. For training, we used the inference pipeline implemented in the BayesFlow 2.0.8 library [Kühmichel et al., 2026] with jax as the backend, utilizing the default training settings. We trained for 300 epochs using the AdamW optimizer [Loshchilov and Hutter, 2019], with a batch size of 64, an initial learning rate of 5×10−4 , and a cosine schedule for the learning rate [Loshchilov and Hutter, 2017]. Furthermore, the prevalence was constrained between 0 and 1 by applying a sigmoid transformation before training. To evaluate the different inference procedures, we computed the absolute difference between the estimated prevalence and the true prevalence using 1000 newly generated simulations. With the NPE, inference is amortized; hence, after training the model once, we sampled 500 posterior samples and additionally evaluated the posterior density for each sample to approximate the mode of the posterior. We then used the mode as a point estimator of the prevalence for comparison with other baseline estimators. Furthermore, we assessed the reliability of the approximate posteriors produced by the NPE using C2ST.

Estimating Time-to-Event Under Missing Disease due to Death Bias We consider time-to-event analysis in a semi-competing risks setting in which death may preclude the observation of dementia onset, as in incidence studies based on the Framingham Heart Study [Satizabal et al., 2016, Binder and Schumacher, 2016]. Let Y denote the time-to-dementia or death outcome and S the selection indicator induced by unobserved dementia onset before death. The target is the posterior p(θ | Y, S; C) where censoring induces the selection probability p(S | Y, Ymis , θ, C, Cmis ), i.e., not missing at random. This structure is embedded in a SBI framework. Illness-Death Multi-State Model Event dynamics were modeled using an illness-death multi-state model with three states: healthy (0), dementia (1), and death (2). Transitions between states followed cause-specific hazard models for the transitions (0 → 1) (dementia onset), (0 → 2) (death without prior dementia), and (1 → 2) (death after dementia), following [Joly et al., 2002, Binder et al., 2017, Binder et al., 2019]. Baseline hazards (h01 , h02 , h12 ) were parameterized by transition-specific scale parameters akl (k → l) and shape parameters κkl governing Weibull hazard functions. The covariate effects β kl of sex and age were included as transition-specific log-linear modifiers of the transition hazards: hkl (t | c) = akl κkl tκkl −1 exp(β ⊤ kl c). 14

(7)

Age was centered before inclusion in the linear predictors. For simulation-based experiments, the transition parameters were drawn independently from the prior distributions. Baseline scales akl and shape parameters κkl were assigned Gamma priors, parameterized through their mean m and coefficient of variation v. Specifically, akl ∼ Gamma(ma = 0.0002993, va = 1) and κkl ∼ Gamma(mκ = 1, vκ = 0.25), corresponding to moderate deviations from exponential hazards. Covariate effects for sex and age were assigned independent Gaussian priors N (0, 12 ). This prior specification induced realistic heterogeneity in event times while remaining weakly informative (Figure S2). Given the sampled parameters and individual-level covariates, complete event trajectories were generated by inverse transform sampling from cause-specific hazard functions. For each individual, competing event times for dementia onset and direct death were simulated, with the earliest event determining the first transition. An additional transition time from dementia to death was simulated for individuals who experienced dementia onset. Visit-Based Censoring and Observation Process The Framingham data were obtained from the Framingham Heart Study, a long-term cohort study initiated in 1948. It comprises an original cohort of 5209 residents of Framingham, Massachusetts, with repeated examinations every two years, and an offspring cohort of 5214 participants initiated in 1971 with examinations approximately every four years. Following [Satizabal et al., 2016] and [Binder et al., 2019], both cohorts were combined for analysis, resulting in comparable sample sizes across four non-overlapping epochs, each spanning a five-year period. To reflect the structure of longitudinal cohort studies, simulated event times were transformed into observed data through a visit-based censoring mechanism extending [Binder et al., 2019]. Each individual was assigned two observation times corresponding to the study visits. Dementia status was only observable at visit times, inducing censoring of dementia onset. If death occurred before dementia was detected at a visit, the dementia status was censored to the last known status. If dementia occurred after the last visit, the dementia status was censored. Administrative censoring was applied at the end of the follow-up. If death or dementia occurred after the study ended, the corresponding event was censored. If the patients were disease-free at the end of the epoch and continued in the next epoch, the last visit time point was set to the end of the epoch. In addition, to reflect random dropout, healthy individuals were assumed to have a 5% probability of leaving the study, independent of covariates and event times. The mean visit times were estimated empirically for each study epoch. For each individual, two visit times were then generated by sampling around these epoch-specific means: the first visit time, corresponding to observations within the first 2.5 years of follow-up, was drawn with variance 100, and the second visit time, corresponding to observations after 2.5 years, was drawn with variance 150. Simulations were conducted separately for each of the four Framingham Heart Study epochs. For each epoch, individual-level covariates were obtained directly from empirical data. For each simulated dataset, the model parameters were sampled from the priors. This procedure resulted in simulated data that closely mimicked the visit structure of the Framingham Heart Study. Neural Posterior Estimation Simulated visit-censored datasets were used to train a neural posterior estimator within a SBI framework. We generated 20,000 training datasets consisting of pairs of censored cohorts Yj and the 12 corresponding parameters θ = (akl , κkl , β kl ) of the multi-state model. Each individual yj,i ∈ Yj was then described as a vector of the time of illness, binary indicator of illness, time of death, binary indicator of death (all possibly censored), sex, age, and epoch ID. We standardized the covariates and time variables before training and log-transformed the scale and shape parameters. Yj was padded with −1, such that it had the same number of individuals (2299) for every epoch. The summary network of the NPE was instantiated as a set transformer to learn permutation-invariant representations [Lee et al., 2019] with a summary dimension twice the number of parameters. As an inference network, we used a consistency model due to its inference speed [Schmitt et al., 2024]. The backbone of the consistency model was chosen as the default 5 conditional multi-layer perceptrons of width 256 with a dropout rate of 0.1. We trained the model using the same settings as before in the BayesFlow library [Kühmichel et al., 2026]. Similar to the previous model, we assessed the reliability of the approximate posteriors produced by the NPE using the C2ST-diagnostic on a validation set with 1000 cohorts. 15

We compared our approach to a full likelihood approach and a naive Cox model [Joly et al., 2002, Binder et al., 2019]. For the Cox approach, separate models were fitted for each transition (0 → 1, 0 → 2, 1 → 2), with competing events treated as independent censoring. In contrast, the likelihood-based approach relies on an illness-death model, where all observed cases of transition, including missing disease by death, are handled through exact likelihood contributions derived from the multi-state process, while penalization enforces smoothness of the hazard functions [Joly et al., 2002]. Baseline hazards were parameterized by splines, and age- and sex-adjusted cumulative hazards were obtained by evaluating the linear predictor at the mean age and sex values, multiplying the baseline hazards by the resulting proportional hazards factor, and integrating the adjusted hazards over time following [Binder et al., 2019]. In addition, a NPE with identical architecture and training configuration was trained on uncensored data subject only to administrative censoring. Posterior calibration was assessed using C2ST on the corresponding validation set. Hence, the second NPE and corresponding classifier were trained exclusively on uncensored data.

Estimating Transmission Rates Under Selection We consider a stochastic household transmission model (described in detail in forthcoming work) with outcome-dependent study inclusion motivated by the PedCovid Study [Delaunay-Moisan et al., 2022]. Let Y denote the observed household-level infection data and S the study inclusion indicator. Inference targets the posterior p(θ | Y, S; C) under a selection mechanism p(S | Y, Ymis , θ, C, Cmis ), that depends on unobserved infection outcomes and thus corresponds to a not missing at random setting. Model Structure and Transmission Hazard Two SARS-CoV-2 variants (Alpha and Omicron) were separately modeled using variant-specific parameterizations. Let τi denote the infection time of individual i ∈ {1, . . . , N }, with τi = ∞ if individual i is never infected. At time t, a susceptible individual i experiences an infection hazard X λi (t) = α + β w(n) κ(t − τj ) µinf µsus µpro (8) j: τj <t

from other household members j, where α is a background hazard, β is the baseline transmission intensity, and w(n) = (n/4)−δ scales the transmission by household size n with exponent δ. The generation-time kernel κ(·) was modeled as a Gamma density, with parameters chosen separately for the Alpha variant (shape=2, rate=0.44) [Chen et al., 2022] and Omicron variant (shape=3.351, rate=1.1098) [an der Heiden and Buchholz, 2022]. The multiplicative terms µinf , µsus , and µpro encode heterogeneity in infectiousness, susceptibility, and protection by vaccination or past exposures. Specifically, µinf denotes the infectivity multiplier of an infectious individual j at time t, defined relative to symptomatic adults and varying by age category (infant <6 years, child 6–11 years, and adult >11 years) and symptom status (symptomatic or asymptomatic). The term µsus represents the susceptibility multiplier of a susceptible individual i, defined relative to adults and depending on age, while µpro captures the protection effects on transmission or acquisition. Priors were specified as weakly informative while reflecting plausible transmission dynamics. The baseline transmission intensity β was assigned a Gamma prior with shape 2.0 and scale 0.5, and the household size scaling exponent δ was assigned a standard normal prior. All multiplicative effects capturing infectiousness, susceptibility, and protection were assigned log-normal priors with a log-mean of 0 and log-standard deviation of 0.7. The background hazard α was fixed to 0.001 for the Alpha variant and 0.01 for the Omicron variant, due to weak identifiability of this parameter under random household selection (Figure S3). Over a short interval [t, t + ∆t), the probability that an individual i becomes infected is Pr i infected in [t, t + ∆t) = 1 − exp −λi (t)∆t .

(9)

Following infection, individuals were assigned as asymptomatic, with a probability of 0.4 for the Alpha variant and 0.3 for the Omicron variant. For symptomatic individuals, the incubation period was modeled by a Gamma distribution, with parameters differing by variant (Alpha: mean=4.42 d, SD=2.30 d; 16

Table 1: Probability of a household member missing a test triggered by a symptomatic member testing positive estimated from the PedCovid data. Household size Variant Alpha Omicron

0.00 0.00

0.10 0.40

0.07 0.24

0.14 0.17

0.10 0.17

0.71 0.14

– 0.13

Omicron: mean=3.09 d, SD=1.64 d) [Galmiche et al., 2023]. The symptom onset times are given by Di = τi + incubationi . Symptomatic individuals underwent testing after symptom onset, with testing delays drawn from a geometric distribution with a success probability of 0.33 for both variants. A test was assumed to yield a positive result between 1 and 15 days after infection. A positive test by a symptomatic individual triggered the testing of other household members with household-size-specific probabilities of missing tests, estimated empirically from the PedCovid data (Table 1). Triggered tests occurred after delays were drawn from geometric distributions (Alpha: p=0.48, Omicron: p=0.46). In addition, all individuals were subjected to background testing for other reasons, with fixed daily testing probabilities of 1/21 for Alpha and 1/14 for Omicron. For asymptomatic individuals, the observation time was defined as the time of the first positive test. The transmission dynamics were simulated using a fixed-step tau-leaping scheme with daily time steps. Household covariates were obtained directly from the PedCovid dataset, and each household was simulated 50 times. At each time step, susceptible and infectious individuals were identified, infection hazards were computed, new infections were sampled, and symptom status and incubation periods were assigned as appropriate. Testing events, inclusion procedures, and follow-up testing were then applied, after which households were selected according to household selection schemes. Household Selection Schemes Household inclusion in the study occurred once a household member tested positive. From all individuals with a positive test result within the age range defined by the selection scheme, a single inclusion case was randomly chosen. The household inclusion date was set to Dk + delay, where the delay followed a Poisson distribution with a mean of 4.8 days for Alpha and 2.7 days for Omicron. Follow-up tests in the household for all members were scheduled at 3, 7, 15, and 45 days after inclusion. All dates were shifted by 30 days to match the earliest infected individuals in the data. All geometric and Poisson delay distributions were fitted by maximum likelihood to the PedCovid data in the original study. From the full set of simulated households, subsets were selected using different selection mechanisms. The target sample sizes were N =128 households for the Alpha variant cohort and N =54 for the Omicron variant cohort. Under random selection, households were sampled uniformly without replacement. Under outcome-dependent selection, households were included only if the designated inclusion case met a specified age criterion (child: <18 years; adult: ≥18 years). The child-based selection scheme replicates the procedure used in the PedCovid Study [Delaunay-Moisan et al., 2022]. Likelihood for MCMC and Random Sampling Given the observed infection and testing histories, the likelihood of the household transmission model can be defined by combining the infection hazard, survival contributions for uninfected individuals, and incubation period densities for symptomatic cases under the assumption of random selection of households. The cumulative hazard for an individual i up to time t was X Λi (t) = α(t − τfirst,i ) + β w(n) Flag (t − τj ) µinf µsus µpro , (10) j: τj <t

where τfirst,i denotes the first infection time in the household of the individual i, and Z x Flag (x) = κ(u) du. 0

(11)

Unobserved infection times were treated as latent variables and constrained by observed testing data with soft quadratic penalties enforcing plausible temporal bounds. Using this likelihood, inference is possible with MCMC under the assumption of random selection. Specifically, we used the No-U-Turn Sampler as implemented in Stan [Carpenter et al., 2017] with 4 chains and 5000 samples after the burn-in. Neural Posterior Estimation Household datasets were used to train a neural posterior estimator within a SBI framework. We generated 32,000 simulated household datasets {(Yj , θ j )}Jj=1 , where each Yj consisted of all observed testing and infection information for the households selected under a given selection mechanism, and θ j = (µinf , µsus , µpro , δ, β) denotes the 11 corresponding transmission model parameters. Each individual in a household was represented by the features of age group, infection status (one-hot encoded), protection status, household size, virus variant, and selection mechanism (one-hot encoded), as well as the date of symptom onset (or the date of the first positive test for asymptomatic individuals), the date of the last negative test, the date of the first or last positive test, and the end of follow-up. To enable batched training, households were padded to a maximum size of eight members using fixed padding values of −1. Posterior inference was performed using a hybrid summary network similar to [Arruda et al., 2026] combining 2 set transformers [Lee et al., 2019]; one to capture within-household temporal structure and generate a latent representation of a household, and another one applied then on the set of representations to keep the permutation invariance property of households. The resulting summaries of size 3 times the number of parameters were passed to a flow matching model [Wildberger et al., 2023] with the default 5 conditional multi-layer perceptrons of width 256 with a dropout rate of 0.1 as backbone. The NPE was trained for 150 epochs on the simulated datasets with a batch size of 64 using the BayesFlow library [Kühmichel et al., 2026]. Inference performance was evaluated using 320 simulated datasets for each of the 3 selection procedures. Posterior inference obtained via the neural posterior estimator was compared with MCMC, which relies on an explicit evaluation of the biased model likelihood. Calibration of posterior distributions was assessed using SBC and C2ST.

Code Availability Code is available at github.com/arrjon/AmortizedSelectionBias.

Acknowledgments J.H. acknowledges financial support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy (EXC 2047 - 390685813, EXC 2151 - 390873048), by the European Union via ERC grant INTEGRATE (grant no 101126146), and by the University of Bonn (via the Schlegel Professorship of J.H.). N.B. acknowledges financial support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Project-ID 499552394 – SFB 1597. S.C. acknowledges financial support from Université de Versailles Saint-Quentin-en-Yvelines, the Inception program (Investissement d’Avenir grant ANR-16-CONV-0005), Institut Pasteur, and the HOME project (ANR 20-CE35-0016). The PedCovid study was funded by a grant from the French Ministry of Health (PHRC) and a grant from the ANR (RA-COVID-19). The authors thank the PedCovid working group, including Sylvie Behillil, Naı̈m Bouazza, Nelly Briand, Agnès Delaunay-Moisan, Flora Donati, Vincent Enouf, Jérémie Guedj, Marianne Leruez-Ville, Lulla Opatowski, Faheemah Padavia, Isabelle SermetGaudelus, Chloé Sturmach, Sylvie van der Werf, and the technical team of the National Reference Center for Respiratory viruses. We acknowledge the Marvin and Unicorn clusters hosted by the University of Bonn. We thank Simon Cauchemez for helpful discussions on bias in household studies.

Author CRediT J.A.: Conceptualization, Formal analysis, Methodology, Software, Validation, Visualization, Writing – original draft, Writing – review & editing. S.C.: Data curation, Methodology, Software, Writing – review & editing. P.S.: Data curation, Formal analysis, Software, Writing – review & editing. A.W.: Investigation, Writing – review & editing. M.H.: Investigation, Writing – review & editing. I.S.: Investigation, Writing – review & editing. N.B.: Data curation, Methodology, Resources, Software, Supervision, Writing – review & editing. L.O.: Conceptualization, Methodology, Supervision, Writing – review & editing. J.H.: Conceptualization, Methodology, Funding acquisition, Project administration, Resources, Supervision, Writing – review & editing.

References [an der Heiden and Buchholz, 2022] an der Heiden, M. and Buchholz, U. (2022). Serial interval in households infected with SARS-CoV-2 variant B.1.1.529 (Omicron) is even shorter compared to Delta. Epidemiology and Infection, 150:e165. [Arruda et al., 2026] Arruda, J., Alamoudi, E., Mueller, R., Vaisband, M., Molkenbur, R., Merrin, J., Kiermaier, E., and Hasenauer, J. (2026). Simulation-based inference of cell migration dynamics in complex spatial environments. npj Systems Biology and Applications. [Arruda et al., 2025] Arruda, J., Bracher, N., Köthe, U., Hasenauer, J., and Radev, S. T. (2025). Diffusion models in simulation-based inference: A tutorial review. arXiv preprint arXiv:2512.20685. [Austin, 2011] Austin, P. C. (2011). An introduction to propensity score methods for reducing the effects of confounding in observational studies. Multivariate Behavioral Research, 46(3):399–424. [Bharti et al., 2026] Bharti, A., Dellaporta, C., Hikida, Y., and Briol, F.-X. (2026). Amortised and provably-robust simulation-based inference. arXiv preprint arXiv:2602.11325. [Binder et al., 2019] Binder, N., Balmford, J., and Schumacher, M. (2019). A multi-state model based reanalysis of the framingham heart study. European Journal of Epidemiology, 34(11):1075–1083. [Binder et al., 2017] Binder, N., Herrnböck, A.-S., and Schumacher, M. (2017). Estimating hazard ratios in cohort data with missing disease information due to death: Estimation in data with missing disease information due to death. Biometrical Journal, 59(2):251–269. [Binder and Schumacher, 2014] Binder, N. and Schumacher, M. (2014). Missing information caused by death leads to bias in relative risk estimates. Journal of Clinical Epidemiology, 67(10):1111–1120. [Binder and Schumacher, 2016] Binder, N. and Schumacher, M. (2016). Letter to ”Incidence of Dementia over Three Decades in the Framingham Heart Study”. The New England Journal of Medicine, 375(1):92–93. [Boyd et al., 2024] Boyd, B. M., Grayling, M., Thorp, S., and Mandel, K. S. (2024). Accounting for selection effects in supernova cosmology with simulation-based inference and hierarchical Bayesian modelling. arXiv preprint arXiv:2407.15923. [Bracher et al., 2025] Bracher, N., Kühmichel, L., Ivanova, D. R., Intes, X., Bürkner, P.-C., and Radev, S. T. (2025). JADAI: Jointly amortizing adaptive design and Bayesian inference. arXiv preprint arXiv:2512.22999. [Braun et al., 2025] Braun, M., Peña, J. M., and Daoud, A. (2025). Flow iv: Counterfactual inference in nonseparable outcome models using instrumental variables. arXiv preprint arXiv:2508.01321. [Bürkner et al., 2025] Bürkner, P.-C., Schmitt, M., and Radev, S. T. (2025). Simulations in statistical workflows. arXiv preprint arXiv:2503.24011. [Burton et al., 2006] Burton, A., Altman, D. G., Royston, P., and Holder, R. L. (2006). The design of simulation studies in medical statistics. Statistics in Medicine, 25(24):4279–4292.

[Carpenter et al., 2017] Carpenter, B., Gelman, A., Hoffman, M. D., Lee, D., Goodrich, B., Betancourt, M., Brubaker, M., Guo, J., Li, P., and Riddell, A. (2017). Stan: A probabilistic programming language. Journal of statistical software, 76:1–32. [Chen et al., 2022] Chen, J., Qiu, Y., Shi, Y., Wu, W., Zheng, E., Xu, L., and Jia, M. (2022). Uncovering the impact of control strategies on the transmission pattern of SARS-CoV-2 - Ruili City, Yunnan Province, China, February–March 2022. China CDC Weekly, 4:1032. [Cochran, 1977] Cochran, W. G. (1977). Sampling techniques. Johan Wiley & Sons Inc. [Copas and Li, 1997] Copas, J. B. and Li, H. (1997). Inference for non-random samples. Journal of the Royal Statistical Society Series B: Statistical Methodology, 59(1):55–95. [Delaunay-Moisan et al., 2022] Delaunay-Moisan, A., Guilleminot, T., Semeraro, M., Briand, N., BaderMeunier, B., Berthaud, R., Morelle, G., Quartier, P., Galeotti, C., Basmaci, R., Benoist, G., Gajdos, V., Lorrot, M., Rifai, M., Crespin, M., M’Sakni, Z., Padavia, F., Savetier-Leroy, C., Lorenzi, M., Maurin, C., Behillil, S., de Pontual, L., Elenga, N., Bouazza, N., Moltrecht, B., van der Werf, S., Leruez-Ville, M., and Sermet-Gaudelus, I. (2022). Saliva for molecular detection of SARS-CoV-2 in pre-school and school-age children. Environmental Microbiology, 24(10):4725–4737. [Efron and Tibshirani, 1994] Efron, B. and Tibshirani, R. J. (1994). An introduction to the bootstrap. Chapman and Hall/CRC. [Elsemüller et al., 2024] Elsemüller, L., Olischläger, H., Schmitt, M., Bürkner, P.-C., Koethe, U., and Radev, S. T. (2024). Sensitivity-aware amortized Bayesian inference. Transactions on Machine Learning Research. [Elsemüller et al., 2025] Elsemüller, L., Pratz, V., von Krause, M., Voss, A., Bürkner, P.-C., and Radev, S. T. (2025). Does unsupervised domain adaptation improve the robustness of amortized Bayesian inference? a systematic evaluation. Transactions on Machine Learning Research. [Forrest et al., 2025] Forrest, I. S., Huang, K.-L., Eggington, J. M., Chung, W. K., Jordan, D. M., and Do, R. (2025). Using large-scale population-based data to improve disease risk assessment of clinical variants. Nature Genetics, 57(7):1588–1597. [Galmiche et al., 2023] Galmiche, S., Cortier, T., Charmet, T., Schaeffer, L., Chény, O., von Platen, C., Lévy, A., Martin, S., Omar, F., David, C., Mailles, A., Carrat, F., Cauchemez, S., and Fontanet, A. (2023). SARS-CoV-2 incubation period across variants of concern, individual factors, and circumstances of infection in France: a case series analysis from the ComCor study. The Lancet Microbe, 4(6):e409–e417. [Gloeckler et al., 2024] Gloeckler, M., Deistler, M., Weilbach, C., Wood, F., and Macke, J. H. (2024). All-in-one simulation-based inference. In Proceedings of the 41st International Conference on Machine Learning, pages 15735–15766. [Heckman, 1979] Heckman, J. J. (1979). Sample selection bias as a specification error. Econometrica: Journal of the Econometric Society, pages 153–161. [Horvitz and Thompson, 1952] Horvitz, D. G. and Thompson, D. J. (1952). A generalization of sampling without replacement from a finite universe. Journal of the American statistical Association, 47(260):663– 685. [Howe et al., 2015] Howe, C. J., Cain, L. E., and Hogan, J. W. (2015). Are all biases missing data problems? Current Epidemiology Reports, 2(3):162–171. [Joly et al., 2002] Joly, P., Commenges, D., Helmer, C., and Letenneur, L. (2002). A penalized likelihood approach for an illness–death model with interval-censored data: application to age-specific incidence of dementia. Biostatistics, 3(3):433–443. [Kaplan et al., 2014] Kaplan, R. M., Chambers, D. A., and Glasgow, R. E. (2014). Big data and large sample size: a cautionary note on the potential for bias. Clinical and Translational Science, 7(4):342–346. [Kawabata et al., 2024] Kawabata, E., Major-Smith, D., Clayton, G. L., Shapland, C. Y., Morris, T. P., Carter, A. R., Fernández-Sanlés, A., Borges, M. C., Tilling, K., Griffith, G. J., et al. (2024). Accounting for bias due to outcome data missing not at random: comparison and illustration of two approaches to probabilistic bias analysis: a simulation study. BMC Medical Research Methodology, 24(1):278. 20

[Khoo et al., 2026] Khoo, S., Prangle, D., Liu, S., and Beaumont, M. (2026). Minimum distance summaries for robust neural posterior estimation. arXiv preprint arXiv:2602.09161. [Kleinbaum et al., 1981] Kleinbaum, D. G., Morgenstern, H., and Kupper, L. L. (1981). Selection bias in epidemiologic studies. American journal of epidemiology, 113(4):452–463. [Kühmichel et al., 2026] Kühmichel, L., Huang, J. M., Pratz, V., Arruda, J., Olischläger, H., Habermann, D., Kucharsky, S., Elsemüller, L., Mishra, A., Bracher, N., Jedhoff, S., Schmitt, M., Bürkner, P.-C., and Radev, S. T. (2026). BayesFlow 2.0: Multi-backend amortized Bayesian inference in Python. arXiv preprint arXiv:2602.07098. [Lee et al., 2019] Lee, J., Lee, Y., Kim, J., Kosiorek, A., Choi, S., and Teh, Y. W. (2019). Set transformer: A framework for attention-based permutation-invariant neural networks. In International conference on machine learning, pages 3744–3753. PMLR. [Leffondré et al., 2013] Leffondré, K., Touraine, C., Helmer, C., and Joly, P. (2013). Interval-censored time-to-event and competing risk with death: is the illness-death model more accurate than the Cox model? International Journal of Epidemiology, 42(4):1177–1186. [Linero and Daniels, 2018] Linero, A. R. and Daniels, M. J. (2018). Bayesian approaches for missing not at random outcome data: the role of identifying restrictions. Statistical Science: a Review Journal of the Institute of Mathematical Statistics, 33(2):198. [Linhart et al., 2023] Linhart, J., Gramfort, A., and Rodrigues, P. (2023). L-C2ST: Local diagnostics for posterior approximations in simulation-based inference. Advances in Neural Information Processing Systems, 36:56384–56410. [Little and Rubin, 2019] Little, R. J. and Rubin, D. B. (2019). Statistical analysis with missing data. John Wiley & Sons. [Lopez-Paz and Oquab, 2017] Lopez-Paz, D. and Oquab, M. (2017). Revisiting classifier two-sample tests. In International Conference on Learning Representations. [Loshchilov and Hutter, 2017] Loshchilov, I. and Hutter, F. (2017). SGDR: Stochastic gradient descent with warm restarts. In International Conference on Learning Representations. [Loshchilov and Hutter, 2019] Loshchilov, I. and Hutter, F. (2019). Decoupled weight decay regularization. In International Conference on Learning Representations. [Lueckmann et al., 2021] Lueckmann, J.-M., Boelts, J., Greenberg, D., Goncalves, P., and Macke, J. (2021). Benchmarking simulation-based inference. In International Conference on Artificial Intelligence and Statistics, pages 343–351. PMLR. [Lueckmann et al., 2017] Lueckmann, J.-M., Goncalves, P. J., Bassetto, G., Öcal, K., Nonnenmacher, M., and Macke, J. H. (2017). Flexible statistical inference for mechanistic models of neural dynamics. Advances in Neural Information Processing Systems, 30. [Manfroi et al., 2024] Manfroi, B., Cuc, B. T., Sokal, A., Vandenberghe, A., Temmam, S., Attia, M., El Behi, M., Camaglia, F., Nguyen, N. T., Pohar, J., et al. (2024). Preschool-age children maintain a distinct memory CD4+ T cell and memory B cell response after SARS-CoV-2 infection. Science Translational Medicine, 16(765):eadl1997. [Olbrich et al., 2021] Olbrich, L., Castelletti, N., Schaelte, Y., Gari, M., Puetz, P., Bakuli, A., Pritsch, M., Kroidl, I., Saathoff, E., Guggenbuehl Noller, J. M., et al. (2021). Head-to-head evaluation of seven different seroassays including direct viral neutralisation in a representative cohort for SARS-CoV-2. Journal of General Virology, 102(10):001653. [Petersen et al., 2021] Petersen, J. M., Ranker, L. R., Barnard-Mayers, R., MacLehose, R. F., and Fox, M. P. (2021). A systematic review of quantitative bias analysis applied to epidemiological research. International Journal of Epidemiology, 50(5):1708–1730. [Radev et al., 2020] Radev, S. T., Mertens, U. K., Voss, A., Ardizzone, L., and Köthe, U. (2020). Bayesflow: Learning complex stochastic models with invertible neural networks. IEEE transactions on neural networks and learning systems. 21

[Radon et al., 2020] Radon, K., Saathoff, E., Pritsch, M., Guggenbühl Noller, J. M., Kroidl, I., Olbrich, L., Thiel, V., Diefenbach, M., Riess, F., Forster, F., et al. (2020). Protocol of a population-based prospective COVID-19 cohort study Munich, Germany (KoCo19). BMC Public Health, 20(1):1036. [Rogan and Gladen, 1978] Rogan, W. J. and Gladen, B. (1978). Estimating prevalence from the results of a screening test. American Journal of Epidemiology, 107(1):71–76. [Rothman et al., 2008] Rothman, K. J., Greenland, S., Lash, T. L., et al. (2008). Modern epidemiology, volume 3. Wolters Kluwer Health/Lippincott Williams & Wilkins Philadelphia. [Rudolph et al., 2023] Rudolph, J. E., Zhong, Y., Duggal, P., Mehta, S. H., and Lau, B. (2023). Defining representativeness of study samples in medical and population health research. BMJ Medicine, 2(1):e000399. [Satizabal et al., 2016] Satizabal, C. L., Beiser, A. S., Chouraki, V., Chêne, G., Dufouil, C., and Seshadri, S. (2016). Incidence of dementia over three decades in the framingham heart study. New England Journal of Medicine, 374(6):523–532. [Schälte et al., 2021] Schälte, Y., Alamoudi, E., and Hasenauer, J. (2021). Robust adaptive distance functions for approximate Bayesian inference on outlier-corrupted data. bioRxiv. [Scharfstein et al., 2003] Scharfstein, D. O., Daniels, M. J., and Robins, J. M. (2003). Incorporating prior beliefs about selection bias into the analysis of randomized trials with missing outcomes. Biostatistics, 4(4):495–512. [Schmitt et al., 2023] Schmitt, M., Bürkner, P.-C., Köthe, U., and Radev, S. T. (2023). Detecting model misspecification in amortized Bayesian inference with neural networks. In Springer, editor, DAGM German Conference on Pattern Recognition, pages 541–557. [Schmitt et al., 2024] Schmitt, M., Pratz, V., Köthe, U., Bürkner, P.-C., and Radev, S. T. (2024). Consistency models for scalable and fast simulation-based inference. Advances in Neural Information Processing Systems, 37:126908–126945. [Simkus and Gutmann, 2025] Simkus, V. and Gutmann, M. U. (2025). Cfmi: Flow matching for missing data imputation. arXiv preprint arXiv:2506.09258. [Smith, 2020] Smith, L. H. (2020). Selection mechanisms and their consequences: understanding and addressing selection bias. Current Epidemiology Reports, 7(4):179–189. [Streeter et al., 2017] Streeter, A. J., Lin, N. X., Crathorne, L., Haasova, M., Hyde, C., Melzer, D., and Henley, W. E. (2017). Adjusting for unmeasured confounding in nonrandomized longitudinal studies: a methodological review. Journal of Clinical Epidemiology, 87:23–34. [Strom, 2019] Strom, B. L. (2019). What is pharmacoepidemiology? Pharmacoepidemiology, pages 1–26. [Talts et al., 2020] Talts, S., Betancourt, M., Simpson, D., Vehtari, A., and Gelman, A. (2020). Validating Bayesian inference algorithms with simulation-based calibration. arXiv preprint arXiv:1804.06788. [Verma et al., 2025] Verma, Y., Bharti, A., and Garg, V. (2025). Robust simulation-based inference under missing data via neural processes. In The Thirteenth International Conference on Learning Representations. [Wang et al., 2024] Wang, Z., Hasenauer, J., and Schälte, Y. (2024). Missing data in amortized simulationbased neural posterior estimation. PLOS Computational Biology, 20(6):e1012184. [Wildberger et al., 2023] Wildberger, J., Dax, M., Buchholz, S., Green, S., Macke, J. H., and Schölkopf, B. (2023). Flow matching for scalable simulation-based inference. Advances in Neural Information Processing Systems, 36:16837–16864. [Yao and Domke, 2023] Yao, Y. and Domke, J. (2023). Discriminative calibration: Check Bayesian computation from simulations and flexible classifier. Advances in Neural Information Processing Systems, 36:36106–36131. [Zaheer et al., 2017] Zaheer, M., Kottur, S., Ravanbakhsh, S., Poczos, B., Salakhutdinov, R. R., and Smola, A. J. (2017). Deep sets. Advances in Neural Information Processing Systems, 30.

Supplementary Material Additional Results A

Absolute error (%)

3 2 1 R1

Unadjusted Round R1

Posterior Density

1.0 0.0

Round R2

Mean C2ST=0.52 p-value=1.00

Simulated round Weighted

Round R3

Mean C2ST=0.53 p-value=0.40

Bias-aware NPE

Round R4

Mean C2ST=0.51 p-value=0.90

Round R5

Mean C2ST=0.52 p-value=0.90

Mean C2ST=0.54 p-value=0.70

1.0

C2ST score (mean per bin)

0.8 0.6

Prevalence (%)

Figure S1: Additional results for the KoCo19 Study. (A) Absolute error of estimated prevalence across 1000 simulated datasets (including selection, but excluding missingness), each with five rounds of the KoCo19 Study using unadjusted counts of infections and inverse-probability weighting. The mode of the posterior is estimated with neural posterior estimation on the dataset, including missingness. (B) Classification of the posterior marginals for the real data. A classifier was trained to distinguish samples from the joint p(θ, y) and p(θ | y)p(y). Table S1: Additional results for the Framingham Study. Contraction of the posterior for the bias-aware NPE (applied to simulated observed data) and the NPE trained on full data (applied to simulated full data). Median and median absolute deviation over 1000 datasets not used during training are calculated. Parameter

Bias-aware NPE

NPE

a01 a02 a12 κ01 κ02 κ12 sex β01 sex β02 sex β12 age β01 age β02 age β12

0.59 (± 0.23) 0.93 (± 0.06) 0.33 (± 0.28) 0.77 (± 0.03) 0.97 (± 0.01) 0.79 (± 0.03) 0.97 (± 0.01) 0.99 (± 0.00) 0.89 (± 0.07) 0.99 (± 0.00) 0.99 (± 0.00) 0.99 (± 0.01)

0.91 (± 0.07) 0.92 (± 0.06) 0.32 (± 0.30) 0.95 (± 0.01) 0.96 (± 0.01) 0.73 (± 0.03) 0.99 (± 0.01) 0.98 (± 0.01) 0.22 (± 0.04) 0.99 (± 0.00) 0.99 (± 0.00) 0.99 (± 0.01)

Epoch 1

Epoch 2

Epoch 3

Epoch 4

10 2

h01

10 3 10 4 10 5 10 6 0

h02

10 1 10 3 10 5

10 2

h12

10 3 10 4 10 5 10 6 2

1 2 3 4 Follow up years since entry in epoch

age

1 0 1 2 2

sex

1 0 1 2

Transition

Naive Cox

Transition Splines IDM

Prior

Transition Bias-aware NPE

Transition

NPE

Figure S2: Additional results for the Framingham Study. Unadjusted hazards and covariate effect estimates of naive Cox, penalized likelihood approach, and the posterior for the bias-aware NPE and the NPE trained on full data (median and 95% credible intervals) are compared on the real data. The 95% quantiles of the prior are shown in grey.

A 3

r = 0.838

2 3

Estimate

r = 0.871

Ground truth

1 0 0

transm pro

r = 0.835

r = 0.889

acq pro

Estimate

Child sus

r = 0.827

0 0

Infant sus

r = 0.151

asym Adult inf

r = 0.531

asym Child inf

r = 0.017

3 4

asym Infant inf 10

r = 0.699

sym Child inf

r = 0.628

1 0

sym Infant inf

r = 0.889

Ground truth

4 2 0 0

Ground truth

Figure S3: Additional results for the PedCovid Study. (A) Recovery of the parameters on simulated data for all three selection mechanisms using bias-aware NPE. The median and median absolute deviation for each dataset are computed. (B) Recovery of the parameters on simulated data only for random selection and a unfixed background hazard α using MCMC.

Record · ID 120537 · SHA-256 8bd9b2817054c951
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.