vincentarelbundock / marginaleffects

R package to compute and plot predictions, slopes, marginal means, and comparisons (contrasts, risk ratios, odds, etc.) for over 100 classes of statistical and ML models. Conduct linear and non-linear hypothesis tests, or equivalence tests. Calculate uncertainty estimates using the delta method, bootstrapping, or simulation-based inference
https://marginaleffects.com
Other
392 stars 43 forks source link

Raise error when `hypotheses()` is called more than once #1102

Closed poliquin closed 2 months ago

poliquin commented 2 months ago

Description of the bug

I am trying to use multiple calls to hypotheses to test for differences in treatment effects across groups and cohorts, and to calculate averages of these effects. The problem is basically this:

  1. Use slopes to get marginal effects. (In my example, treatment effects by two covariates with two levels: sex and cohort.)
  2. Use hypotheses with the object returned by Step 1 and a matrix to estimate/test two linear combinations of the marginal effects. (In my example, whether the treatment effect varies in one of the covariates: sex differences in treatment effects within cohorts.)
  3. Use hypotheses again on the object returned by Step 2 to estimate/test the average of the estimates from step 2. (In my example, the simple average of the cohort differences in treatment effects by sex.)

Step 3 fails because the second call to hypotheses expects a hypothesis with length equal to the number of terms returned by slopes in Step 1 (4 in my example). This is unexpected because Step 3 is being called with the object from Step 2, which has fewer dimensions (it is length 2).

Example

library(dplyr)
library(marginaleffects)
library(tibble)
set.seed(5431)

test_data <- tibble(
  id = 1:200,
  male = rep(0:1, each = 100),
  treat = rep(0:1, times = 100),
  cohort = rep(c(1, 2, 2, 1), times = 50),
  y = treat * cohort + treat * male + rnorm(200)
)

test_data$male <- factor(test_data$male, levels = 0:1)
test_data$cohort <- factor(test_data$cohort, levels = 1:2)
m <- lm(y ~ treat : cohort : male + male, data = test_data)

effects <- slopes(
  m,
  newdata = datagrid(treat = 0:1, male = 0:1, cohort = 1:2),
  variables = "treat",
  by = c("male", "cohort")
)

# how treatment effect differs for male by cohort
hyp <- rbind(diag(-1, 2), diag(1, 2))
colnames(hyp) <- c("cohort1", "cohort2")

male_diff <- hypotheses(effects, hyp = hyp)
male_diff
#> 
#>     Term Estimate Std. Error    z Pr(>|z|)    S 2.5 % 97.5 %
#>  cohort1     1.09      0.356 3.05  0.00229  8.8 0.388   1.79
#>  cohort2     1.29      0.356 3.61  < 0.001 11.7 0.589   1.99
#> 
#> Columns: term, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high 
#> Type:  response

# want to calculate an average of the above cohort effects...

# ========== ERROR: This fails with an error about hypothesis length ==========
hypotheses(male_diff, hypothesis = c(0.5, 0.5))
#> Error in get_hypothesis(out, hypothesis, by = by, newdata = original, :
#>   Assertion on 'hypothesis' failed: Must have length 4, but has length 2.

# this works and doesn't complain about length
hypotheses(as_tibble(male_diff), vcov = vcov(male_diff), hypothesis = c(0.5, 0.5))
#> 
#>    Term Estimate Std. Error    z Pr(>|z|)    S 2.5 % 97.5 %
#>  custom     1.19      0.291 4.08   <0.001 14.4 0.617   1.76
#> 
#> Columns: term, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high

Created on 2024-04-30 with reprex v2.0.2

Session info

R version 4.2.1 (2022-06-23)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS 14.4.1

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] grDevices graphics  utils     stats     methods   base     

other attached packages:
[1] tibble_3.2.1             marginaleffects_0.19.0.2 dplyr_1.1.4              magrittr_2.0.3          

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.12       fansi_1.0.4       utf8_1.2.3        R6_2.5.1          backports_1.4.1   lifecycle_1.0.3   jsonlite_1.8.4    pillar_1.9.0     
 [9] rlang_1.1.3       cli_3.6.1         data.table_1.15.4 checkmate_2.3.1   vctrs_0.6.5       generics_0.1.3    glue_1.6.2        compiler_4.2.1   
[17] pkgconfig_2.0.3   insight_0.19.10   tidyselect_1.2.0   

Further Comments

I realize I can get around this in various ways, but I found the behavior unexpected since I'm passing an object with two terms to the second call to hypotheses. The above example is simple, but my true use case is estimating expanded TWFE difference-in-differences models with cohort/time heterogeneity (Wooldridge 2021) in which there are a lot of coefficients. I then want to use a bunch of sequential calls to hypotheses to generate several quantities of interest on the way to calculating an ATT. For example, going from model estimates to cohort-by-covariate-specific treatment effects, and then using the cohort-specific effects in another call to hypotheses to look for differences in effects within cohort, and then another call to average across cohorts and eventually arrive at a single estimate of the ATT. I would love to accomplish this with a simple series of calls to hypotheses:

est <- hypotheses(my_slopes, ...)
est <- hypotheses(est, ...)
est <- hypotheses(est, ...)
vincentarelbundock commented 2 months ago

Ooof. I will look into this when I find some time, but here's some background for now.

Some days I regret introducing hypotheses(), and I would strongly discourage calling it sequentially like you are trying to do. That involves many lots of catching and modifying calls, quosures, and esoteric things. Doing this in a "safe" way would require a lot of behind the scenes tracking and a ton of code complexity.

The hypotheses() function will always work in the future, but please be aware that I may eventually raise an error when users try to call it twice on an object.

In my own personal work, always try to define my hypotheses directly in the hypothesis argument, rather than as a second call. I recommend taking a look at the specification of complex hypotheses using functions here (requires the dev version of marginaleffects:

https://marginaleffects.com/vignettes/hypothesis.html#functions

BTW, have you looked at the etwfe package? It uses marginaleffects under the hood to implement the Wooldridge (2021) strategy:

https://grantmcdermott.com/etwfe/

poliquin commented 2 months ago

Thanks for the quick reply! Also, I wasn't aware of etwfe, so I'm glad I mentioned what I was doing. The package implements exactly the test I need. And I'll try to avoid multiple calls to hypotheses...

vincentarelbundock commented 2 months ago

Great!

FYI, the new version on Github now raises an error when a user tries to call hypotheses() twice.