Add unconditional variance for averaged marginal effects by pedrohcgs · Pull Request #1737 · vincentarelbundock/marginaleffects

pedrohcgs · 2026-07-08T01:37:04Z

Hi Vincent,

I am opening this as a draft because it overlaps with the motivation in #1240 and the prototype in #1711, but it takes a different implementation path. Rather than adding a standalone post-processing helper, this tries to make unconditional variance a regular vcov option inside the existing marginaleffects workflow.

The practical motivation is my own workflow with etwfe. I have been using etwfe for applied work where I need averaged/aggregated marginal effects, and the usual conditional standard errors miss an important part of the uncertainty because the empirical covariate distribution is treated as fixed. This PR is an attempt to make that workflow reliable directly in marginaleffects, instead of carrying around local patches or ad hoc post-processing code.

Tagging @grantmcdermott because this came up for me through etwfe, and he may have useful thoughts on whether the API and scope make sense for that use case.

What this adds

The new user-facing entry points are:

vcov = "unconditional"
vcov = unconditional(~cluster)

The feature accounts for sampling variation in the empirical covariate distribution when computing averaged or aggregated predictions, comparisons, and slopes.

I tried to keep this additive and conservative. The implementation reuses the existing prediction/comparison plans, aggregation logic, Jacobian machinery, and result finalization instead of building a separate parallel path.

Current scope

Supported workflows include averaged or aggregated:

predictions
comparisons
slopes
one-way clustered unconditional inference
valid subsets of original model data
counterfactual grids that preserve row mappings

The code intentionally rejects cases where the empirical-distribution influence function is not well-defined, or where important uncertainty components are not represented.

Currently unsupported:

unit-level effects without aggregation
raw hypotheses() calls on model objects
multiple-imputation objects
survey-design models
posterior-draw models
mixed-effects models
multi-equation and multinomial models
unsupported coxph prediction types
fixed-effect model cases where fixed-effect uncertainty is not represented

Validation

The new tests cover average predictions, comparisons, and slopes; robust and clustered unconditional variance; degrees-of-freedom behavior; transformations; by; weights; offsets; custom hypotheses; scalar and pairwise comparisons; counterfactual grids; valid subsets of original model data; high-dimensional covariate settings; optional fixest, etwfe, survival, AER, survey, and mice paths; and explicit errors for unsupported model/data structures.

Local checks before opening the draft:

focused unconditional variance tests pass
R CMD build passes
full R CMD check passes, including vignettes and PDF manual

I do not expect this to be the final shape of the API or internals. The main thing I would like feedback on is whether integrating this through vcov is the right direction, and whether the current conservative set of supported cases is a reasonable starting point.

vincentarelbundock · 2026-07-08T10:05:25Z

Pedro,

Thanks so much for this. I've been wanting to get this for a long time but never got to it.

Your API is much better than the other PR. This should definitely go through the vcov argument.

I'm not sure we want to export a unconditional function, since that's a generic word that may be used by other packages. But that's a small detail.

I'm traveling now but will review as soon as I find some time.

Thanks!

grantmcdermott · 2026-07-09T13:43:27Z

This looks great! I'm also out traveling / on vacation, so won't be able to review carefully for a couple of weeks. But happy to do so when I get back if there's still need.

vincentarelbundock · 2026-07-13T14:15:39Z

Thanks again for the implementation @pedrohcgs !

I made some changes:

Reorganized the code and split into separate files math vs. bookeeping.
Renamed unconditional() to vcovUnconditional() and mimicked the sandwich interface.
Moved from blacklist to whitelist approach to be more conservative in the models and estimands we support.
Added some closed-form expressions from Hansen-Overgaard instead of jackknife for a few commands.
More tests.

No rush for @grantmcdermott, but I'd really love your feedback on this, especially on r/R/vcov_unconditional.R

@ngreifer also expressed interest in the other thread, so tagging here in case you want to check this out.

And @snhansen, I'd be especially grateful for a review by you, since this implements several things from (or related to) your paper.

Edit: one more thing that would be nice would be a .do file that produces "golden" values that we could dump to a CSV and check against in the test suite.

ngreifer · 2026-07-13T14:41:26Z

This is incredible and much-needed work, thank you so much @pedrohcgs for taking it on and @vincentarelbundock for pushing it through. I'm eager to validate this against the implementation in adrftools.

For multiply imputed data, I think it makes sense to just estimate the unconditional variance in each dataset and pool the final SEs together.

For survey objects, I use this code in adrftools. It's a simple matter of extracting the influence function from a svyglm object.

vincentarelbundock · 2026-07-13T18:05:47Z

Thanks a lot for this, @ngreifer , super helpful!

I added support for multiple imputation and the survey package. I (read: codex) also created some comparisons to your adrftools package.

Things look pretty good at 1e-6 tolerance.

Let me know if you think of other cases to add, or if you spot any issues:

https://github.com/pedrohcgs/marginaleffects/blob/c9102405b5784be29409356e2740738ad3386230/r/inst/tinytest/test-vcov_unconditional-adrftools.R

ngreifer · 2026-07-13T18:54:11Z

Looks great to me! I don't think there should be any trouble supporting multinomial/ordinal models, as long as they have methods of extracting the influence functions. I know multinom_weightit and ordinal_weightit objects from WeightIt work in this framework with adrftools, so I expect them to work here, too. glm_weightit objects inherit from glm so they should already be supported.

Also just FYI, in survey is the function svypredmeans(), which also generates predicted means and a corrected variance, but it's one of the incomplete corrected variances that isn't robust to model misspecification. Maybe worth noting for the few users who will notice differences between the results. The marginaleffects unconditional variance should be preferred.

grantmcdermott · 2026-07-31T20:49:47Z

No rush for @grantmcdermott, but I'd really love your feedback on this, especially on r/R/vcov_unconditional.R

Finally took a quick pass at this. As you know, my main interest is seeing how well this plays with etwfe. I think the answer is "nicely", but please note that I haven't vetted my results against known benchmarks. (@vincentarelbundock did you ever get around to adding a canonical .do file of known results?)

In the absence of known results, the thing I'm most concerned about is making sure that there's no "double counting" of the VCOV via etwfe's consecutive two-step calls: etwfe() -> emfx(). To that end, I tested passing through (and not) a cluster variable into the first etwfestep (powered by fixest::feols). It looks like the output is stable, so long as I call vcovUnconditional(cluster = ~...) in the second emfx step.

library(etwfe)
library(marginaleffects)
data("mpdta", package = "did")

mod = etwfe(
  fml  = lemp ~ lpop,
  tvar = year,
  gvar = first.treat,
  data = mpdta,
  vcov = ~countyreal  # clustered vcov
)
mod2 = etwfe(
  fml  = lemp ~ lpop,
  tvar = year,
  gvar = first.treat,
  data = mpdta,
  # vcov = ~countyreal  # iid
)

emfx(mod, vcov = vcovUnconditional(cluster = ~countyreal)) # ATT from clustered 1st step
#> 
#>  .Dtreat Estimate Std. Error     z Pr(>|z|)    S   2.5 %  97.5 %
#>     TRUE  -0.0506      0.013 -3.89   <0.001 13.3 -0.0761 -0.0251
#> 
#> Term: .Dtreat
#> Type: response
#> Comparison: TRUE - FALSE
emfx(mod2, vcov = vcovUnconditional(cluster = ~countyreal)) # ATT from iid 1st step
#> 
#>  .Dtreat Estimate Std. Error     z Pr(>|z|)    S   2.5 %  97.5 %
#>     TRUE  -0.0506      0.013 -3.89   <0.001 13.3 -0.0761 -0.0251
#> 
#> Term: .Dtreat
#> Type: response
#> Comparison: TRUE - FALSE

So, that's good news.

On the other hand, I wonder whether we couldn't do something smart whereby, if the user passed a clustered VCOV in the first stage, then vcovUnconditional() would automatically recognize and pass this cluster variable down? It might be too much effort and peculiar to fixest... For the record, we get the same result for both models here, which I assume is "H0"?

emfx(mod, vcov = vcovUnconditional()) # ATT from clustered 1st step
#> 
#>  .Dtreat Estimate Std. Error     z Pr(>|z|)   S 2.5 % 97.5 %
#>     TRUE  -0.0506     0.0404 -1.25     0.21 2.3 -0.13 0.0286
#> 
#> Term: .Dtreat
#> Type: response
#> Comparison: TRUE - FALSE
emfx(mod2, vcov = vcovUnconditional()) # ATT from iid 1st step
#> 
#>  .Dtreat Estimate Std. Error     z Pr(>|z|)   S 2.5 % 97.5 %
#>     TRUE  -0.0506     0.0404 -1.25     0.21 2.3 -0.13 0.0286
#> 
#> Term: .Dtreat
#> Type: response
#> Comparison: TRUE - FALSE

P.S. For posterity, confirming that vcovUnconditional() is yielding (slightly) different results from the default emfx() behaviour.

emfx(mod)
#> 
#>  .Dtreat Estimate Std. Error     z Pr(>|z|)    S   2.5 %  97.5 %
#>     TRUE  -0.0506     0.0125 -4.05   <0.001 14.3 -0.0751 -0.0261
#> 
#> Term: .Dtreat
#> Type: response
#> Comparison: TRUE - FALSE

And also in event study form:

emfx(mod, "event")
#> 
#>  event Estimate Std. Error     z Pr(>|z|)    S   2.5 %   97.5 %
#>      0  -0.0332     0.0134 -2.48    0.013  6.3 -0.0594 -0.00702
#>      1  -0.0573     0.0171 -3.34   <0.001 10.2 -0.0910 -0.02373
#>      2  -0.1379     0.0308 -4.48   <0.001 17.0 -0.1982 -0.07753
#>      3  -0.1095     0.0323 -3.39   <0.001 10.5 -0.1729 -0.04620
#> 
#> Term: .Dtreat
#> Type: response
#> Comparison: TRUE - FALSE
emfx(mod, "event", vcov = vcovUnconditional(cluster = ~countyreal))
#> 
#>  event Estimate Std. Error     z Pr(>|z|)    S   2.5 %  97.5 %
#>      0  -0.0332     0.0135 -2.46   0.0141  6.2 -0.0597 -0.0067
#>      1  -0.0573     0.0180 -3.19   0.0014  9.5 -0.0925 -0.0222
#>      2  -0.1379     0.0339 -4.07   <0.001 14.4 -0.2042 -0.0715
#>      3  -0.1095     0.0323 -3.39   <0.001 10.5 -0.1728 -0.0463
#> 
#> Term: .Dtreat
#> Type: response
#> Comparison: TRUE - FALSE

grantmcdermott · 2026-07-31T21:10:20Z

Update: It turns out I still have access to a Stata license through work. I just tested against jwdid (guide) and the resulting SEs are very close... but not exactly the same as I'm getting.

. qui jwdid lemp lpop, ivar(county) time(year) gvar(first_treat) method(regress)

. estat simple, vce(unconditional)
                           (Std. err. adjusted for 500 clusters in countyreal)
------------------------------------------------------------------------------
             |            Unconditional
             | Coefficient  std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
      simple |   -.050627   .0130756    -3.87   0.000     -.076317   -.0249371
------------------------------------------------------------------------------

. estat event, vce(unconditional)
                           (Std. err. adjusted for 500 clusters in countyreal)
------------------------------------------------------------------------------
             |            Unconditional
             | Coefficient  std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
   __event__ |
          0  |  -.0332122   .0136055    -2.44   0.015    -.0599432   -.0064812
          1  |  -.0573456   .0180591    -3.18   0.002    -.0928269   -.0218644
          2  |  -.1378704   .0340633    -4.05   0.000    -.2047956   -.0709451
          3  |  -.1095395   .0324765    -3.37   0.001    -.1733469    -.045732
------------------------------------------------------------------------------

Compare my results from above (coercing to regular data.frame for more decimals):

> emfx(mod, vcov = vcovUnconditional(cluster = ~countyreal)) |> as.data.frame() 
     term     contrast .Dtreat    estimate  std.error statistic      p.value  s.value    conf.low   conf.high
    <chr>        <chr>   <lgl>       <num>      <num>     <num>        <num>    <num>       <num>       <num>
1 .Dtreat TRUE - FALSE    TRUE -0.05062703 0.01299948 -3.894542 9.838462e-05 13.31121 -0.07610555 -0.02514851

> emfx(mod, 'event', vcov = vcovUnconditional(cluster = ~countyreal)) |> as.data.frame()
     term     contrast event    estimate  std.error statistic      p.value   s.value    conf.low    conf.high
    <chr>        <chr> <num>       <num>      <num>     <num>        <num>     <num>       <num>        <num>
1 .Dtreat TRUE - FALSE     0 -0.03321220 0.01352628 -2.455383 1.407346e-02  6.150879 -0.05972323 -0.006701176
2 .Dtreat TRUE - FALSE     1 -0.05734565 0.01795403 -3.194027 1.403028e-03  9.477240 -0.09253489 -0.022156405
3 .Dtreat TRUE - FALSE     2 -0.13787039 0.03386512 -4.071162 4.677921e-05 14.383773 -0.20424480 -0.071495972
4 .Dtreat TRUE - FALSE     3 -0.10953946 0.03228746 -3.392631 6.922473e-04 10.496425 -0.17282172 -0.046257194

pedrohcgs · 2026-07-31T23:02:59Z

I will test this here with Stata, too. I have played enough with jwdid enough. Plus, I can also validate this if it is numerically the same as csdid when they must match.

…

---------------------------------------------------- *Pedro H. C. Sant'Anna* *https://psantanna.com <https://psantanna.com>*

---------------------------------------------------- Warning: This email may contain confidential or privileged information intended only for the use of the individual or entity to whom it is addressed. If you are not the intended recipient, please understand that any disclosure, copying, distribution, or use of the contents of this email is strictly prohibited. ----------------------------------------------------

On Fri, Jul 31, 2026 at 17:10 Grant McDermott ***@***.***> wrote: *grantmcdermott* left a comment (vincentarelbundock/marginaleffects#1737) <#1737 (comment)> Update: It turns out I still have access to a Stata license through work. I just tested against jwdid (guide <https://friosavila.github.io/app_metrics/app_metrics8.html>) and the resulting SEs are very close... but not *exactly* the same as I'm getting. . qui jwdid lemp lpop, ivar(county) time(year) gvar(first_treat) method(regress) . estat simple, vce(unconditional) (Std. err. adjusted for 500 clusters in countyreal)------------------------------------------------------------------------------ | Unconditional | Coefficient std. err. t P>|t| [95% conf. interval]-------------+---------------------------------------------------------------- simple | -.050627 .0130756 -3.87 0.000 -.076317 -.0249371------------------------------------------------------------------------------ . estat event, vce(unconditional) (Std. err. adjusted for 500 clusters in countyreal)------------------------------------------------------------------------------ | Unconditional | Coefficient std. err. t P>|t| [95% conf. interval]-------------+---------------------------------------------------------------- __event__ | 0 | -.0332122 .0136055 -2.44 0.015 -.0599432 -.0064812 1 | -.0573456 .0180591 -3.18 0.002 -.0928269 -.0218644 2 | -.1378704 .0340633 -4.05 0.000 -.2047956 -.0709451 3 | -.1095395 .0324765 -3.37 0.001 -.1733469 -.045732------------------------------------------------------------------------------ Compare my results from above (coercing to regular data.frame for more decimals): > emfx(mod, vcov = vcovUnconditional(cluster = ~countyreal)) |> as.data.frame() term contrast .Dtreat estimate std.error statistic p.value s.value conf.low conf.high <chr> <chr> <lgl> <num> <num> <num> <num> <num> <num> <num>1 .Dtreat TRUE - FALSE TRUE -0.05062703 0.01299948 -3.894542 9.838462e-05 13.31121 -0.07610555 -0.02514851 > emfx(mod, 'event', vcov = vcovUnconditional(cluster = ~countyreal)) |> as.data.frame() term contrast event estimate std.error statistic p.value s.value conf.low conf.high <chr> <chr> <num> <num> <num> <num> <num> <num> <num> <num>1 .Dtreat TRUE - FALSE 0 -0.03321220 0.01352628 -2.455383 1.407346e-02 6.150879 -0.05972323 -0.0067011762 .Dtreat TRUE - FALSE 1 -0.05734565 0.01795403 -3.194027 1.403028e-03 9.477240 -0.09253489 -0.0221564053 .Dtreat TRUE - FALSE 2 -0.13787039 0.03386512 -4.071162 4.677921e-05 14.383773 -0.20424480 -0.0714959724 .Dtreat TRUE - FALSE 3 -0.10953946 0.03228746 -3.392631 6.922473e-04 10.496425 -0.17282172 -0.046257194 — Reply to this email directly, view it on GitHub <#1737?email_source=notifications&email_token=ABE734Y77WBS5YGGLWGDECD5HUDNHA5CNFSNUABFM5UWIORPF5TWS5BNNB2WEL2JONZXKZKDN5WW2ZLOOQXTKMJUG42DSNBRGA42M4TFMFZW63VHNVSW45DJN5XKKZLWMVXHJLDGN5XXIZLSL5RWY2LDNM#issuecomment-5147494109>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABE734YT53VGI4RZFN4LEYD5HUDNHAVCNFSNUABFKJSXA33TNF2G64TZHM2DAMZWGI2TOMJVHNEXG43VMU5TIOBTGM2TKNZWGE42C5QC> . You are receiving this because you were mentioned.Message ID: ***@***.***>

Add unconditional variance estimates

0927694

vincentarelbundock added 5 commits July 8, 2026 14:06

vcov="unconditional": blacklist -> whitelist

498eebf

unconditional flag for simplicity

799861c

vcov_unconditional: df + keep survival warning

8c7298a

simplify and consolidate unconditional vcov validation

783fd6a

split math and admin

b96f35f

vincentarelbundock added 12 commits July 11, 2026 23:20

Merge branch 'main' into codex/unconditional-variance

22a8a9b

more sanitation simplification

53aa53d

df argument moves from vcov_unconditional to top level function

7a1d98a

unconditional: comment + reorder + hansen and overgaard exact

1b074d8

references

ee6e2a2

unconditional() -> vcovUnconditional()

3b2c9bc

vcovUnconditional: stricter fixest

f21ad34

misc code review

f4a0c0d

remove HC2-HC5 which were added via yolo

fc671ae

restore golden tests

0325277

strengthen unconditional variance validation and tests

89c4628

document unconditional HC1 convention

2de6e4d

vincentarelbundock added 3 commits July 13, 2026 13:36

unconditional-adrftools test

8853a5f

vcovUnconditional() supports multiple imputation

e4b7d54

Support svyglm in unconditional variance

c910240

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add unconditional variance for averaged marginal effects - #1737

Add unconditional variance for averaged marginal effects#1737
pedrohcgs wants to merge 21 commits into
vincentarelbundock:mainfrom
pedrohcgs:codex/unconditional-variance

pedrohcgs commented Jul 8, 2026 •

edited

Loading

Uh oh!

vincentarelbundock commented Jul 8, 2026

Uh oh!

grantmcdermott commented Jul 9, 2026

Uh oh!

vincentarelbundock commented Jul 13, 2026 •

edited

Loading

Uh oh!

ngreifer commented Jul 13, 2026

Uh oh!

vincentarelbundock commented Jul 13, 2026

Uh oh!

ngreifer commented Jul 13, 2026

Uh oh!

grantmcdermott commented Jul 31, 2026

Uh oh!

grantmcdermott commented Jul 31, 2026

Uh oh!

pedrohcgs commented Jul 31, 2026 via email

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

Conversation

pedrohcgs commented Jul 8, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

What this adds

Current scope

Validation

Uh oh!

vincentarelbundock commented Jul 8, 2026

Uh oh!

grantmcdermott commented Jul 9, 2026

Uh oh!

vincentarelbundock commented Jul 13, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

ngreifer commented Jul 13, 2026

Uh oh!

vincentarelbundock commented Jul 13, 2026

Uh oh!

ngreifer commented Jul 13, 2026

Uh oh!

grantmcdermott commented Jul 31, 2026

Uh oh!

grantmcdermott commented Jul 31, 2026

Uh oh!

pedrohcgs commented Jul 31, 2026 via email

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

pedrohcgs commented Jul 8, 2026 •

edited

Loading

vincentarelbundock commented Jul 13, 2026 •

edited

Loading