This blogpost documents the challenges we faced when attempting to run a mediation analysis using the mediation package with a binomial mediator and outcome while incorporating sampling weights into the analysis with the help of the survey package. Issues and solutions are summarized below. Click the issue number to jump straight to the relevant section of the post.

Summary: Issues and Solutions

Issue 1: Correctly applying the Fragile Families weights

Description: We need to specify jackknife estimation of standard errors for survey to properly handle our replicate weights.

Solution: Besides specifying the weights and repweights arguments, we also have to specify combined_weights = TRUE, type = "JKn", scales = 1, rscales = 1, and mse = TRUE.

Issue 2: Mediation with “multiple-trial binomial mediator” runs for weighted but not for unweighted models?

Description: When using multiple-trial binomial GLMs for the mediator and outcome models, mediate() complains: “weights on outcome and mediator models not identical”.

Solution: The problem occurs when we specify the GLMs with a proportion as the outcome and the total number of trials given in the weights argument, if the number of trials is different for both models. We can overcome this by specifying the outcome as cbind(successes, failures). However, according to the package documentation, mediate() should not accept multiple-trial binomial mediators in the first place, leaving us with some doubts about whether all the estimates are calculated correctly in this case.

Issue 3: Using a binary mediator and a binomial outcome — works with cbind(successes, failures) but not with proportion as outcome and weights argument specified

Description: When using a single trial binomial GLM for the mediator model and a multiple-trial binomial GLM for the outcome model, mediate() complains: “weights on outcome and mediator models not identical”.

Solution: As in Issue 2, the issue is solved by using the cbind(successes, failures) specification for the outcome model.

Issue 4: Ordinal mediator and binomial outcome don’t work together

Description: When the mediator model is ordinal (fit with MASS::polr(), mediate() will not work with a binomial outcome model either way it’s specified.

Solution: No solution was found, other than using an ordinal mediator model with a binary outcome model.

Issue 5: mediate() is not compatible with svyolr() models

Description: It seems that mediate() simply won’t work with models fit with survey::svyolr().

Solution: No solution was found.

Issue 6: Options and choices to boot — To robustSE or not to robustSE?

Description: mediate() can work via boostrap or quasi-Bayesian approximation, and has an option to calculate “robust standard errors” for the latter. Which options work with survey-weighted GLMs?

Solution: Bootstrapping should probably not be used. robustSE should be set to FALSE. If marginaleffects is used for comparisons, vcov should be set to vcov(model) or NULL, as the sandwich package does not correctly compute robust standard errors for models fit with survey.

Issue 7: Bounds of the confidence intervals are rounded inconsistently

Description: When you apply summary() to an object produced by mediate(), the upper bound of the confidence intervals is rounded to only two decimal places, whereas there are more decimal places for the lower bound.

Solution: We can modify the function print.summary.mediate() to get around this problem, although the current solution then messes with how the p-values are reported.

Background

Code

# load the packages we will use
library(tidyverse)
library(here)
library(mediation)
library(survey)
library(srvyr)
library(gt)
library(MASS)
library(marginaleffects)


# directory for the post
post_folder <- "2023-12-21_mediation-weighted-glm"

# number of simulations for mediation analysis
n_sims <- 1000

As part of a research program seeking to better understand intimate partner violence (IPV) through an intersectional lens, we wanted to use the Future of Families & Child Wellbeing Study ¹ dataset to run a mediation analysis. The data come from a survey with a complex design, and baseline and replicate weights are provided for each participant in the “national” sample.² These weights should allow us to make estimates nationally representative. Our mediator and outcome variables could reasonably be construed as binomial outcomes, and this was our preferred approach to fitting the necessary models. However, we encountered several difficulties along the way, and attempt to document and make sense of them here.

For demonstrative purposes we created a synthetic dataset using the synthpop package. This dataset mirrors the characteristics of our data without it being possible to identify any real individuals from it.

Code

# get synthetic dataset
dat <- readRDS(here::here("posts", post_folder, "data_reduced_synth.RDS")) |> 
  # create other variables we will need
  dplyr::mutate(
    informal_support_prop = informal_support/3,
    informal_support_max = 3,
    informal_support_binary = ifelse(informal_support == 3, 1, 0),
    informal_support_ordinal = forcats::as_factor(informal_support) |> 
      forcats::fct_relevel("0", "1", "2", "3"),
    IPV_prop = IPV/24,
    IPV_max = 25,
    IPV_binary = ifelse(IPV > 0, 1, 0),
    IPV_ordinal = forcats::as_factor(IPV) |> 
      forcats::fct_expand("0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19", "20", "21", "22", "23", "24") |> 
      forcats::fct_relevel("0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19", "20", "21", "22", "23", "24")
  ) |> 
  # move all the non-weight variables to the beginning
  dplyr::relocate(
    !dplyr::starts_with("m1"),
    .before = "m1natwt"
  )

We have three variables that will be used for mediation, race (“White” or “Black”), informal_support (an integer from 0 to 3), and IPV (an integer from 0 to 24). We manually create _prop and _max variables for informal support and IPV, so that we can use the proportion as an outcome in our GLMs and use the weights argument to specify the number of trials. We also create informal_support_binary, which classifies women as having high access to informal support (takes the value 1) when they score 3 out of 3, and low access (takes the value 0) when they score less than 3, as well as IPV_binary, which classifies women as having experienced some form of IPV (takes the value 1) if they score more than 0 and as not having experienced IPV (takes the value 0) if they score exactly 0. These binarized variables will become useful later on. Finally, we create the ordered factors informal_support_ordinal and IPV_ordinal to be used in ordinal models later on.

Each person also has a baseline weight, m1natwt, and 33 replicate weights, m1natwt_rep1 and so on. Behold:

Code

head(dat) |> 
  gt::gt()

race	informal_support	IPV	informal_support_prop	informal_support_max	informal_support_binary	informal_support_ordinal	IPV_prop	IPV_max	IPV_binary	IPV_ordinal	m1natwt	m1natwt_rep1	m1natwt_rep2	m1natwt_rep3	m1natwt_rep4	m1natwt_rep5	m1natwt_rep6	m1natwt_rep7	m1natwt_rep8	m1natwt_rep9	m1natwt_rep10	m1natwt_rep11	m1natwt_rep12	m1natwt_rep13	m1natwt_rep14	m1natwt_rep15	m1natwt_rep16	m1natwt_rep17	m1natwt_rep18	m1natwt_rep19	m1natwt_rep20	m1natwt_rep21	m1natwt_rep22	m1natwt_rep23	m1natwt_rep24	m1natwt_rep25	m1natwt_rep26	m1natwt_rep27	m1natwt_rep28	m1natwt_rep29	m1natwt_rep30	m1natwt_rep31	m1natwt_rep32	m1natwt_rep33
Black	3	18	1	3	1	3	0.75000000	25	1	18	41.13614	37.85516	41.24734	42.44724	44.39487	43.73881	44.45506	40.57142	40.78370	41.43095	34.34425	44.46661	0.00000	0.00000	37.94593	41.50523	43.69954	41.31135	39.22795	42.50257	42.02170	46.58911	41.82341	45.53830	42.27204	42.28968	43.25393	40.67841	0.00000	42.64928	42.28931	42.82657	43.13633	44.10038
Black	3	4	1	3	1	3	0.16666667	25	1	4	19.97570	22.96839	20.49259	20.21818	21.76786	20.79500	20.41886	21.54646	19.36423	19.80069	19.47093	19.58197	20.42065	19.78305	20.23575	21.36903	20.48008	20.91415	20.30889	21.36875	19.87386	21.06265	20.41880	23.67257	20.81135	21.42606	21.32215	19.58435	20.62730	21.86095	20.91533	20.08396	20.84812	20.38860
Black	3	0	1	3	1	3	0.00000000	25	0	0	12.24003	12.68406	11.59396	12.10748	12.00069	13.18730	11.88648	12.80008	0.00000	12.55978	11.89597	12.56892	12.57526	12.47836	11.46254	12.59188	12.65551	11.98436	12.49196	12.09257	13.70510	13.18920	12.23701	12.59980	12.74326	13.78494	12.39277	13.22096	11.71815	13.12761	12.33747	11.22335	12.79405	13.22079
Black	3	0	1	3	1	3	0.00000000	25	0	0	121.76894	129.03618	131.50407	124.49722	125.08351	115.39091	124.86470	0.00000	131.66747	121.04597	126.42019	125.73327	120.78219	132.49487	114.30053	127.66466	124.08061	125.39723	130.57684	134.18793	122.83195	0.00000	124.68690	119.24700	128.60410	128.66008	124.02149	126.49405	123.43700	141.58299	135.33669	126.14429	131.47270	125.73920
White	3	11	1	3	1	3	0.45833333	25	1	11	85.96230	87.22457	83.79171	91.27752	80.11185	88.11595	88.67043	90.75494	91.45559	85.68980	86.45096	92.47426	81.41688	99.84825	79.96129	85.80882	87.77487	89.85811	84.30115	95.09169	89.03872	90.58102	90.92682	102.89832	93.30130	94.10433	86.57666	88.45363	82.06447	89.72732	97.44034	91.89236	94.17548	0.00000
Black	3	2	1	3	1	3	0.08333333	25	1	2	98.04827	101.53748	0.00000	99.79310	97.24787	102.89133	105.71425	105.91802	102.19711	100.52937	103.41740	95.16078	97.53283	98.28993	103.28411	105.16373	102.85249	97.12798	107.91775	102.52666	107.20479	100.78193	95.96693	97.10217	101.69630	99.66428	101.72647	96.58415	95.98499	100.43622	0.00000	105.47424	96.22710	95.11009

We’ll need to run two models — the first predicting the mediator, the second predicting the outcome — and feed these models into mediation::mediate().

Issue 1: Correctly applying the Fragile Families weights

We fit the models using survey::svyglm() in order to be able to take the weights into account. However, svyglm() will not work with our dataframe directly — we need to create a “survey design” object by telling survey specifically how to handle the weights in our data. To do this correctly, we refer to the FFCW’s brief guide as well as other resources, like this CrossValidated question and the answers to it. Let’s build up the code step by step.

First, we pick the function that will create our survey design object. As we are dealing we replicate weights, we use the survey::svrepdesign() function and look at its documentation. We ignore the variables argument, as we will want all the variables that are already in our data to be included.

Next, we see the repweights and weights arguments, which ask us to specify the replicate weights and the sampling weights, respectively. Our FFCW guide says (emphasis in the original):

When setting the weights in your statistical package, you will need to indicate the name of the basic weight, as well as the replicate weights.

…

As an example, to weight the baseline mother data to be representative of births in the 77 cities, you would apply the weight m1natwt (and replicate weights m1natwt_rep1- m1natwt_rep33).

To do this in survey::svrepdesign(), we specify weights = ~m1natwt and repweights = "m1natwt.rep[1-9]+".

The data argument is straightforward, we just need to specify the dataframe we wish to use, dat.

Now the fun begins. How do we set the type argument, and do we need anything else? The FFCW guide is made for Stata users, and shows the following general command for specifying the survey design:

svyset [pweight=BASICWEIGHT], jkrw(REPLICATES, multiplier(1)) vce(jack) mse

So far we’ve incorporated the “pweight=BASICWEIGHT” and the “REPLICATES” part, but the Stata code reveals some more intricacies. We can see based on “vce(jack)” and “jkrw” that Jackknife variance estimation is used, and that corresponds to setting type = "JKn".

Next we set combined.weights = TRUE, which does not seem to have an equivalent in the Stata code, but seems to be correct approach according to the function documentation.³

Stata’s “multiplier(1)” means we need to set scales = 1 as well as rscales = 1, as explained by Lumley in this answer.

And finally, “mse” has its survey equivalent in the argument mse = TRUE.

Putting it all together yields:

survey::svrepdesign(
  data = dat,
  weights = ~m1natwt,
  repweights = "m1natwt.rep[1-9]+",
  combined.weights = FALSE,
  type = "JKn",
  scales = 1,
  rscales = 1,
  mse = TRUE
)

Because we like to stay close to the tidyverse, the code we actually run uses the srvyr package instead. This way we get to avoid specifying our sampling weight as a formula (recall “~m1natwt”) and our replicate weights as a regular expression:

Code

dat_weights <- dat |> 
  srvyr::as_survey_rep(
    repweights = dplyr::contains("_rep"),
    weights = m1natwt,
    combined_weights = FALSE,
    type = "JKn",
    scales = 1,
    rscales = 1,
    mse = TRUE
)

Running the weighted models and mediation

With the weighted dataframe in hand, we can use survey::svyglm() to fit our GLMs predicting informal_support and IPV.

Our informal support variable is the sum of three items for which each mother could get a 1 or a 0 (indicating that she did or did not have access to that form of support). As such, it seems sensible to model it as a binomial outcome — each mother has a certain number of “successes” out of three possible trials.

Code

model_m_1 <- survey::svyglm(
  formula = informal_support_prop ~ race,
  weights = informal_support_max,
  design = dat_weights,
  family = "binomial"
)

summary(model_m_1)


Call:
survey::svyglm(formula = informal_support_prop ~ race, weights = informal_support_max, 
    design = dat_weights, family = "binomial")

Survey design:
Called via srvyr

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)    3.689      7.167   0.515    0.610
raceBlack     -1.732      7.069  -0.245    0.808

(Dispersion parameter for binomial family taken to be 1299.299)

Number of Fisher Scoring iterations: 6

We get the following warning many times:

Warning: non-integer #successes in a binomial glm!

This is survey’s way of asking us to specify the family as “quasibinomial” instead of “binomial”, since our outcome is a proportion (weighted by the number of trials) and not binary. However, if we do that, mediation will not work with our models.

Our IPV variable is the sum of twelve items for which each mother could get a score of 0, 1, or 2 (indicating how often she experienced a given form of violence from her partner). As such, it also seems sensible to model it as a binomial outcome — each mother has a certain number of “successes” out of 24 possible “trials”.⁴

Code

model_y_1 <- survey::svyglm(
  formula = IPV_prop ~ race + informal_support_prop,
  weights = IPV_max,
  design = dat_weights,
  family = "binomial"
)

summary(model_y_1)


Call:
survey::svyglm(formula = IPV_prop ~ race + informal_support_prop, 
    weights = IPV_max, design = dat_weights, family = "binomial")

Survey design:
Called via srvyr

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)
(Intercept)            -1.4968     1.3423  -1.115    0.274
raceBlack               0.1207     1.2946   0.093    0.926
informal_support_prop  -0.7482     1.3249  -0.565    0.576

(Dispersion parameter for binomial family taken to be 291.4434)

Number of Fisher Scoring iterations: 5

(Here again we suppressed the warning that would have otherwise popped up.)

Finally, our mediation:

Code

med_1 <- mediation::mediate(
  model.m = model_m_1, # binomial svyglm, with weights argument
  model.y = model_y_1, # binomial svyglm, with weights argument
  sims = n_sims,
  treat = "race",
  mediator = "informal_support_prop"
)

weights taken as sampling weights, not total number of trials

Warning in mediation::mediate(model.m = model_m_1, model.y = model_y_1, :
treatment and control values do not match factor levels; using White and Black
as control and treatment, respectively

Code

summary(med_1)


Causal Mediation Analysis 

Quasi-Bayesian Confidence Intervals

                         Estimate 95% CI Lower 95% CI Upper p-value
ACME (control)           -0.03074     -0.40904         0.08    0.93
ACME (treated)           -0.00935     -0.17430         0.15    0.93
ADE (control)             0.02669     -0.37848         0.41    0.91
ADE (treated)             0.04808     -0.15878         0.45    0.91
Total Effect              0.01734     -0.52797         0.47    0.94
Prop. Mediated (control)  0.02231     -3.41339         2.60    0.92
Prop. Mediated (treated)  0.03808     -2.52514         2.72    0.92
ACME (average)           -0.02005     -0.27917         0.10    0.93
ADE (average)             0.03738     -0.25409         0.41    0.91
Prop. Mediated (average)  0.03019     -2.97625         2.68    0.92

Sample Size Used: 1955 


Simulations: 1000

A small but statistically significant mediation effect! The warning about weights not being taken as total number of trials makes us nervous, but then again, we did give the models sampling weights, so maybe it’s ok?

Issue 2: Mediation with “multiple-trial binomial mediator” runs for weighted but not for unweighted models?

We knew that we would eventually perform a sensitivity analysis where all models would be run without using the sampling weights. Running the previous steps with the glm() command resulted in an issue:

Code

model_m_2 <- glm(
  formula = informal_support_prop ~ race,
  weights = informal_support_max,
  data = dat,
  family = "binomial"
)

model_y_2 <- glm(
  formula = IPV_prop ~ race + informal_support_prop,
  weights = IPV_max,
  data = dat,
  family = "binomial"
)

med_2 <- mediation::mediate(
  model.m = model_m_2, # binomial glm, with weights argument
  model.y = model_y_2, # binomial glm, with weights argument
  sims = n_sims,
  treat = "race",
  mediator = "informal_support_prop"
)

Error in mediation::mediate(model.m = model_m_2, model.y = model_y_2, : weights on outcome and mediator models not identical

It seems like mediation::mediate() is taking the weights we pass to the GLM function (which correspond to total number of trials, not to sampling weights) as sampling weights, and complaining when it notices that these weights are not the same for the mediator and outcome models. Vasco perhaps clumsily posted this as an issue on the mediation package Github, but at least for now there has been no answer. We can work around this by using an alternative way of specifying the same GLM model — instead of using a proportion as an outcome and passing the number of trials through the weights argument, we can instead use a vector of “successes” and “failures” (which corresponds to the total number of trials minus the number of successes) as the outcome:

Code

model_m_2.1 <- glm(
  formula = cbind(informal_support, informal_support_max - informal_support) ~ race,
  data = dat,
  family = "binomial"
)

model_y_2.1 <- glm(
  formula = cbind(IPV, IPV_max - IPV) ~ race + informal_support_prop,
  data = dat,
  family = "binomial"
)

med_2.1 <- mediation::mediate(
  model.m = model_m_2.1, # binomial glm, with cbind() approach
  model.y = model_y_2.1, # binomial glm, with cbind() approach
  sims = n_sims,
  treat = "race",
  mediator = "informal_support_prop"
)

summary(med_2.1)


Causal Mediation Analysis 

Quasi-Bayesian Confidence Intervals

                         Estimate 95% CI Lower 95% CI Upper p-value    
ACME (control)            0.00294      0.00157         0.00  <2e-16 ***
ACME (treated)            0.00324      0.00172         0.00  <2e-16 ***
ADE (control)             0.01295      0.00763         0.02  <2e-16 ***
ADE (treated)             0.01325      0.00785         0.02  <2e-16 ***
Total Effect              0.01619      0.01086         0.02  <2e-16 ***
Prop. Mediated (control)  0.17910      0.09765         0.31  <2e-16 ***
Prop. Mediated (treated)  0.19892      0.11054         0.33  <2e-16 ***
ACME (average)            0.00309      0.00165         0.00  <2e-16 ***
ADE (average)             0.01310      0.00776         0.02  <2e-16 ***
Prop. Mediated (average)  0.18901      0.10409         0.32  <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sample Size Used: 1955 


Simulations: 1000

Success! But at the cost of suspicion. The documentation for the mediate()function warns us that (emphasis added):

As of version 4.0, the mediator model can be of either ‘lm’, ‘glm’ (or `bayesglm’), ‘polr’ (or `bayespolr’), ‘gam’, ‘rq’, `survreg’, or `merMod’ class, corresponding respectively to the linear regression models, generalized linear models, ordered response models, generalized additive models, quantile regression models, parametric duration models, or multilevel models.. For binary response models, the ‘mediator’ must be a numeric variable with values 0 or 1 as opposed to a factor. Quasi-likelihood-based inferences are not allowed for the mediator model because the functional form must be exactly specified for the estimation algorithm to work. The ‘binomial’ family can only be used for binary response mediators and cannot be used for multiple-trial responses. This is due to conflicts between how the latter type of models are implemented in glm and how ‘mediate’ is currently written.

Thus, until the package author clarifies whether this implementation results in correct estimates, we are cautious about trusting them.

Sanity check: comparing results with linear probability models

As a sanity check, to become a little more confident that mediate is working as intended even though it’s not supposed to work well with a multiple-trial binomial mediator, we compare the results of our model that did run (med_2.1) with a mediation run on linear probability models (basically, using lm() instead of glm() even though the outcomes are proportions.

Code

model_m_2.2 <- lm(
  formula = informal_support_prop ~ race,
  data = dat
)

model_y_2.2 <- lm(
  formula = IPV_prop ~ race + informal_support_prop,
  data = dat
)

med_2.2 <- mediation::mediate(
  model.m = model_m_2.2, # lm
  model.y = model_y_2.2, # lm
  sims = n_sims,
  treat = "race",
  mediator = "informal_support_prop"
)

summary(med_2.2)


Causal Mediation Analysis 

Quasi-Bayesian Confidence Intervals

               Estimate 95% CI Lower 95% CI Upper p-value    
ACME            0.00314      0.00146         0.01  <2e-16 ***
ADE             0.01372      0.00110         0.03   0.036 *  
Total Effect    0.01686      0.00401         0.03   0.008 ** 
Prop. Mediated  0.18518      0.07060         0.72   0.008 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sample Size Used: 1955 


Simulations: 1000

The results seem fairly similar, but this provides limited comfort.

Bonus finding: you can specify a mediator that is not technically the variable that your mediator model predicts

Note that in the models used for med_2.1 there is a potential discrepancy. The mediator model is formulated as formula = cbind(informal_support, informal_support_max - informal_support) ~ race, while the outcome model is formulated as formula = cbind(IPV, IPV_max - IPV) ~ race + informal_support_prop. So the outcome model is using informal_support_prop (our intended mediator) as a predictor of IPV, but the mediator model uses the variables informal_support and informal_support_max to internally generate the proportion of successes. In the call to mediation::mediate() we can then specify mediator = "informal_support_prop", and all seems to work as intended. What if we used informal_support as predictor in the outcome model instead? We can then specify mediator = "informal_support" and end up with something like this:

Code

model_m_2.3 <- glm(
  formula = cbind(informal_support, informal_support_max - informal_support) ~ race,
  data = dat,
  family = "binomial"
)

model_y_2.3 <- glm(
  formula = cbind(IPV, IPV_max - IPV) ~ race + informal_support,
  data = dat,
  family = "binomial"
)

med_2.3 <- mediation::mediate(
  model.m = model_m_2.3, # binomial glm, with cbind() approach
  model.y = model_y_2.3, # binomial glm, with cbind() approach
  sims = n_sims,
  treat = "race",
  mediator = "informal_support"
)

summary(med_2.3)


Causal Mediation Analysis 

Quasi-Bayesian Confidence Intervals

                         Estimate 95% CI Lower 95% CI Upper p-value    
ACME (control)           0.001150     0.000595         0.00  <2e-16 ***
ACME (treated)           0.001260     0.000660         0.00  <2e-16 ***
ADE (control)            0.017324     0.010286         0.02  <2e-16 ***
ADE (treated)            0.017434     0.010369         0.02  <2e-16 ***
Total Effect             0.018584     0.011436         0.03  <2e-16 ***
Prop. Mediated (control) 0.061515     0.031134         0.12  <2e-16 ***
Prop. Mediated (treated) 0.067573     0.034355         0.13  <2e-16 ***
ACME (average)           0.001205     0.000625         0.00  <2e-16 ***
ADE (average)            0.017379     0.010328         0.02  <2e-16 ***
Prop. Mediated (average) 0.064544     0.032554         0.12  <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sample Size Used: 1955 


Simulations: 1000

Looking at the last three rows (Average Causal Mediated Effect (average), Average Direct Effect (average), and Proportion Mediated (average)), the results seem different than the ones estimated in med_2.1 (which used informal_support_prop as mediator) and med.2.2 (which used linear proability models). But looking closer, you’ll notice that the estimates related to the mediator, namely of the Average Causal Mediated Effect and the Proportion Mediated just seem to be divided by three!

My assumption about how that came to be: When we use informal_support_prop, which can only take values on the \((0, 1)\) interval, the ACME estimates correspond to the change in probability of IPV associated with going from 0 to 1 on the mediator and thus span the full range of the mediator. When we use informal_support, which takes values in the set \(\{0, 1, 2, 3\}\), the estimates also correspond to the increase in probability of IPV associated with an increase of 1 on the mediator (e.g., from 0 to 1), but this is only \(\frac{1}{3}\) of the total possible way. Thus, applying this effect three times would be equivalent to moving from 0 to 3 on the mediator scale, the same as moving from 0 to 1 when the mediator is a proportion.⁵

Issue 3: Using a binary mediator and a binomial outcome — works with cbind(successes, failures) but not with proportion as outcome and weights argument specified

Next, we attempted to run everything with the same binomial outcome and a binary mediator.

Code

model_m_3 <- glm(
  formula = informal_support_binary ~ race,
  data = dat,
  family = "binomial"
)

model_y_3 <- glm(
  formula = IPV_prop ~ race + informal_support_binary,
  weights = IPV_max,
  data = dat,
  family = "binomial"
)

med_3 <- mediation::mediate(
  model.m = model_m_3, # binary glm
  model.y = model_y_3, # binomial glm, with weights argument
  sims = n_sims,
  treat = "race",
  mediator = "informal_support_binary"
)

Error in mediation::mediate(model.m = model_m_3, model.y = model_y_3, : weights on outcome and mediator models not identical

Once again we find that mediate() complains about different weights on mediator and outcome models, even though, according to the documentation, it works with a binary mediator and a binomial outcome. We then try using the vector of successes and failures specification for the outcome model:

Code

model_y_3.1 <- glm(
  formula = cbind(IPV, IPV_max - IPV) ~ race + informal_support_binary,
  data = dat,
  family = "binomial"
)

med_3.1 <- mediation::mediate(
  model.m = model_m_3,   # binary glm 
  model.y = model_y_3.1, # binomial glm, with cbind() approach
  sims = n_sims,
  treat = "race",
  mediator = "informal_support_binary"
)

summary(med_3.1)


Causal Mediation Analysis 

Quasi-Bayesian Confidence Intervals

                         Estimate 95% CI Lower 95% CI Upper p-value    
ACME (control)            0.00344      0.00166         0.01   0.002 ** 
ACME (treated)            0.00379      0.00181         0.01   0.002 ** 
ADE (control)             0.01240      0.00714         0.02  <2e-16 ***
ADE (treated)             0.01275      0.00734         0.02  <2e-16 ***
Total Effect              0.01619      0.01055         0.02  <2e-16 ***
Prop. Mediated (control)  0.21219      0.10668         0.38   0.002 ** 
Prop. Mediated (treated)  0.23359      0.11912         0.40   0.002 ** 
ACME (average)            0.00361      0.00174         0.01   0.002 ** 
ADE (average)             0.01257      0.00724         0.02  <2e-16 ***
Prop. Mediated (average)  0.22289      0.11285         0.39   0.002 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sample Size Used: 1955 


Simulations: 1000

This approach runs, but we are still left wondering whether we can trust the estimates.

Issue 4: Ordinal mediator and binomial outcome don’t work together

Because we would prefer not to binarize the social support variable in order to keep all the information that we can, we also explored using an ordinal mediator model. mediate() works with ordinal models fit with MASS::polr(), so we attempted that.

Code

model_m_4 <- MASS::polr(
  formula = informal_support_ordinal ~ race,
  data = dat,
  Hess = TRUE
)

model_y_4 <- glm(
  formula = cbind(IPV, IPV_max - IPV) ~ race + informal_support_ordinal,
  data = dat,
  family = "binomial"
)

med_4 <- mediation::mediate(
  model.m = model_m_4, # ordinal polr
  model.y = model_y_4, # binomial glm, with cbind() approach
  sims = n_sims,
  treat = "race",
  mediator = "informal_support_ordinal"
)

Error in UseMethod("isSymmetric"): no applicable method for 'isSymmetric' applied to an object of class "c('double', 'numeric')"

The first issue arises, but it’s easy to fix. Again, from the documentation of the mediate function:

The quasi-Bayesian approximation (King et al. 2000) cannot be used if ‘model.m’ is of class ‘rq’ or ‘gam’, or if ‘model.y’ is of class ‘gam’, ‘polr’ or ‘bayespolr’. In these cases, either an error message is returned or use of the nonparametric bootstrap is forced. Users should note that use of the nonparametric bootstrap often requires significant computing time, especially when ‘sims’ is set to a large value.

We can get around this by setting the argument boot to TRUE.

Code

med_4 <- mediation::mediate(
  model.m = model_m_4, # ordinal polr
  model.y = model_y_4, # binomial glm, with cbind() approach
  sims = n_sims,
  treat = "race",
  mediator = "informal_support_ordinal",
  boot = TRUE
)

Error in cbind(IPV, IPV_max - IPV): object 'IPV' not found

Now we get a different problem. It seems the boot package is having trouble with the way we set up our outcome model — hopefully using the other implementation will take care of that.

Code

model_y_4.1 <- glm(
  formula = IPV_prop ~ race + informal_support_ordinal,
  weights = IPV_max,
  data = dat,
  family = "binomial"
)

med_4 <- mediation::mediate(
  model.m = model_m_4,   # ordinal polr
  model.y = model_y_4.1, # binomial glm, with weights argument
  sims = n_sims,
  treat = "race",
  mediator = "informal_support_ordinal",
  boot = TRUE
)

Error in mediation::mediate(model.m = model_m_4, model.y = model_y_4.1, : weights on outcome and mediator models not identical

Finally, we reach a dead-end. It seems once again the mediate() function is taking the prior weights we fed into our glm() as sampling weights and complaining that our ordinal model doesn’t have the same sampling weights. If we are willing to throw away some information, we can make it all work with a binary outcome, thusly:

Code

model_y_4.2 <- glm(
  formula = IPV_binary ~ race + informal_support_ordinal,
  data = dat,
  family = "binomial"
)

med_4 <- mediation::mediate(
  model.m = model_m_4,   # ordinal polr
  model.y = model_y_4.2, # binary glm
  sims = n_sims,
  treat = "race",
  mediator = "informal_support_ordinal",
  boot = TRUE
)

summary(med_4)


Causal Mediation Analysis 

Nonparametric Bootstrap Confidence Intervals with the Percentile Method

                          Estimate 95% CI Lower 95% CI Upper p-value
ACME (control)            0.003298    -0.000859         0.01    0.10
ACME (treated)            0.003204    -0.000878         0.01    0.10
ADE (control)             0.012612    -0.026731         0.05    0.54
ADE (treated)             0.012518    -0.026356         0.05    0.54
Total Effect              0.015816    -0.022595         0.06    0.43
Prop. Mediated (control)  0.208514    -2.101207         1.87    0.49
Prop. Mediated (treated)  0.202555    -2.134302         1.89    0.49
ACME (average)            0.003251    -0.000860         0.01    0.10
ADE (average)             0.012565    -0.026544         0.05    0.54
Prop. Mediated (average)  0.205535    -2.117755         1.88    0.49

Sample Size Used: 1955 


Simulations: 1000

This runs, although the bootstrap adds significantly to computing time and we had to binarize IPV to achieve it.

Issue 5: mediate() is not compatible with svyolr() models

With this partial success, we checked to see if we could also run this latest mediation on the weighted data.

Code

model_m_5 <- survey::svyolr(
  formula = informal_support_ordinal ~ race,
  design = dat_weights
)

Error in optim(s0, fmin, gmin, method = "BFGS", ...): initial value in 'vmmin' is not finite

This throws an error, which this Stack Overflow post can help us solve. We can help the function by providing plausible start values for the coefficients and the intercepts/thresholds. We have one coefficient and three thresholds (from 0 to 1, from 1 to 2, and from 2 to 3), and we know they should be relatively close to zero as well as that the thresholds should be in ascending order. With the syntax start = c(coefficients, thresholds) we can make it work:

Code

model_m_5 <- survey::svyolr(
  formula = informal_support_ordinal ~ race,
  design = dat_weights,
  start = c(0, -1, 0, 1)
)

summary(model_m_5)

Call:
svyolr(formula = informal_support_ordinal ~ race, design = dat_weights, 
    start = c(0, -1, 0, 1))

Coefficients:
             Value Std. Error    t value
raceBlack -1.62481   6.985481 -0.2325982

Intercepts:
    Value   Std. Error t value
0|1 -4.5463  7.1876    -0.6325
1|2 -3.7384  7.1644    -0.5218
2|3 -3.0236  7.0091    -0.4314

Ok, moving on.

Code

model_y_5 <- survey::svyglm(
  formula = IPV_binary ~ race + informal_support_ordinal,
  design = dat_weights,
  family = "binomial"
)

med_5 <- mediation::mediate(
  model.m = model_m_5, # ordinal svyolr
  model.y = model_y_5, # binomial svyglm, with weights argument
  sims = n_sims,
  treat = "race",
  mediator = "informal_support_ordinal",
  boot = TRUE
)

Error in mediation::mediate(model.m = model_m_5, model.y = model_y_5, : weights on outcome and mediator models not identical

Not this again 😭. We knew that the weights used for the svyglm() and svyolr() models were the same, so there must be something else happening in the background. Vasco did some investigating and again turned to Stack Overflow. We can see what weights are stored in the model objects generated by svyglm() and svyolr() using the model.frame() function. In the following two code blocks we also use dplyr::slice_head() to display just the first line of results output by model.frame().

Code

model.frame(model_m_5) |> dplyr::slice_head() |> gt::gt()

informal_support_ordinal	race	(weights)
3	Black	41.13614

Code

model.frame(model_y_5) |> dplyr::slice_head() |> gt::gt()

IPV_binary	race	informal_support_ordinal	(weights)
1	Black	3	0.1142592

The first weight for our ordinal model appears to be 41.13614, while the binomial model uses 0.1142592. We can extract all the weights for both models and compare them. Let’s look at their ratio.

Code

weight_comparison <- tibble::tibble(
  binomial_weights = model.frame(model_y_5)$`(weights)`,
  ordinal_weights = model.frame(model_m_5)$`(weights)`,
  weight_quotient = binomial_weights / ordinal_weights
)

weight_comparison |> dplyr::slice_head(n = 5) |> gt::gt()

binomial_weights	ordinal_weights	weight_quotient
0.11425923	41.13614	0.002777587
0.05548425	19.97570	0.002777587
0.03399776	12.24003	0.002777587
0.33822387	121.76894	0.002777587
0.23876779	85.96230	0.002777587

Would you look at that! The ratio seems to be constant. We can check the entire weight_quotient column just to be extra sure:

Code

weight_comparison |> dplyr::count(weight_quotient) |> gt::gt()

weight_quotient	n
0.002777587	218
0.002777587	1733
0.002777587	4

There seem to be three different values that all look the same — perhaps this will make more sense once we learn more about floating point arithmetic, but for now we’ll say the ratio between the weights is always the same. On Stack Overflow, Thomas Lumley, maintainer of the survey package, explains in a reply to my post:

Some functions in the survey package rescale the weights to have unit mean because the large weights (thousands or tens of thousands) in some national surveys can cause convergence problems.

Since the results aren’t affected at all, the mediate package could probably work around this fairly easily. However, there’s a rescale=FALSE option to svyglm to not rescale the weights, provided for this sort of purpose.

If you then have convergence problems in svyglm you could manually rescale the weights to have unit mean before doing any of the analyses.

Good to know!

Let’s set rescale = FALSE in the svyglm() model and try our luck:

Code

model_y_5.1 <- survey::svyglm(
  formula = IPV_binary ~ race + informal_support_ordinal,
  design = dat_weights,
  family = "binomial",
  rescale = FALSE
)

Error in is.data.frame(pwts): object 'pwts' not found

This is perplexing, but someone else on Stack Overflow a couple of years ago encountered a similar problem, where object “.survey.prob.weights” was not found when rescale was set to FALSE. After looking at the source code and the accepted answer by user xilliam to the Stack Overflow post, we try to assign the prior weights contained in the survey design object to an object named pwts in the global environment, and that makes svyglm() happy again.

Code

pwts <- dat_weights$pweights

model_y_5.1 <- survey::svyglm(
  formula = IPV_binary ~ race + informal_support_ordinal,
  design = dat_weights,
  family = "binomial",
  rescale = FALSE
)

# remove pwts from the global environment
rm(list = "pwts")

However, we can’t be sure if this correctly solves the problem or introduces new errors. Vasco posted about this issue on Stack Overflow in the hopes that perhaps Thomas Lumley can clarify this for us.

Now we try to run the mediation.

Code

med_5 <- mediation::mediate(
  model.m = model_m_5,   # ordinal svyolr
  model.y = model_y_5.1, # binary svyglm
  sims = n_sims,
  treat = "race",
  mediator = "informal_support_ordinal",
  boot = TRUE
)

Error in MASS::polr(formula, data = df, ..., Hess = TRUE, model = FALSE, : formal argument "data" matched by multiple actual arguments

Somehow it is not working again. What if the outcome model was ordinal as well, just for fun?

Code

model_y_5.2 <- survey::svyolr(
  formula = IPV_ordinal ~ race + informal_support_ordinal,
  design = dat_weights,
  start = c(0, 0, 0, 0, seq(from = -1, to = 1, length.out = 24))
)

med_5 <- mediation::mediate(
  model.m = model_m_5,   # ordinal svyolr
  model.y = model_y_5.2, # ordinal svyolr
  sims = n_sims,
  treat = "race",
  mediator = "informal_support_ordinal",
  boot = TRUE
)

Error in MASS::polr(formula, data = df, ..., Hess = TRUE, model = FALSE, : formal argument "data" matched by multiple actual arguments

Ok, we give up here.

Issue 6: Options and choices to boot — To robustSE or not to robustSE?

For now, we’ve ignored an optional setting to mediate(): the robustSE argument defaults to FALSE, and, according to the documentation, “If ‘TRUE’, heteroskedasticity-consistent standard errors will be used in quasi-Bayesian simulations. Ignored if ‘boot’ is ‘TRUE’ or neither ‘model.m’ nor ‘model.y’ has a method for vcovHC in the sandwich package.” So we have three options:

boot = TRUE and, automatically, robustSE = FALSE
boot = FALSE and robustSE = TRUE, and
boot = FALSE and robustSE = FALSE.

And we know that, when chosen, the heteroskedasticity-consistent standard errors will be calculated with the help of the sandwich package.

According to the mediation package vignettes, the authors recommend setting boot = TRUE: “In general, as long as computing power is not an issue, analysts should estimate confidence intervals via the bootstrap with more than 1000 resamples, which is the default number of simulations.” So we could do that. Unfortunately, computing power is an issue, as we wanted to run a multiverse analysis that meant we would use the mediate() function a bunch of times. To figure out whether to set robustSE to TRUE or FALSE, it could be useful to see which option gives use results that are closer to the bootstrap option. Let’s use a binary mediator and outcome, just so we are sure there aren’t other issues at play.

For convenience, we’ll create a function called compare_meds() that takes in different mediation models and outputs a table.

Code

get_meds <- function(med_model){
  
  boot <- as.character(med_model[["boot"]])
  
  robustSE <- as.character(med_model[["robustSE"]])
  
  results_te <- tibble::tibble(
    effect = "total_effect",
    estimate = med_model[["tau.coef"]],
    ci.low = med_model[["tau.ci"]][[1]],
    ci.high = med_model[["tau.ci"]][[2]],
    p.value = med_model[["tau.p"]]
  )
  
  results_de <- tibble::tibble(
    effect = "direct_effect",
    estimate = med_model[["z.avg"]],
    ci.low = med_model[["z.avg.ci"]][[1]],
    ci.high = med_model[["z.avg.ci"]][[2]],
    p.value = med_model[["z.avg.p"]]
  )
  
  results_ie <- tibble::tibble(
    effect = "indirect_effect",
    estimate = med_model[["d.avg"]],
    ci.low = med_model[["d.avg.ci"]][[1]],
    ci.high = med_model[["d.avg.ci"]][[2]],
    p.value = med_model[["d.avg.p"]]
  )
  
  results <- dplyr::bind_rows(results_te, results_de, results_ie)
  
  output <- results |> 
    dplyr::mutate(boot = boot, robustSE = robustSE) |> 
    dplyr::relocate(boot, robustSE)
  
  return(output)
}

compare_meds <- function(med_models, effect){
  if (class(med_models) != "list") {
    stop("Argument should be a list of models")
  }
  
  med_models_df <- med_models  |> 
    purrr::map(.f = ~ get_meds(.x)) |> 
    purrr::list_rbind() |> 
    dplyr::filter(effect == !!effect)
  
  return(med_models_df)
}

First, we compare unweighted models. To reduce simulation-dependent variability, we set sims = 3000.

Code

model_m_6 <- glm(
  formula = informal_support_binary ~ race,
  data = dat,
  family = "binomial"
)

model_y_6 <- glm(
  formula = IPV_binary ~ race + informal_support_binary,
  data = dat,
  family = "binomial"
)

med_6.1 <- mediation::mediate(
  model.m = model_m_6, # binary glm
  model.y = model_y_6, # binary glm
  sims = 3000,
  boot = TRUE,
  treat = "race",
  mediator = "informal_support_binary"
)

med_6.2 <- mediation::mediate(
  model.m = model_m_6, # binary glm
  model.y = model_y_6, # binary glm
  sims = 3000,
  robustSE = TRUE,
  treat = "race",
  mediator = "informal_support_binary"
)

med_6.3 <- mediation::mediate(
  model.m = model_m_6, # binary glm
  model.y = model_y_6, # binary glm
  sims = 3000,
  robustSE = FALSE,
  treat = "race",
  mediator = "informal_support_binary"
)

mod_list <- list(med_6.1, med_6.2, med_6.3)

compare_meds(mod_list, effect = "total_effect") |> gt::gt()

boot	robustSE	effect	estimate	ci.low	ci.high	p.value
TRUE	FALSE	total_effect	0.01503999	-0.02630048	0.05645876	0.4880000
FALSE	TRUE	total_effect	0.01479701	-0.02676844	0.05651676	0.4833333
FALSE	FALSE	total_effect	0.01458587	-0.02668057	0.05480012	0.4733333

Code

compare_meds(mod_list, effect = "direct_effect") |> gt::gt()

boot	robustSE	effect	estimate	ci.low	ci.high	p.value
TRUE	FALSE	direct_effect	0.01067146	-0.02968020	0.05288181	0.6073333
FALSE	TRUE	direct_effect	0.01072508	-0.03080966	0.05218602	0.6020000
FALSE	FALSE	direct_effect	0.01049483	-0.03054998	0.05206773	0.6073333

Code

compare_meds(mod_list, effect = "indirect_effect") |> gt::gt()

boot	robustSE	effect	estimate	ci.low	ci.high	p.value
TRUE	FALSE	indirect_effect	0.004368532	-1.682579e-04	0.009016640	0.05933333
FALSE	TRUE	indirect_effect	0.004071932	-7.255864e-05	0.009293128	0.05400000
FALSE	FALSE	indirect_effect	0.004091046	-1.968454e-04	0.009081814	0.06200000

We can see that the results are very similar, and in fact running the mediations several times reveals that the difference seems to be due to simulation variability. Given this, Vasco would personally rather specify robustSE = FALSE, since it does not make sense to correct for heteroskedasticity in a logistic regression model which cannot suffer from heteroskedasticity in the first place.

What about if we use weighted models?

Code

model_m_6.1 <- survey::svyglm(
  formula = informal_support_binary ~ race,
  design = dat_weights,
  family = "binomial"
)

model_y_6.1 <- survey::svyglm(
  formula = IPV_binary ~ race + informal_support_binary,
  design = dat_weights,
  family = "binomial"
)

med_6.4 <- mediation::mediate(
  model.m = model_m_6.1, # binary svyglm
  model.y = model_y_6.1, # binary svyglm
  sims = 3000,
  boot = TRUE,
  treat = "race",
  mediator = "informal_support_binary"
)

med_6.5 <- mediation::mediate(
  model.m = model_m_6.1, # binary svyglm
  model.y = model_y_6.1, # binary svyglm
  sims = 3000,
  robustSE = TRUE,
  treat = "race",
  mediator = "informal_support_binary"
)

med_6.6 <- mediation::mediate(
  model.m = model_m_6.1, # binary svyglm
  model.y = model_y_6.1, # binary svyglm
  sims = 3000,
  robustSE = FALSE,
  treat = "race",
  mediator = "informal_support_binary"
)

mod_list <- list(med_6.4, med_6.5, med_6.6)

compare_meds(mod_list, effect = "total_effect") |> gt::gt()

boot	robustSE	effect	estimate	ci.low	ci.high	p.value
TRUE	FALSE	total_effect	0.044763016	0.03783175	0.04736893	0.0000000
FALSE	TRUE	total_effect	-0.008917991	-0.23020520	0.17866819	0.9513333
FALSE	FALSE	total_effect	-0.032263015	-0.53714755	0.31789234	0.9733333

Code

compare_meds(mod_list, effect = "direct_effect") |> gt::gt()

boot	robustSE	effect	estimate	ci.low	ci.high	p.value
TRUE	FALSE	direct_effect	0.027850175	0.02773439	0.02806805	0.0000000
FALSE	TRUE	direct_effect	0.006836401	-0.23961806	0.15778901	0.8206667
FALSE	FALSE	direct_effect	-0.026254065	-0.51121904	0.29904897	0.9340000

Code

compare_meds(mod_list, effect = "indirect_effect") |> gt::gt()

boot	robustSE	effect	estimate	ci.low	ci.high	p.value
TRUE	FALSE	indirect_effect	0.01691284	0.009863525	0.01956228	0.0000000
FALSE	TRUE	indirect_effect	-0.01575439	-0.102836657	0.02885321	0.8973333
FALSE	FALSE	indirect_effect	-0.00600895	-0.173359841	0.18032563	0.9540000

While the estimates are similar, the confidence intervals are much tighter when boot or robustSE are set to TRUE, and suddenly all estimates are significant! What gives?

First, an answer by Thomas Lumley (maintainer of the survey package) on Stack Overflow suggests that bootstrapping doesn’t work well with these survey-weighted models, although his answer does not fully map to our situation. In another answer, this time about applying robust standard errors with the sandwich package⁶, Lumley is clearer: “For svyglm models, vcov() already produces the appropriate sandwich estimator, and I don’t think the”sandwich” package knows enough about the internals of the object to get it right.” A few years prior, Achim Zeileis, maintainer of the sandwich package, had given a different perspective which also boils down to do not use sandwich on models fit with survey.

There you have it, boot = FALSE and robustSE = FALSE is the way to go.

Bonus: this goes for using `marginaleffects` as well

Often we are not interested in mediation, but rather would like a clean and interpretable way to report results from our GLMs. One could not be faulted for preferring to do this on the predicted probabilities scale rather than dealing with odds-ratios and other headaches. For this, the marginaleffects package can be of enormous help.

For example, we can examine the predicted probability of experiencing IPV for Black and White women:

Code

marginaleffects::avg_predictions(
  model = model_y_6.1, # binary svyglm
  variables = "race"
) |> gt::gt()

race	estimate	std.error	statistic	p.value	conf.low	conf.high
White	0.6986166	0.09146286	7.638254	2.201870e-14	0.5193527	0.8778805
Black	0.7264507	0.30101376	2.413347	1.580675e-02	0.1364746	1.3164269

And then compare the probabilities themselves, for which we again have the option of “correcting” the standard errors via the vcov argument. Look at the difference in the confidence intervals and p-values spit out when we don’t or do specify vcov = "HC":

Code

marginaleffects::avg_comparisons(
  model = model_y_6.1, # binary svyglm
  variables = "race"
) |> gt::gt()

term	contrast	estimate	std.error	statistic	p.value	conf.low	conf.high
race	Black - White	0.02783414	0.2652147	0.1049494	0.9164159	-0.4919772	0.5476455

Code

marginaleffects::avg_comparisons(
  model = model_y_6.1, # binary svyglm
  variables = "race",
  vcov = "HC"
) |> gt::gt()

term	contrast	estimate	std.error	statistic	p.value	conf.low	conf.high
race	Black - White	0.02783414	0.1108251	0.2511537	0.8016953	-0.1893791	0.2450473

The estimates are the same, but the standard error (wrongly) becomes minuscule if we apply the “correction” for heteroskedasticity without realizing that sandwich does not work properly on models fit with survey.

Issue 7: Bounds of the confidence intervals are rounded inconsistently

Did you notice? Each time we use summary() on a mediate() object, the lower bound of the 95% CI is reported with 5 decimal places, while the upper bound is reported with 2 decimal places. This also bothered Stack Overflow user sateayam over 4 years ago, and luckily user Roland answered with a fix:

Code

trace(mediation:::print.summary.mediate, 
      at = 11,
      tracer = quote({
        printCoefmat <- function(x, digits) {
          p <- x[, 4]
          x[, 1:3] <- sprintf("%.6f", x[, 1:3])
          x[, 4] <- sprintf("%.2f", p)
          print(x, quote = FALSE, right = TRUE)
        } 
      }),
      print = FALSE)

[1] "print.summary.mediate"

Code

mediation:::print.summary.mediate(summary(med_6.6))


Causal Mediation Analysis 

Quasi-Bayesian Confidence Intervals

                          Estimate 95% CI Lower 95% CI Upper p-value
ACME (control)           -0.006650    -0.178268     0.208017    0.95
ACME (treated)           -0.005368    -0.187307     0.144949    0.95
ADE (control)            -0.026895    -0.518915     0.302295    0.93
ADE (treated)            -0.025613    -0.509389     0.299770    0.93
Total Effect             -0.032263    -0.537148     0.317892    0.97
Prop. Mediated (control)  0.079758    -4.342051     5.819054    0.75
Prop. Mediated (treated)  0.031101    -2.418510     3.520780    0.75
ACME (average)           -0.006009    -0.173360     0.180326    0.95
ADE (average)            -0.026254    -0.511219     0.299049    0.93
Prop. Mediated (average)  0.055430    -3.405108     4.643884    0.75

Sample Size Used: 1955 


Simulations: 3000

Code

untrace(mediation:::print.summary.mediate)

Unfortunately, this solution also changes the way p-values are presented. Look again at the results from med_6.5 and compare with the following:

Code

trace(mediation:::print.summary.mediate, 
      at = 11,
      tracer = quote({
        printCoefmat <- function(x, digits) {
          p <- x[, 4]
          x[, 1:3] <- sprintf("%.6f", x[, 1:3])
          x[, 4] <- sprintf("%.2f", p)
          print(x, quote = FALSE, right = TRUE)
        } 
      }),
      print = FALSE)

[1] "print.summary.mediate"

Code

mediation:::print.summary.mediate(summary(med_6.5))


Causal Mediation Analysis 

Quasi-Bayesian Confidence Intervals

                          Estimate 95% CI Lower 95% CI Upper p-value
ACME (control)           -0.017761    -0.126301     0.034654    0.90
ACME (treated)           -0.013748    -0.092601     0.026874    0.90
ADE (control)             0.004830    -0.239519     0.165314    0.82
ADE (treated)             0.008843    -0.239709     0.160117    0.82
Total Effect             -0.008918    -0.230205     0.178668    0.95
Prop. Mediated (control)  0.173146    -9.108290     8.817234    0.65
Prop. Mediated (treated)  0.110865    -6.130542     6.182164    0.65
ACME (average)           -0.015754    -0.102837     0.028853    0.90
ADE (average)             0.006836    -0.239618     0.157789    0.82
Prop. Mediated (average)  0.142006    -7.661888     7.648164    0.65

Sample Size Used: 1955 


Simulations: 3000

Code

untrace(mediation:::print.summary.mediate)

The interim solution is to print the summary twice, once with the original function for expected p-value behavior, and another time with the edited function for confidence-interval consistency. Someone more well-versed might want to take a shot at adapting the code so that the behavior is maintained for the p-values⁷ and fixed for the confidence intervals.

Conclusion

Mediation analysis is always tricky business⁸, but it feels especially discouraging when the software you need doesn’t seem to want to cooperate. We hope this post will at least shorten someone else’s journey into getting their mediation and/or survey models to behave.

Session Info

Code

sessionInfo()

R version 4.2.2 (2022-10-31)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Ventura 13.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] grid      stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] marginaleffects_0.11.1 gt_0.8.0               srvyr_1.2.0           
 [4] survey_4.1-1           survival_3.4-0         mediation_4.5.0       
 [7] sandwich_3.0-2         mvtnorm_1.1-3          Matrix_1.5-1          
[10] MASS_7.3-58.1          here_1.0.1             lubridate_1.9.2       
[13] forcats_1.0.0          stringr_1.5.0          dplyr_1.1.0           
[16] purrr_1.0.1            readr_2.1.4            tidyr_1.3.0           
[19] tibble_3.1.8           ggplot2_3.4.1          tidyverse_2.0.0       

loaded via a namespace (and not attached):
 [1] sass_0.4.4        jsonlite_1.8.4    splines_4.2.2     Formula_1.2-5    
 [5] yaml_2.3.6        pillar_1.8.1      backports_1.4.1   lattice_0.20-45  
 [9] glue_1.6.2        digest_0.6.31     checkmate_2.1.0   minqa_1.2.5      
[13] colorspace_2.1-0  htmltools_0.5.4   lpSolve_5.6.18    pkgconfig_2.0.3  
[17] scales_1.2.1      tzdb_0.3.0        lme4_1.1-33       timechange_0.2.0 
[21] htmlTable_2.4.1   generics_0.1.3    ellipsis_0.3.2    withr_2.5.0      
[25] nnet_7.3-18       cli_3.6.1         magrittr_2.0.3    evaluate_0.19    
[29] fansi_1.0.4       nlme_3.1-160      foreign_0.8-83    tools_4.2.2      
[33] data.table_1.14.6 hms_1.1.2         mitools_2.4       lifecycle_1.0.3  
[37] munsell_0.5.0     cluster_2.1.4     compiler_4.2.2    rlang_1.0.6      
[41] nloptr_2.0.3      rstudioapi_0.14   htmlwidgets_1.6.2 base64enc_0.1-3  
[45] rmarkdown_2.19    boot_1.3-28       gtable_0.3.1      DBI_1.1.3        
[49] R6_2.5.1          gridExtra_2.3     zoo_1.8-11        knitr_1.41       
[53] fastmap_1.1.0     utf8_1.2.3        Hmisc_5.1-0       rprojroot_2.0.3  
[57] insight_0.19.1    stringi_1.7.8     Rcpp_1.0.10       vctrs_0.5.2      
[61] rpart_4.1.19      tidyselect_1.2.0  xfun_0.35

Footnotes

Previously called the “Fragile Families and Child Wellbeing Study”.↩︎
For more information about the weights, have a look at the methodology document linked here.↩︎
To be 100% sure, Vasco also posted this question to StackOverflow, hoping for clarification.↩︎
Another option, which we won’t go into, would be to take the responses to the items themselves and use a multi-level ordinal model with item as the random grouping factor.↩︎
Of course, the effects mentioned here are not really the effects of moving from X to Y on the mediator, but the mediated effect that changing from Race = “White” to Race = “Black” has on the probability of IPV through its effect on the mediator.↩︎
Which Lumley apparently has also worked on!↩︎
Perhaps taking into account user2554330’s insight in another answer about how p-values are stored and printed.↩︎
See, e.g., That’s a Lot to Process! Pitfalls of Popular Path Models by Rohrer and colleagues.↩︎

Reuse

CC BY-NC 4.0

Citation

BibTeX citation:

@online{brazão2023,
  author = {Brazão, Vasco and Devendran, Priya},
  title = {Mediation Analysis with Weighted {GLMs} in {R:} {What} Could
    Possibly Go Wrong?},
  date = {2023-12-21},
  langid = {en}
}

For attribution, please cite this work as:

Brazão, V., & Devendran, P. (2023, December 21). Mediation analysis with weighted GLMs in R: What could possibly go wrong?

Background

Issue 1: Correctly applying the Fragile Families weights

Running the weighted models and mediation

Issue 2: Mediation with “multiple-trial binomial mediator” runs for weighted but not for unweighted models?

Sanity check: comparing results with linear probability models

Bonus finding: you can specify a mediator that is not technically the variable that your mediator model predicts

Issue 3: Using a binary mediator and a binomial outcome — works with cbind(successes, failures) but not with proportion as outcome and weights argument specified

Issue 4: Ordinal mediator and binomial outcome don’t work together

Issue 5: mediate() is not compatible with svyolr() models

Issue 6: Options and choices to boot — To robustSE or not to robustSE?

Bonus: this goes for using marginaleffects as well

Issue 7: Bounds of the confidence intervals are rounded inconsistently

Conclusion

Session Info

Footnotes

Reuse

Citation

Bonus: this goes for using `marginaleffects` as well