Mediation analysis with weighted GLMs in R: What could possibly go wrong?

Issues (and some solutions) when working with R packages {mediation} and {survey}

statistics
mediation
survey
glm
fragile families
Authors

Vasco Brazão

Priya Devendran

Published

December 21, 2023

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 Study1 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.2 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.3

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”.4

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.5

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:

  1. boot = TRUE and, automatically, robustSE = FALSE
  2. boot = FALSE and robustSE = TRUE, and
  3. 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 package6, 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-values7 and fixed for the confidence intervals.

Conclusion

Mediation analysis is always tricky business8, 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

  1. Previously called the “Fragile Families and Child Wellbeing Study”.↩︎

  2. For more information about the weights, have a look at the methodology document linked here.↩︎

  3. To be 100% sure, Vasco also posted this question to StackOverflow, hoping for clarification.↩︎

  4. 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.↩︎

  5. 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.↩︎

  6. Which Lumley apparently has also worked on!↩︎

  7. Perhaps taking into account user2554330’s insight in another answer about how p-values are stored and printed.↩︎

  8. See, e.g., That’s a Lot to Process! Pitfalls of Popular Path Models by Rohrer and colleagues.↩︎

Reuse

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?