Split out of #637, which establishes the underlying fact. This issue covers the alternative estimator; correcting inference for the current one is a sibling issue (#645).
What is missing
The methodology article describes the model the method is meant to fit: biomarkers are conditionally independent given the time since infection, so their response densities multiply inside a single integral over that shared latent time.
$$
p(Y_1=y_1, Y_2=y_2) = \int_t p(Y_1=y_1 \mid T=t), p(Y_2=y_2 \mid T=t), p_\lambda(T=t), dt
$$
The code does not compute this. log_likelihood() loops over antigen_isos, calls f_dev() once per biomarker, and sums, so each biomarker gets its own integral over its own latent time. Conditional independence given t does not survive integrating t out, so one integral over a shared t is not the product of two separate integrals -- and #637 shows the identity holds exactly on real data, confirming it is the product that is being computed.
That the code implements the weaker model is a real modelling gap, independent of the inference problem in #645. Even a perfectly calibrated standard error for the current estimator is a standard error for a different model than the one documented.
Why the joint version should be better, and why that needs measuring
Using the shared latent time is expected to be more efficient: it exploits the fact that a person has one infection history rather than one per assay, so a high IgA reading and a high IgG reading in the same person reinforce each other's evidence about t rather than counting as two independent observations.
Expected does not mean established. The joint estimator is only better if the conditional-independence-given-t assumption holds; where biomarkers are correlated beyond what the shared t explains, it can be worse than the composite one, which is more robust to exactly that misspecification. So this issue should produce a comparison, not just an implementation -- point estimates, standard errors, and coverage against simulated data where the truth is known. sim_pop_data_multi() plus analyze_sims() already provide that harness, and the methodology article's validation section already exercises it for the single-biomarker case.
Where the work is
The integration happens in C: f_dev0() (R/f_dev.r) marshals one biomarker's values and parameters -- y, nu, eps, y.low, y.high, y1, alpha, d -- into .C("negloglik", ...). A joint likelihood needs negloglik (src/serocalc.c) to accept several biomarkers' data and parameter sets per subject and perform a single quadrature over t with the per-biomarker densities multiplied inside the integrand.
Two things worth settling early, since both affect the interface rather than just the internals:
- Curve-parameter draws must be coupled across biomarkers. The posterior draws are indexed by
iter; the joint model needs the same subject's biomarkers evaluated under a coherent draw, not independently sampled ones. Whether the existing iter alignment already guarantees this needs checking before the C work starts.
- Cost. One integral per subject over a shared
t, with a product of densities inside, is not obviously cheaper or dearer than several separate integrals; it should be benchmarked, since est_seroincidence_by() runs this in an optimizer loop across many strata.
Relationship to the planned multivariate work
This is not the multivariate model described under #sec-multivariate in the methodology article (Kwan Ho Lee's Stan work, UCD-SERG/shigella#13). That work relaxes conditional independence given t, allowing genuine covariance among biomarkers. This issue implements conditional independence given t correctly -- the model the package already claims.
They stack rather than compete: this is the model the current documentation describes, and that work is the generalization of it. Doing this one first also gives that work a like-for-like baseline to be compared against.
Backwards compatibility
Keep the current summed-marginals estimator available as an explicit option, and do not remove it. It is not merely legacy: as noted above, an independence likelihood is more robust to residual between-biomarker correlation than the joint one, so it remains a defensible choice on its own merits as well as for reproducing previously published estimates.
Suggested shape, to be settled in review: an argument on log_likelihood() / est_seroincidence() / est_seroincidence_by() selecting how biomarkers are combined, with the current behaviour reachable by name and numerically identical to today's output. Whether the default changes is a separate call from whether the option exists.
Note the interaction with #645: the two estimators need different variance estimators. The composite one needs the sandwich; the joint one, being a genuine likelihood, has a valid H^-1 (up to any sampling-design clustering, which is a separate level). Whatever selects the estimator should therefore also select the matching default variance, or it becomes possible to request a mismatched pair.
Acceptance
Split out of #637, which establishes the underlying fact. This issue covers the alternative estimator; correcting inference for the current one is a sibling issue (#645).
What is missing
The methodology article describes the model the method is meant to fit: biomarkers are conditionally independent given the time since infection, so their response densities multiply inside a single integral over that shared latent time.
The code does not compute this.
log_likelihood()loops overantigen_isos, callsf_dev()once per biomarker, and sums, so each biomarker gets its own integral over its own latent time. Conditional independence giventdoes not survive integratingtout, so one integral over a sharedtis not the product of two separate integrals -- and #637 shows the identity holds exactly on real data, confirming it is the product that is being computed.That the code implements the weaker model is a real modelling gap, independent of the inference problem in #645. Even a perfectly calibrated standard error for the current estimator is a standard error for a different model than the one documented.
Why the joint version should be better, and why that needs measuring
Using the shared latent time is expected to be more efficient: it exploits the fact that a person has one infection history rather than one per assay, so a high IgA reading and a high IgG reading in the same person reinforce each other's evidence about
trather than counting as two independent observations.Expected does not mean established. The joint estimator is only better if the conditional-independence-given-
tassumption holds; where biomarkers are correlated beyond what the sharedtexplains, it can be worse than the composite one, which is more robust to exactly that misspecification. So this issue should produce a comparison, not just an implementation -- point estimates, standard errors, and coverage against simulated data where the truth is known.sim_pop_data_multi()plusanalyze_sims()already provide that harness, and the methodology article's validation section already exercises it for the single-biomarker case.Where the work is
The integration happens in C:
f_dev0()(R/f_dev.r) marshals one biomarker's values and parameters --y,nu,eps,y.low,y.high,y1,alpha,d-- into.C("negloglik", ...). A joint likelihood needsnegloglik(src/serocalc.c) to accept several biomarkers' data and parameter sets per subject and perform a single quadrature overtwith the per-biomarker densities multiplied inside the integrand.Two things worth settling early, since both affect the interface rather than just the internals:
iter; the joint model needs the same subject's biomarkers evaluated under a coherent draw, not independently sampled ones. Whether the existingiteralignment already guarantees this needs checking before the C work starts.t, with a product of densities inside, is not obviously cheaper or dearer than several separate integrals; it should be benchmarked, sinceest_seroincidence_by()runs this in an optimizer loop across many strata.Relationship to the planned multivariate work
This is not the multivariate model described under
#sec-multivariatein the methodology article (Kwan Ho Lee's Stan work, UCD-SERG/shigella#13). That work relaxes conditional independence givent, allowing genuine covariance among biomarkers. This issue implements conditional independence giventcorrectly -- the model the package already claims.They stack rather than compete: this is the model the current documentation describes, and that work is the generalization of it. Doing this one first also gives that work a like-for-like baseline to be compared against.
Backwards compatibility
Keep the current summed-marginals estimator available as an explicit option, and do not remove it. It is not merely legacy: as noted above, an independence likelihood is more robust to residual between-biomarker correlation than the joint one, so it remains a defensible choice on its own merits as well as for reproducing previously published estimates.
Suggested shape, to be settled in review: an argument on
log_likelihood()/est_seroincidence()/est_seroincidence_by()selecting how biomarkers are combined, with the current behaviour reachable by name and numerically identical to today's output. Whether the default changes is a separate call from whether the option exists.Note the interaction with #645: the two estimators need different variance estimators. The composite one needs the sandwich; the joint one, being a genuine likelihood, has a valid
H^-1(up to any sampling-design clustering, which is a separate level). Whatever selects the estimator should therefore also select the matching default variance, or it becomes possible to request a mismatched pair.Acceptance
lambda: bias, standard error, and interval coverage.