Marginal structural model projections#
What this solves#
You have four arms, or six regimens, or a dose you want to read as a trend. Reporting every counterfactual mean separately is honest and unreadable. You want one summary: a slope in the dose, or an effect that varies with one baseline characteristic.
A marginal structural model in cleverly is a projection. It is the best approximation of the
true counterfactual surface within a working model you declared. It is not a claim that the working
model is the true causal response surface, and the estimand is well defined whether or not the
working model fits.
your situation |
what this method buys |
what it costs |
|---|---|---|
many arms, or many regimens |
one coefficient vector instead of one mean per level |
the coefficients mean what the working model says they mean. Read them as a projection |
a dose you want to summarise as a trend |
a slope, with an influence curve and an interval |
a model linear in the arm reads the arm as a dose, so non-numeric labels are refused |
effect modification by a baseline variable |
an interaction term in the working design |
the design must be full rank on the realized cells |
repeated treatment over time |
the same projection over regimen and horizon cells |
the node fluctuations are pooled. Under a link, one in-sample round of the alternation is a whole backward pass |
The projection is a fourth parameter axis. msm= cannot be combined with interventions= or
policies=, because one fluctuation solves one set of score equations, and a fit reporting
parameters from two axes would put two of them under one heading.
A worked applied analysis is in the MSM projections tutorial. It checks the saturated case against the per-arm report.
The algorithm as implemented#
For design vector \(\phi(a,V)\), nonnegative projection weight \(h(a,V)\), and the identity link,
With \(M=E\sum_a h\phi\phi^\top\) and \(r=E\sum_a h\phi\,\bar Q_P(a,W)\), the coefficient is
\(\beta=M^{-1}r\). The working design and the weights are known functions, and \(M\) must be full rank.
The package cannot inspect a callable, so it reads a declaration for each one. design_kind
declares the design, and weights_kind declares the weights. Variations gives both
declarations and their refusals.
With the identity link the clever covariate is \(h(a,V)\,\phi(a,V)/g(a\mid W)\), one column per term, so the score equation is one per coefficient rather than one per arm. The counterfactuals are still the arms. What changed is the summary. A saturated working model, with one indicator per arm, reproduces the per-arm report exactly, at the point estimate and at the influence curve.
On a continuous dose, a trapezoid sum over the declared grid doses= replaces the sum over arms.
The clever covariate then divides by the conditional treatment density at the observed dose and
at each grid dose. build_submodel in src/cleverly/estimators/targeting.py builds it.
Theory: Neugebauer and van der Laan (2007), Rosenblum and van der Laan (2010), and chapter 12 of
van der Laan and Rose (2011). Petersen et al. (2014) gives the longitudinal construction. See the
reference list.
Implementation:
msm.py,
longitudinal/msm.py,
and
fluctuation/submodel.py.
Under a link#
Three things change inside, and each is a place where the obvious generalisation is wrong.
The clever covariate reads the coefficient. It becomes
\(h(a,V)\,(dm/d\eta)\,\phi(a,V)/g(a\mid W)\), and \(dm/d\eta\) is a function of \(m(a,V;\beta)\). So the
fluctuation and the projection are solved together, alternating until both settle. It converges
fast, because the coefficient reaches the covariate only through that smooth factor. The shift falls
by a factor of \(10^{-3}\) to \(10^{-4}\) per round, and res.fluctuations["msm"].projection carries
the per-round trace.
The matrix the influence curve is premultiplied by is no longer the Gram matrix. It is
and the second term vanishes only where the working model fits, which is exactly what a projection does not promise. Dropping it is wrong in a way that no saturated-model check can see.
The remainder is second order without being zero. With the identity link a correct mechanism drives the remainder to exactly zero, and that exactness is the linearity of the estimating equation in the coefficient rather than a stronger form of double robustness. Under a link what is left is quadratic in the coefficient error. The other half is untouched: a correct outcome regression still gives exactly zero, under every link.
A saturated working model still reproduces the per-arm report through the link. The inverse link of the coefficient is the counterfactual mean to machine precision, with the influence curves related by the delta method. That is what says a link is a reparameterisation of the same counterfactual means rather than a second estimator.
The longitudinal projection#
The longitudinal MSM projects regimen, horizon, and cause-specific means onto a pooled working model. Three things differ, and each is again a place where the obvious generalisation is wrong.
The node fluctuation is pooled across the regimens. At one node the design’s columns get their rank by summing over the arms within a row, because a unit contributes its design at the arm it received. A regimen is a plan and not a value some unit took, so there is nothing to sum over within a row. Whenever the working model has no effect modifier, the design is constant down the rows, which makes a per-regimen design rank one and collapses its score equations into one. Each node therefore solves a single fluctuation over the regimens stacked, with one shared coefficient. The backward recursion is lockstep: outer over the nodes, inner over the regimens, one update, all carried forward together.
A saturated working model reproduces the per-regimen report, and not bit for bit on the
estimate. One indicator per regimen makes the stacked design exactly block diagonal, and each block
carries the loss weight the plain recursion uses. The pooled Newton convergence test and line search
are taken over all the stacked rows, so the two can stop on different iterates. On a law the sample
realises exactly, no step is taken at all and the agreement is exact. Elsewhere the gap is of the
order of the targeting tolerance tol. At the default tol=1e-10 and five folds it measured
1.1e-10, when one stacked solve stopped one Newton step away from the separate solves.
Cross-fitting pools the update. Above one fold, each outer fold runs the untargeted backward recursion of every regimen and horizon cell on its training rows. The fit stitches the held-out predictions into one out-of-fold initial fit per node and cell. One stacked fluctuation per node then targets every follower of every live cell. Its loss weight carries the out-of-fold cumulative mechanism and the projection weight. No fold solves a fluctuation, and no fold computes a coefficient.
The base result is the cross-fitted construction of Díaz, Williams, Hoffman and Schenck (2023), Section 5.2, page 852, with Theorem 3, page 853, for each cell. The arXiv v4 algorithm fits the tilting model “using all the data points in the sample”, with the cross-fitted initial nuisances. The stacked solve does not solve each cell’s own score, so no cell estimate is a Section 5.2 TMLE on its own. The coefficient error has three parts: the delta-method image of each cell’s stitched second-order remainder, the mean of the coefficient’s curve, and an empirical-process term. Each cell remainder has the product-rate form of Theorem 3.
Each fold’s initial fit is fixed given its training rows, and that controls the empirical-process
term. The targeted fit adds a fixed-dimension index: one fluctuation coefficient per term at each
node, and the projection coefficients. This is the fixed-dimension stack of the
natural-extension verdicts. The inherited conditions are the
Theorem 3 rates and bounded density ratios in every cell, a fixed number of cells, and a full-rank
realized design. With id=, the cluster rules apply as well.
evidence |
what it shows |
|---|---|
|
a saturated model reproduces the cross-fitted per-regimen report on survival, competing-risk, dynamic and three-level grids. The coefficient curve is the stacked delta-method curve. Every node’s stacked score is solved. Eight mutations fail, including a fold-local fluctuation |
|
the saturated reduction at one and five folds, under the identity and the logit link |
repeated-sampling properties on the projection law of the ordinary study, comparator-backed and paired with a projection of four cross-fitted R |
Survival, competing-risk, weighted and clustered projections have the fast-tier exact identities
only. A cross-fitted lmtp 1.5.4 regimen fit differs from this construction twice. It fluctuates
on each training fold, and it targets each regimen with its own scalar fluctuation. Upstream
commit 9996b04 classifies that training-fold update as a bug.
Under a link, one in-sample round of the alternation is a whole backward pass. The coefficient enters the covariate through the derivative of the inverse link, so each targeted regression moves with it. A targeted regression is the previous node’s regression target, so every earlier node’s learner is refit. The fixed point is stated over the whole pass, and it costs four or five passes in practice. The mechanism is free of the coefficient and is fitted once. Above one fold, a round re-solves the pooled update only, because the stitched initial fit carries no coefficient.
The projection weight and the observation weight are different objects and must stay so. The first says how the regimens are traded off inside the projection. The second tilts the population the projection is taken over. Both multiply the node’s loss weight, and only the observation weight multiplies the finished influence curve row by row. The projection is solved on the raw outcome scale, because a coefficient vector has no single scale to map back with.
Under SDR, the projection solves one estimating equation over the transforms. The
sequentially doubly robust estimator
fits no fluctuation. Each regimen and horizon cell \(c\) has its own first transform
\(\varphi_{1,c}\), and project_sdr_cells solves the working model’s equation over the stacked
cells. The coefficient curve is
with \(x\) the design, \(V\) the baseline covariates, and \(M\) the derivative of the equation in
\(\beta\). The conditional form of Lemma 1 of Díaz, Williams, Hoffman and Schenck (2023) at \(H_1\)
makes each cell’s transform unbiased given \(V\). No published result states this composition, and
the natural-extension verdicts record its steps. The identity link
always has a solution. A transform can leave \([0, 1]\), so a logit or log projection can have no
finite solution. The fit raises CapabilityError when that projection does not converge.
tests/unit/test_sdr_exact_laws.py checks the four families on the exact law, and
tests/unit/test_sdr_double_robustness.py checks the identity projection in every robustness
pattern.
Variations#
option |
what it does |
|---|---|
|
the working design vector, including effect modifiers and interactions |
|
the declaration that a |
projection weights |
how the arms or regimens are traded off inside the projection |
|
the declaration that a |
|
the clever covariate is free of the coefficient, and a correct mechanism drives the remainder to exactly zero |
|
the covariate reads the coefficient, so the fluctuation and the projection alternate. |
|
a model linear in the arm. Refused on non-numeric labels, because it would read the sort order as a dose scale nobody chose |
longitudinal |
the same projection over regimen, horizon, and cause cells, with rank checked on the actual realized design. Above one fold, the fold recursions stay untargeted and one stacked fluctuation per node pools the update |
|
each fold solves its own coefficient, since the coefficient is something the covariate reads. This removes coupling between folds, and the rows inside a fold still fit both the coefficient and the fluctuation used for that fold. The pooled score is exactly zero because each fold’s is zero at its own coefficient. This is a package extension and not the common-update CV-TMLE of Zheng and van der Laan |
point-treatment |
point-treatment fold targeting is supported. A longitudinal fit never targets inside a fold: above one fold it pools the update over every follower |
Baseline strata#
With strata=, a stratum’s coefficients \(\beta_s\) are the projection on the law given \(S=s\).
The projection loss separates by stratum, so \((\beta_s)_s\) are the coefficients of one working
model with the expanded design \((I(S=s)\varphi)_s\). The table gives where each estimate comes
from.
link |
stratum coefficients |
marginal coefficients |
|---|---|---|
identity, finite arms or a dose grid |
one block \(I(S=s)H/P_n(S=s)\) per stratum, in one pooled fluctuation |
the same fluctuation; its score is the weighted sum of the blocks |
log, logit |
a nested fluctuation of the blocks \(I(S=s)H(\beta_s)/P_n(S=s)\), alternating with the projection on each stratum’s rows, recorded on |
the unstratified alternation, unchanged |
The two links differ because a linked block reads its own stratum’s coefficients. One fluctuation cannot be the stratum fluctuation for every stratum and the marginal fluctuation at once, so a linked fit carries two. Each estimate is asymptotically linear with its own curve, so the stacked-curve covariance of any marginal and stratum coefficients stays valid.
Before any learner, the fit restricts the evaluated design to each stratum’s rows and reruns the
shipped rank rule. A design singular inside a stratum refuses with the stratum named first. A
term that is constant there, such as the stratum column itself, and a term that copies another
there, such as a:S in stratum S = 1, are both singular. A link does not change the rank,
because \(dm/d\eta>0\).
tests/unit/test_stratified_msm_exact.py checks the linked coefficients and curves against the
expanded-design fit, and the curve against a complex-step Gateaux derivative of \(\beta_s\).
tests/unit/test_stratified_continuous_msm_exact.py checks the dose case against the expanded
design under the identity and logit links. The
stratified identity-link MSM study and the
stratified incremental MSM study measure the
stratum coefficients against exact laws.
The package cannot inspect a callable, so it reads the declarations of the design and of the projection weight. The table gives each place that checks both declarations.
checked by |
when |
|---|---|
the |
when you declare the model |
the |
before the first learner and before the first call to the design. A copied or modified model can carry a declaration that the fit refuses. The |
|
before it recomputes an estimate. Each sweep that calls it, such as |
|
before they call the design or the weight. A direct call and the simulated-confounding replay both meet this check |
Each refusal raises CapabilityError. The
scope and refusals page defines each reason code.
refused |
code |
reason |
|---|---|---|
a |
|
the package cannot tell a fixed design from a design computed from the sample. The message asks for |
|
|
for a population-law target that centres at \(E_P[W]\), the curve needs the derivative through that mean. The reported curve omits it. Conditioning on a centre learned from these same rows defines a different target, with no validated interval here |
a |
|
the package cannot tell a fixed weight from a weight computed from the sample. The message asks for |
|
|
the weight is a functional of \(P\), so the influence curve needs a further term for the pathwise derivative through the estimated mechanism. The reported curve does not have it. This is the same argument that gives an incremental intervention its own axis |
an array in place of a |
|
an array is one evaluation of the weight, and nothing shows that the sample did not set it |
a continuous dose with missing outcomes ( |
|
the continuous-dose covariate that the package builds must then also divide by the response or intermediate mechanism at each grid dose. The package does not predict either mechanism at a grid dose. |
A malformed input raises DataError. The table gives the four cases.
input |
example |
|---|---|
a |
an array or a float. The message gives the point-treatment signature and the longitudinal signature |
a |
|
a |
an array or |
|
the declaration describes a weight that the model does not have |
The simulated-confounding replay reports each of these inputs as CapabilityError. Its function
validate_fixed_replay converts each DataError that its checks raise into a CapabilityError.
A longitudinal MSM stores evaluated arrays and evidence that its source declarations passed. A projection built by hand has no such evidence, and it refuses truncation replay. The evidence survives a save and a load.
Except under the rule below, every call that recomputes an estimate from a TMLE result with an
undeclared function refuses with the undeclared message. retarget() checks the declarations first,
so each sweep that calls it refuses, for example truncation_curve(). refute() refits through
fit(), which checks them too. The simulated-confounding replay refuses in MSMSet.evaluate,
before the design or the weight runs.
One rule reads a design with no declaration as known. MSM.linear leaves design_kind as None.
When the design has the exact type that MSM.linear builds, the model reads as
design_kind="known". Its result still recomputes when its weight is uniform or declared
"known". That design is a fixed function of the arm and the named covariates.
A result with any other undeclared design refuses every recomputation. The from_linear flag
is not a declaration, because a user can set it on any design. A subclass of the MSM.linear design
type is not a declaration either.
test file |
what it pins |
|---|---|
|
for the weight: the stored interval, the refusals, the evaluator checks, a declared control, and a mutation that removes the check |
|
for the design: the same five items, the exact-type rule, a malformed model at each fit entry, the replay refusal before any call to the design or the weight, and an |
An LTMLE MSM result holds the evaluated design and weight arrays, not the MSM object. Its
truncation_curve() reuses those arrays, runs no user function, and reports point estimates only.
It therefore needs no declaration check. tests/unit/test_longitudinal_truncation_refit.py pins the
columns of an LTMLE truncation curve, and it replays an MSM result from the stored arrays.
dataclasses.replace(model, design=new) keeps the design_kind of model, as it keeps
weights_kind. Pass design_kind again when you replace the design function. A model from
MSM.linear holds None, so replace refuses a new design that has no design_kind. The test
test_replacing_the_shorthand_design_drops_its_declaration pins this.
A weight computed from the sample and declared "known" still fits, because the declaration is
your statement. tests/unit/test_msm_projection_weights.py measures the cost on an exact law, with
the arm share as the weight. The reported curve is the efficient influence function of the
fixed-weight coefficient. The table gives its standard error over that of the estimated-weight
coefficient. RM13 records the defect.
coefficient |
standard-error ratio |
|---|---|
|
0.742 |
|
0.913 |
|
1.000, because the share term vanishes for this coefficient on this design |
A design computed from the sample and declared "known" also still fits.
tests/unit/test_msm_design_declaration.py measures the cost on an exact law with the design
\([1, a, W-c]\). Here \(c\) is the sample mean of \(W\), which is 0.7. The reported curve is the efficient
influence function of the fixed-centre coefficient. The table gives its standard error over that of
the estimated-centre coefficient. RM27 records
the defect.
coefficient |
standard-error ratio |
|---|---|
|
0.893 |
|
1.000, because the centre term moves only the intercept |
|
1.000, because the centre term moves only the intercept |
The omitted term is \(\beta_W(W-c)\) in the intercept curve. It changes the variance of the intercept
by \(\beta_W^2\,\mathrm{Var}(W)+2\beta_W\,\mathrm{Cov}(D,W)\), where \(D\) is the reported curve. That
change can have either sign, and it vanishes at \(\beta_W=0\). The class
TestAKnownDesignKeepsItsInterval fits two controls, a design with the fixed centre 0.7 and
MSM.linear. Both keep their intervals. The code cannot detect the false "known" declaration: the
fixed-centre curve and the witness curve are bitwise equal.
A one-shot non-identity-link fit is also refused. The derivative of the inverse link depends on the coefficient, so a single pass would report a standard error for an equation it did not solve. The link is supported. What is refused is skipping the alternation it needs.
Validation issues special to this method#
A saturated oracle proves nothing about a projection. A saturated working model fits, which is exactly what a projection does not promise. A check against the per-arm means therefore cannot see a component that vanishes when the model fits, such as the projection weight or a link’s curvature term. The oracle’s working model is therefore deliberately not saturated: three coefficients against six covariate-and-arm cells at one node, and three coefficients against twelve regimen cells longitudinally. Both choices are asserted on the law itself, so they are shown to be load-bearing rather than claimed to be.
Uniform weights are equally disqualifying. With a constant weight the design is orthogonal and the coefficient collapses to the marginal ATE identically, so code that reported the ATE under the name of an MSM coefficient would pass every check. The oracle’s weights are deliberately not uniform.
A comparison is not analytic. Under a link the oracle solves its own normal equations by a fixed number of Newton steps with no convergence test, because a functional that branched on a convergence test could not be differentiated by a complex step at all. A separate check requires that doubling the step count moves neither the value nor the curve.
Coefficient parity with R ltmleMSM is not an acceptance check. That package uses a
quasibinomial working-model projection. cleverly declares an outcome-scale weighted least-squares
projection. Raw coefficient parity would compare different estimands.
Seven mutations were applied to the longitudinal projection and the tests watched. Three passed on the first try and exposed real gaps. That is why deliberate mutation is part of the development guidance rather than an optional extra.
where to read the evidence |
what is there |
|---|---|
the Gateaux comparison, the nonzero remainder, the saturated reduction, the link-specific submodel tests, and the continuous-dose density-ratio scores |
|
the non-saturated, nonuniform projection law, the exact pooled-design and loss-weight checks, and the rank-refusal tests |
|
truth, R |
|
projected joint R |
|
projected joint cross-fitted R |
The registered rows cover fixed projections. In sample they cover the identity link. Under cross-fitting they cover the identity link and the calibration of a logit link. The continuous-dose test still uses a linear truth, and a nonlinear continuous-dose Gateaux oracle remains absent.