2
$\begingroup$

Background

Looking into the MatchIt articles made me realize that using Bootstrap with BCa is a better practice for assessing uncertainty estimation (since: "For nonlinear models (e.g., logistic regression), the delta method is only an approximation subject to error"). Using Estimating Effects After Matching article for clustered bootstrap, I got satisfactory results. Still, I would like to utilize the full functionality of the marginaleffects package, for instance - using the "by =" to get stratified results, etc.

The marginaleffects documentation shows how can marginaleffects::avg_comparisons can be used in combination with the boot package using marginaleffect::inferences.

The problem

When I try to use bootstrap after running marginaleffects::avg_comparisons with clusters G-comp according to MatchIt article, I get this error:

avg_comparisons(fit, variables = "A", vcov = ~subclass, newdata = subset(md, A == "1"), wts = "weights", comparison = "lnratioavg", transform = exp ) %>% inferences(method = "boot", R = 999, conf_type = "bca") # Error: Assertion failed. One of the following must apply: # * checkmate::check_choice(x): Must be element of set {'exp','ln'}, but is not atomic scalar # * checkmate::check_function(x): Must be a function, not 'list'

Also, running the code after removing the transform = exp or/and comparison = "lnratioavg" arguments shows results that are very different than the estimate and CI computed using delta method or bootstrapping after manually computing G-comp. as shown in the article above (which are very similar in my instance).

To avoid from posting very long code, I'm attaching a reproduced example with downloadable code

My questions:

  1. Is using marginaleffects with boot considered good practice in a matched cohort, and is it an appropriate and equal substation for the manual g-computation in the MatchIt articles?
  2. What is the proper code for doing it?
  3. How to run bootstrap with by = argument, and is this a good practice?

================

UPDATE - the reproduced example above has been updated with @Noah's solution.

Improve this question
$\endgroup$

1 Answer 1

Reset to default
2
$\begingroup$

Supplying vcov = ~subclass does not tell inferences() that you want clustering; it tells avg_comparisons() that you want to use a cluster-robust delta method-based standard error, which you don't want. You want the cluster bootstrap, in which clusters (matched pairs or strata) are resampled in each bootstrap sample.

marginaleffects does not support the cluster bootstrap when using inferences() with method = "boot" or method = "rsample" because the arguments passed to inferences() are passed directly to the corresponding bootstrapping functions in those packages, which also do not support the cluster bootstrap by default. method = "fwb" does support the cluster bootstrap through its cluster argument, but that uses the fractional weighted bootstrap, which has not been investigated in this scenario and isn't guaranteed to work well with weighted regression models supplied to avg_comparisons() (for example, the wts argument can't be used).

So, your options are the following:

  • Petition the boot or rsample authors to allow the cluster bootstrap to be used by supplying an argument to the corresponding function. This will likely not go well if these packages were not designed with this in mind.
  • Investigate whether fwb will work for your purposes by doing a manual assessment outside of inferences(). That is, make sure the fractional weighted cluster bootstrap produces results that are valid for use after matching, since this has not been investigated.
  • Manually code the cluster bootstrap as explained in the MatchIt vignette. Instead of computing all the quantities internally manually, you can instead just call avg_comparisons() within each bootstrap iteration and extract the effect estimates to be used in computing the confidence intervals. Remember to set vcov = FALSE in the call to avg_comparisons() since you will only be extracting the effect estimates, not the standard errors. This will run much more slowly than coding it yourself, but gives you full access to the infrastructure of avg_comparisons() up to the point the effects are estimated. An example demonstrating how to do this is in the WeightIt documentation.

Good on you for taking the time to figure out the best method instead of taking he easy way out and just using the cluster-robust SE.

Here is some example code, building on the liked file in the OP.

#Unique pair IDs pair_ids <- levels(md$subclass) #Unit IDs, split by pair membership split_inds <- split(seq_len(nrow(md)), md$subclass) cluster_boot_fun <- function(pairs, i) { #Extract units corresponding to selected pairs ids <- unlist(split_inds[pairs[i]]) #Subset md with block bootstrapped indices boot_md <- md[ids,] #Fit outcome model fit <- glm(Y_B ~ A * (X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8 + X9), data = boot_md, weights = weights, family = quasibinomial()) ## G-computation ## est <- avg_comparisons(fit, variables = "A", vcov = FALSE, newdata = subset(boot_md, A == "1"), wts = "weights", comparison = "ratioavg", by = "X_ORDINAL") setNames(est$estimate, paste(est$contrast, est$X_ORDINAL, sep = "|")) } cluster_boot_out <- boot(pair_ids, cluster_boot_fun, R = R) cluster_boot_out
Improve this answer
$\endgroup$
5
  • $\begingroup$ I used the last option, and the Rmarkdown page is now updated with the solution. Also, I used the newdata argument in the marginaleffects::avg_comparisons function to manually run a stratified bootstrap. A comparison to with calculation using the delta method looks promising. Can you please approve this approach for stratifying? I must say, I admire your work and comprehensive answers, thank you! $\endgroup$ Commented Jul 1, 2023 at 14:25
  • 1
    $\begingroup$ Your code is not correct; I'll provide an example in my answer. $\endgroup$ Commented Jul 1, 2023 at 16:26
  • $\begingroup$ Thank you, the code in the Rmarkdown is updated accordingly. $\endgroup$ Commented Jul 1, 2023 at 17:14
  • $\begingroup$ Followup problem: In rare instances, we get a bootstrap data copy that doesn't contain all X_ORDINAL levels. The result is, for example, that the G-comp gives only 3 rows of est instead of 4 as a response to the cluster_boot_fun, and then the boot function fails with an error (boot::boot line 124: ...t.star[r, ] <- res[[r]] gives: number of items to replace is not a multiple of replacement length). Can I avoid this? I already tried to expand the factor beforehand, doesn't help. $\endgroup$ Commented Dec 18, 2023 at 14:30
  • 1
    $\begingroup$ This is a fundamental problem with the basic bootstrap and a reason to use consider an alternative method. I think at this point just using the cluster-robust SE should be fine. $\endgroup$ Commented Dec 18, 2023 at 16:30

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.