kosukeimai / mediation Goto Github PK

When trying to run the following code for non-parametric mediation the function "mediate.binary" no longer works. Has the name of this function been updated or changed? Below is the replication code for 'figure 3' that I am having issues with
Imai, Keele, & Tingley (2010):

#Initialize set of models

a<-lm(depress2 ~ treat + depress1 +age+educ+sex+nonwhite+econ_hard+marital+occp+income, data=job)
b <- lm(job_seek ~ treat + depress1 +age+educ+sex+nonwhite+econ_hard+marital+occp+income, data=job)
c <- gam(depress2 ~ treat +s(job_seek, bs="cr")+ depress1+age+educ+sex+nonwhite+econ_hard+marital+occp+income , data=job)

#No interaction assumption

mod.1 <- mediate.binary(b, c, sims=sims, boot=TRUE, T="treat", M="job_seek")
summary(mod.1)

Error in mediate.binary(b, c, sims = sims, boot = TRUE, T = "treat", M = "job_seek") :
could not find function "mediate.binary"

The demo for code here for gam interactions no longer runs and gives error "Error in predict.gam(new.fit.Y, type = "response", newdata = pred.data.t) : 0 not in original fit". This appears to be an error that happens when predict.gam is given unknown factor levels.

Strangely if I comment out control as below it runs or if I set "control=treat" it runs. However I get different results in with these two methods. The commenting out method gives a summary result that is closer to the old vignette:

model.m <- lm(job_seek ~ treat + depress1 + econ_hard + sex + age + occp + marital
                 + nonwhite + educ + income, data = jobs)

model.y <- gam(depress2 ~ treat + s(job_seek, by = treat) + s(job_seek, by = control)
               + depress1 + econ_hard + sex + age + occp + marital + nonwhite + educ + income,
               data = jobs)

out.5 <- mediation::mediate(model.m, model.y, sims = 1000, boot = TRUE,
                            treat = "treat", mediator = "job_seek",
                            #control = "control"
                       )
summary(out.5)

However the result isn't close enough for me to feel confident that this is correct (e.g. within .01).

Both the paper referenced in mediate's documentation for gam and the old vignette imply that the by variables should be indicator variables. The current jobs data does not have this form: treat is binary but control is a factor. If I create a new variable jobs$control.1 <- 1 - jobs$treat I can then set the by options sensibly, but again I don't get the values from the old vignette. I've also tried just using the spline by factor control since the gam documentation suggests that creates an interaction by itself anyways. Altogether I've explored a number of alternatives, and the closest to the old vignette seems to be the "comment out control" listing above.

How can I perform a mediation analysis with a continuous treatment. I want to esimate indirect effect when exposure increase 1 unit instead of specifying control.value and treat.value.

Can you please create a conda wheel so I can install this R package using anaconda? Thanks!

The mediate function does not work when package lmerTest is loaded, because lmerTest changes the class of the model from lmerMod to lmerModLmerTest.
The function returns the error: mediator model is not yet implemented

Can the function be upgraded to be compatible with lmerTest package?
Thanks

Hi,

I was wondering if there was a potential bug with how headers are reported out for the summary of a mediation object?

I have a 4-level ordinal outcome (No (1), Mild (2), Moderate (3), Severe (4)) analyzed as a factor, a continuous mediator & a binary treatment.

I ran a glm for my mediator ~ treatment & a polr for my outcome ~ treatment + mediator.

When I run the mediation analysis using the 4-level factor outcome where the levels are coded as text, and run summary(), it reports out this way:

Causal Mediation Analysis

Nonparametric Bootstrap Confidence Intervals with the Percentile Method

            Pr(Y=Mild) Pr(Y=Moderate) Pr(Y=No) Pr(Y=Severe)

ACME (control) -0.01801 0.00874 0.00447 0.00480
2.5% -0.02936 0.00262 0.00134 0.00139
97.5% -0.00561 0.01455 0.00773 0.00842
p-value 0.00200 0.00200 0.00200 0.00200

But when I change my outcome to a 4 level factor outcome where the levels are coded as numbers (1 (No), 2 (Mild), 3 (Moderate), 4 (Severe), it reports out this way:

Causal Mediation Analysis

Nonparametric Bootstrap Confidence Intervals with the Percentile Method

             Pr(Y=1) Pr(Y=2) Pr(Y=3) Pr(Y=4)

ACME (control) -0.01801 0.00874 0.00447 0.00480
2.5% -0.02936 0.00262 0.00134 0.00139
97.5% -0.00561 0.01455 0.00773 0.00842
p-value 0.00200 0.00200 0.00200 0.00200

Is this expected behavior or have I made a mistake?

Thanks!

-Stephanie

When calling summary on a mediate object, the underlying code uses the printCoefmat method. By default, this assumes the results of the second-last column (i.e. the upper-confidence interval here) are a test statistic, which is reported with less precision than other values such the lower confidence interval. That is inconsistent between the bounds, but can also lead to rounding to 0, which can lead to inconsistencies between the confidence intervals and the p-values. For instance

> dat <- read.csv("https://mspeekenbrink.github.io/sdam-jasp-companion/data/redist2015.csv")
> mod2 <- lm(SC_mean ~ income, data=dat)
> mod3 <- lm(redist ~ SC_mean + income, data=dat)
> set.seed(20201027)
> med <- mediate(model.m = mod2, model.y = mod3, sims=1000, boot = TRUE, boot.ci.type = "bca", treat = "income", mediator = "SC_mean")
Running nonparametric bootstrap

> summary(med)

Causal Mediation Analysis 

Nonparametric Bootstrap Confidence Intervals with the BCa Method

               Estimate 95% CI Lower 95% CI Upper p-value    
ACME           -0.00229     -0.00463         0.00   0.014 *  
ADE            -0.00292     -0.00650         0.00   0.078 .  
Total Effect   -0.00521     -0.00801         0.00  <2e-16 ***
Prop. Mediated  0.43953      0.12895         1.29   0.014 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Sample Size Used: 301 


Simulations: 1000 

It would be better if the upper confidence interval had the same precision as the lower CI. An easy fix (I believe) is to indicate that there is no test statistic, replacing all calls to printCoefmat to

printCoefmat(smat, tst.ind=NULL)

Removing the digits=3 argument that is used at the moment also allows the digits to be controlled by the general options, rather than hardcoded. E.g.

> x <- med
> smat <- c(x$d1, x$d1.ci, x$d1.p)
> smat <- rbind(smat, c(x$z0, x$z0.ci, x$z0.p))
> smat <- rbind(smat, c(x$tau.coef, x$tau.ci, x$tau.p))
> smat <- rbind(smat, c(x$n0, x$n0.ci, x$n0.p))
> rownames(smat) <- c("ACME", "ADE",
+                     "Total Effect", "Prop. Mediated")
> clp <- "95"
> colnames(smat) <- c("Estimate", paste(clp, "% CI Lower", sep=""),
+                     paste(clp, "% CI Upper", sep=""), "p-value")
> printCoefmat(smat, tst.ind=NULL)
                  Estimate 95% CI Lower 95% CI Upper p-value    
ACME           -0.00228812  -0.00462947  -0.00046455   0.014 *  
ADE            -0.00291772  -0.00649948   0.00024205   0.078 .  
Total Effect   -0.00520584  -0.00801499  -0.00287072  <2e-16 ***
Prop. Mediated  0.43953008   0.12894917   1.29021831   0.014 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Hi,

I just re-ran an old analysis which functioned without problems in November 2020 but returns this error now Error in eval(extras, data, env) : object 'Call.M' not found.

I did not change the code at all.

I already tried reinstalling the package, but that didn't help.

Could you help me with this issue?

Best regards
Leo

I am conducting a mediation analysis with multiply imputed 30 data sets using "amelidiate {mediation}".

variables: mediators= m,outcome=y, treatment=rep(e,30),covariates=z, control_val=a0, treat_val=a1,

medout <- mediations(datasets, treatment, mediators, outcome, covariates,
families=c("binomial","gaussian"), control.value = control_val,
treat.value = treat_val, interaction=TRUE, conf.level=.95, sims=1000, set.seed(1))

med_pool <- amelidiate(medout)
summary(med_pool)

The package provides pooled CI for all estimates but does not provide pooled p-values. Would there any way to calculate pooled p-values?

Thank you so much!

I'm having trouble with mediation:: medsens()

med.fit <- glm(ADM_ROTSUP_0_4S ~ fgrupo, data = dados, family="gaussian")

out.fit <- lm(VAS.0.8 ~ ADM_ROTSUP_0_4S +fgrupo + VAS.0 + ADM_ROTSUP_0_INICIAL + Idade..anos. + 
Tempo.de.dor..meses. + fsexo, data = dados)

med.out <- mediate(med.fit, out.fit, boot = TRUE, treat = "fgrupo", mediator = "ADM_ROTSUP_0_4S", robustSE = TRUE, sims = 1000)

The code works fine until med.out, but when I try to run a sensitivity analysis using medsens I get:

sens.out <- medsens(med.out, rho.by = 0.1, sims = 1000, effect.type = "indirect")
summary(sens.out)
ERROR: subscript out of bounds 

Someone might help me?

Thanks!

Mediate returns NaNs for the CIs if I run a model with ordered outcomes and the bca option.
Here a reproducible example:

library(mediate)
library(MASS)

dv <- sample(1:4, size=200, replace=TRUE)
m <- sample(1:4, size=200, replace=TRUE)
t <- sample(0:1, size=200, replace=TRUE)

ex.dat <- data.frame(dv,m,t)
ex.dat$dv <- as.factor(ex.dat$dv)
ex.dat$m <- as.factor(ex.dat$m)

med.m <- polr(m ~ t, data = ex.dat)
summary(med.m)

out.m <- polr(dv ~ t + m, data = ex.dat)
summary(out.m)

out.ord <- mediate(med.m, out.m, treat = "t", 
                    mediator = "m", sims = 100, 
                    boot = T, boot.ci.type = "bca")
summary(out.ord)

It produces an output, but the CIs are missing.

It returns the following warning:

Warning messages:
1: In qnorm(z.inv) : NaNs generated
2: In qnorm(z.inv) : NaNs generated
3: In qnorm(z.inv) : NaNs generated
4: In qnorm(z.inv) : NaNs generated
5: In qnorm(z.inv) : NaNs generated

Hi!

Thank you very much for providing the mediation package to R, I really enjoy working with it! I was wondering, however, if it is somehow possible to extend the mediation package to also accept GAMM models, in one way or another.

Being aware that "mediation" does not handle "gamm" models directly, I tried to circumvent the problem by implementing my design with the mgcv "gam" function by using bs="re" flag as specification (as suggested by S. Wood , https://stat.ethz.ch/R-manual/R-devel/library/mgcv/html/random.effects.html). Unfortunately, it appears that using this approach the "predict.gam" (as currently used in "mediation") fails, and I get the following error message:

"Error in predict.gam(new.fit.M, type = "response", newdata = pred.data.t) :
403, 407, 418, 421,..."

In advance, thank you for your time!

Best regards
René

Hello,
Thank you for developing such a great package. It is highly useful in implementing many different kinds of mediation analyses, and I'm getting a lot of help for my projects.
But recently, I ran into trouble while using the "ivmediate" function.

The codes I ran was like this:
`############################################## IV Mediation with mediation package ##############################

library(mediation)
library(car)
library(fastDummies)

treatment = 'rh_adi_bi_0y'
instrument = 'section8_0y'
mediator = 'discount_rate_1y_z'
outcome = 'totalscore_pps_1y_z'

data<-dummy_cols(data, select_columns = c('sex_0y'))
data$sex_0y<-as.factor(data$sex_0y)
data$married_0y<-as.factor(data$married_0y)
data$race_g<-as.factor(data$race_g)
data$parent_identity_0y<-as.factor(data$parent_identity_0y)
data$gender_identity_0y<-as.factor(data$gender_identity_0y)
data$foreign_born_0y<-as.factor(data$foreign_born_0y)
data$foreign_born_family_0y<-as.factor(data$foreign_born_family_0y)
data$gay_parent_0y<-as.factor(data$gay_parent_0y)
data$gay_youth_0y<-as.factor(data$gay_youth_0y)
data$race_ethnicity_0y<-as.factor(data$race_ethnicity_0y)

fmla1 <- as.formula("rh_adi_bi_0y~section8_0y+nihtbx_totalcomp_uncorrected_0y_z+iqeur2+
vol_z+age_0y_z+bmi_0y_z+income_0y_z+high_educ_0y_z+parent_age_0y_z+history_ratio_0y_z+vol_z+
sex_0y_F+married_0y+parent_identity_0y+gender_identity_0y+
foreign_born_family_0y+foreign_born_0y+gay_parent_0y+gay_youth_0y+race_ethnicity_0y")

fmla2 <- as.formula("discount_rate_1y_z~rh_adi_bi_0y+section8_0y+nihtbx_totalcomp_uncorrected_0y_z+iqeur2+
vol_z+age_0y_z+bmi_0y_z+income_0y_z+high_educ_0y_z+parent_age_0y_z+history_ratio_0y_z+vol_z+
sex_0y_F+married_0y+parent_identity_0y+gender_identity_0y+
foreign_born_family_0y+foreign_born_0y+gay_parent_0y+gay_youth_0y+race_ethnicity_0y")

fmla3 <-as.formula("totalscore_pps_1y_z ~ discount_rate_1y_z*(rh_adi_bi_0y+section8_0y)+nihtbx_totalcomp_uncorrected_0y_z+iqeur2+
vol_z+age_0y_z+bmi_0y_z+income_0y_z+high_educ_0y_z+parent_age_0y_z+history_ratio_0y_z+vol_z+
sex_0y_F+married_0y+parent_identity_0y+gender_identity_0y+
foreign_born_family_0y+foreign_born_0y+gay_parent_0y+gay_youth_0y+race_ethnicity_0y")

model_treatment <- glm(fmla1, data = data, family = binomial)
model_mediator <- lm(fmla2, data = data)
model_outcome <- lm(fmla3, data = data)

ivmediation_result <- ivmediate(model_treatment, model_mediator, model_outcome, sims = 500,
enc = 'section8_0y', treat='rh_adi_bi_0y', mediator='discount_rate_1y_z')

summary(ivmediation_result)
plot(ivmediation_result)`

and the error messages that I got were like this:

Error in integrate(LACME1.integrand, -Inf, Inf, i = i, stop.on.error = F) : evaluation of function gave a result of wrong type In addition: Warning messages: 1: Using formula(x) is deprecated when x is a character vector of length > 1. Consider formula(paste(x, collapse = " ")) instead. 2: Using formula(x) is deprecated when x is a character vector of length > 1. Consider formula(paste(x, collapse = " ")) instead. 3: Using formula(x) is deprecated when x is a character vector of length > 1. Consider formula(paste(x, collapse = " ")) instead. 4: Using formula(x) is deprecated when x is a character vector of length > 1. Consider formula(paste(x, collapse = " ")) instead. 5: glm.fit: algorithm did not converge 6: Using formula(x) is deprecated when x is a character vector of length > 1. Consider formula(paste(x, collapse = " ")) instead. 7: Using formula(x) is deprecated when x is a character vector of length > 1. Consider formula(paste(x, collapse = " ")) instead.

I think I've done the same thing as shown in the example below, but this just doesn't seem to work for some reason.
https://www.rdocumentation.org/packages/mediation/versions/4.5.0/topics/ivmediate

`a <- lm(comply ~ treat + sex + age + marital + nonwhite + educ + income,
data = jobs)
b <- glm(job_dich ~ comply + treat + sex + age + marital + nonwhite + educ + income,
data = jobs, family = binomial)
c <- lm(depress2 ~ job_dich * (comply + treat) + sex + age + marital + nonwhite + educ + income,
data = jobs)

out <- ivmediate(a, b, c, sims = 50, boot = FALSE,
enc = "treat", treat = "comply", mediator = "job_dich")

summary(out)
plot(out)`

Have I missed something, or is there some problem with the ivmediate function that makes ends up with this error?

I'll really appreciate it if you can help.
Thanks in advance!

I am having some trouble using your mediate() function with a Generalized Additive Model(GAM) and see that your documentation specifies it can be done, but does not explain the appropriate syntax. Specifically, I am trying to run a causal mediation analysis with a GAM with an interaction between a continuous variable (time) and two factor categorical variable (illness group). The model aims to asses the relationship of how gene expression (continuous variable) differs over time (also a continuous variable) by a two factor categorical variable (Illness group) and whether this is mediated by viral load (a continuous variable), which differs overtime between the two illness groups. The model syntax for the GAMs is here:

model.mediate <- gam(viral_load~illnessgroup+ s(time) +s(time, by = illnessgroup), data=data, method = "REML")

model.outcome <- gam( gene_exp~illnessgroup + viral_load + s(time) +s(time, by = illnessgroup) ,data=data, method = "REML")
mediation_result <- mediate(model.m = model.mediate, model.y = model.outcome, mediator = "viral_load", treat = "illnessgroup", treat.value = "Severe", control.value = "Not_Severe", boot = TRUE, control = "Not_Severe")

This gives me an output with the average ACME, ADE, Prop. mediated and Total Effects but I was expecting to see separate statistics for the treatment (Severe) and control (Not Severe) conditions. This is what we would expect if using a linear model rather than a GAM. Any advice on how to appropriately run the mediate() function with a GAM model using an interaction for 2 comparison groups to get the output for both the treatment and control groups would be much appreciated.
Thank you,
Basilin Benson

Hello,

I am trying to run the mediate function using a mediator model of class "lmer" (from lme4) and an outcome model of class "gamm", due to nested random effects. However, only a class "gam" is compatible. Is it possible to add "gamm" as an allowable model class?

Thanks!
Ashley

Hello,

I am performing a mediation analysis with mixed models in the presence of a (potential) treatment * mediator interaction. I would like to test if this interaction is significant. Is there anything wrong with using the same logic in this code to preform this test? My model.m and model.y that I provide to the mediate function are both lmer objects, and test.TMint says that it is not implemented for these models quite yet.

So my question is, is there anything theoretically wrong with this:

mod.m <- lmer(M ~ X + (1|id), data = data)
mod.y <- lmer(Y ~ M * X + (1|id), data = data)

med.results <- mediate(mod.m, mod.y, mediator = "M", treat = "X" )

acme_diff = med.results$d1.sims - med.results$d0.sims

quantile(acme_diff, probs = c(0.025, 0.975))

I believe there must be an error which is causing the sign of the reported ACME effect to change depending on the contrast coding matrix for X, when both X->M (a) and M->Y (b) are negative, regardless of how control.value and treat.value are specified in the call to mediate.

In the experimental data during the analysis of which I encountered this issue, X represents a two-contrast factor with english and dutch as the levels. I am interested in setting english as the control and testing for the effect of dutch compared with english. Note that this implies changing the factor contrast coding from the default, because E comes after D, and would usually default as the treatment level.

I simulate these data below:

simulate_data = function(x1){
  # x2, x3, and x4 in a matrix, these will be modified to meet the criteria
  x2 <- scale(matrix( rnorm(length(x1)), ncol=1))
  
  # put all into 1 matrix for simplicity
  x12 <- cbind(scale(x1),x2)
  
  # find the current correlation matrix
  c1 <- var(x12)
  
  # cholesky decomposition to get independence
  chol1 <- solve(chol(c1))
  
  newx <-  x12 %*% chol1 
  
  # create new correlation structure (zeros can be replaced with other rvals)
  newc <- matrix(c(1, -0.5,
                   -0.5, 1), 
                 ncol=2)
  
  chol2 <- chol(newc)
  
  finalx <- newx %*% chol2 * sd(x1) + mean(x1)
  
  finalx[,2]
  }


data = data.frame("X"=c(rep('B_english',300),rep('A_dutch',300)),
                  "M"=scale(c(rnorm(300,mean=150, sd = 50),
                                           rnorm(300,mean=100,sd=50))
                    ))
data$Y = simulate_data(data[,'M'])

#set contrasts
data$X_rev = factor(data$X)
contrasts(data$X_rev) = contr.treatment(2,base=2)
colnames(contrasts(data$X_rev)) = 'A_dutch'

data$X = factor(ifelse(data$X == 'B_english','A_english','B_dutch'))

> contrasts(data$X)
          B_dutch
A_english       0
B_dutch         1
> contrasts(data$X_rev)
          A_dutch
A_dutch         1
B_english       0

The only difference between data$X and data$X_rev is the letter they start with. In both of the contrast matrices, the level representing dutch is coded as the treatment contrast, so I would expect that this difference should not result in any substantive changes to model outputs. However this does not appear to be true.

When X = data$X, the sign for ACME is correct (i.e. positive since a and b are both negative):

### X = data$X
med.fit = lm(M ~  X,data = data)
out.fit = lm(Y ~ M + X,data = data)
med.out = mediate(med.fit,out.fit,treat='X',mediator='M',sims=100,
                  control.value = 'A_english', treat.value = 'B_dutch')
>summary(med.out)
Causal Mediation Analysis 
Quasi-Bayesian Confidence Intervals
               Estimate 95% CI Lower 95% CI Upper p-value    
ACME             0.4568       0.3632         0.55  <2e-16 ***
ADE             -0.0764      -0.2247         0.07    0.38    
Total Effect     0.3804       0.2186         0.54  <2e-16 ***
Prop. Mediated   1.2186       0.8517         1.98  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Sample Size Used: 600 
Simulations: 100 

However, in the case of X = data$X_rev, the ACME is reported as negative when it should be positive:

med.fit.rev = lm(M ~  X_rev,data = data)
out.fit.rev = lm(Y ~ M + X_rev,data = data)
med.out.rev = mediate(med.fit.rev,out.fit.rev,treat='X_rev',mediator='M',sims=100,
                  control.value = 'B_english', treat.value = 'A_dutch')
> summary(med.out.rev)
Causal Mediation Analysis 
Quasi-Bayesian Confidence Intervals
               Estimate 95% CI Lower 95% CI Upper p-value    
ACME            -0.4670      -0.5757        -0.35  <2e-16 ***
ADE              0.0787      -0.0637         0.19    0.24    
Total Effect    -0.3883      -0.5335        -0.22  <2e-16 ***
Prop. Mediated   1.2161       0.8620         1.83  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Sample Size Used: 600 
Simulations: 100 

This would seem to be undesirable and/or confusing behaviour. I'm sorry if I missed somewhere in the docs where it specifies that factor contrasts must be coded to the R default (i.e. in alphabetical order). I would suggest that this point either be made more salient in the vignette, or that this behaviour be corrected.

Using the mediate function (please see my code below), I receive the following error message. I would like to estimate a multilevel mediation for the following 2x5 mixed design. While the IV is manipulated between-subjects, the DV is a repeated measure (binary). In particular, each participant made 5 binary choices for 5 different choice sets (within-subject). The mediator (i.e., perceived expertise) is a repeatedly measured continuous variable, which was measured after each choice.

Do you have any idea what could be wrong with my code?

Error in if (xhat == 0) out <- 1 else { :
missing value where TRUE/FALSE needed

med.fit <- lmer(Expertise ~ Condition + (1|SubjectID), data = data)
out.fit <- glmer(Choice ~ Expertise + Condition + (1 + Expertise|SubjectID),
family = binomial(link = "logit"), data = data)

med.out <- mediate(med.fit, out.fit, treat = "Condition", mediator = "Expertise",
sims = 5000)

summary(med.out)

The data looks like the following:
SubjectID I ChoiceSet I Condition I Expertise I Choice
1 I 1 I 0 I 5 I 1
1 I 2 I 0 I 4 I 1
1 I 3 I 0 I 4 I 0
1 I 4 I 0 I 4 I 0
1 I 5 I 0 I 3 I 0
2 I 1 I 1 I 3 I 1
2 I 2 I 1 I 3 I 1
2 I 3 I 1 I 4 I 1
2 I 4 I 1 I 4 I 1
2 I 5 I 1 I 5 I 0

I recently updated R and my packages and the below code, which was previously working, is no longer working:

fit.dv <- lme4::lmer(opti3 ~ no_recalled_notes + N + condition + log(attempt) + log(attempt):N + (1|unique_melody_name), data = na.omit(main2))

fit.mediator <- lme4::lmer(no_recalled_notes ~ N + condition + log(attempt) + log(attempt):N + (1|unique_melody_name), data = na.omit(main2))

results <- mediation::mediate(fit.mediator, fit.dv, treat = c('N', 'condition'), mediator = 'no_recalled_notes')

The error is:

Error in if(!INT & isGam.y) { : the condition has length > 1

I'm running two mediations on the same data, once using prior (sampling) weights when fitting the models (through the survey package), and once without those prior weights (just using the glm() function). My outcome is binary. My mediator is a proportion out of 3, so I also pass 3 as the weights argument (or rather, the name of a column that is equal to 3 for every subject). When using the models fit with survey:svyglm() inside mediate(), there is no issue; however, when using the glm() models, mediate() complains that "weights on outcome and mediator models not identical".

After checking the documentation, I see that you're not supposed to use a multiple-trial binomial mediator, as that is not implemented yet. My question is, why does it work when the models are fit through survey::svyglm()? And can I trust those results at all?

[Error in Surv(time, event) : cannot find object "time"] in survreg()

Here I am doing mediation analysis among X, Y, and m.
Here model.m is linear model and model.y is survival model.
The code is as follow:

model.m <- lm(m ~ X, data = data)
model.y <- survreg(Surv(time, event) ~ X+m, data = data)
mediation <- mediate(model.m, model.y, treat = "X", mediator = "m" , boot = T, outcome = "event")

However, error occurred as Error in Surv(time, event) : cannot find object "time.

This has caused some problems from some users.

Can be replicated from the package example file (below). I suggest using the same rounding for upper as we do for lower.

summary(contcont)

Causal Mediation Analysis

Quasi-Bayesian Confidence Intervals

           Estimate 95% CI Lower 95% CI Upper p-value  

ACME -0.0183 -0.0397 0.00 0.04 *
ADE -0.0450 -0.1417 0.03 0.40
Total Effect -0.0634 -0.1636 0.01 0.16
Prop. Mediated 0.2686 -0.9015 2.54 0.20

This works fine with boot=F, but throws an error when boot=T.
Seems related to #29

library(mediation)
#> mediation: Causal Mediation Analysis
#> Version: 4.5.0
test <- function(f1, f2, boot){
    mtcars$mpg <- mtcars$mpg/max(mtcars$mpg)
    m_med <- lm(f1, mtcars)
    m_out <- lm(f2, mtcars)
    mediate(m_med, m_out, treat = "wt", mediator = 'am', boot = boot)
}

x <- test('am ~ wt', 'mpg ~ wt + am', F)
#success!
x <- test('am ~ wt', 'mpg ~ wt + am', T)
#> Running nonparametric bootstrap
#> Error in stats::model.frame(formula = f1, data = structure(list(am = c(1, : object 'f1' not found

Is it possible to use the mediate function for multi-level outcome variables?

I have what seems a simple case that cannot properly fit.

I have a dependent variable that is measured 10 times per subject. However, the mediator is only once after the intervention (there are two groups: Control and Interventions)

I thought that the following would work:

model_m <- lm(mediator ~ group, data=unique (df) )
model_dv <- lmer(dv ~ group + mediator  + (1 | subj), data = df)
mediate(model_m, model_dv, treat='group', mediator='mediator', boot=F)

To be sure: model_m is fit using lm because the mediator is only measured once (note that data=unique(df)), whereas model_dv is fitted using lmer because it is measured 10 times (within subject).

It triggers this error, which makes sense: one group has 10 times more rows.

Instead, I tried the following:

model_m <- lmer(mediator ~ group + (1|subj), data=df )
model_dv <- lmer(dv ~ group + mediator  + (1 | subj), data = df)
mediate(model_m, model_dv, treat='group', mediator='mediator', boot=F)

This works except that model.m is incorrect because it is fitted in 10 times the measurements that actually were taken. As a consequence, there are several warnings displayed (e.g., model failed to converge which makes sense), *the estimates are the same as the ones obtained with lm but the the standard errors are smaller.

Is it OK to still the model_m (fitted with lmer) even though some of the statistics are wrong? In other words, which statistics of the lmerMod object are used inside the mediate function? (I can then check if they are all the same and it is just about the format of the data)

I want to do a mediation analysis, the outcome model is a survival model, but the document of this packages shows the outcome model can be survreg object instead of coxph object. Cox ph model can be used as outcome model?

Hi there,
thanks for the package. It is insightful and helpful.
My outcome variable has different characteristics for the bottom and upper quantiles. Unfortunately, multimed function doesn't have tau.y argument unlike mediate function. Do you have any idea how to deal with this issue?

Sincerely,

Hi!

I love the package, thanks for it (and the theoretical work underlying it)

Is it possible to extract a standard path diagram from the output of mediation()?
I believe reviewers will require these standard graphs because the readers are used to those.

Thanks for the help!
Kind regards
Jana

As title.

Hi, I know there is a function mediations() to perform "multiple" mediation models, and each of them is with a single mediator.

However, if I want to perform a mediaiton model with multiple mediators (like multiple mechanisms in a model, and they are uncorrelated with each other), how could I do that in this package? Thanks!

Hi! Thank you so much for this wonderful package. I was wondering if it would be possible to add support for models that require cluster-robust estimates, specifically lm_robust models from the estimatr package?

Hello,
I'm running a mediation analysis and trying to increase the number of simulations, but I run into an error when my number of simulations reaches n ~= 700.

M2  <- lm( mediator~ temperature, dsite )
M3  <- glmer( rate~  temperature+ mediator + (1|Site), d, family="binomial" )

med <- mediation::mediate( M2,M3, sims=700, treat="temperature", mediator="mediator" )
summary(med)

This produces the following error: Error in if (xhat == 0) out <- 1 else { : missing value where TRUE/FALSE needed

Note that the mediator and treatment here are both continuous variables, so the medation model is a linear model. The mediator model also comes from data at a higher level than the outcome model, so the datasets are slightly different. Interesting though that no error is produced when the number of simulations is lower (e.g., 500). Can you help me understand why I might be getting these errors?

Here is the output of a working mediation model with a low number of simulations:

Causal Mediation Analysis 

Quasi-Bayesian Confidence Intervals

Mediator Groups: 

Outcome Groups: Site 

Output Based on Overall Averages Across Groups 

                          Estimate 95% CI Lower 95% CI Upper p-value    
ACME (control)            4.49e-03     2.28e-06         0.03  <2e-16 ***
ACME (treated)            3.80e-03     1.77e-06         0.02  <2e-16 ***
ADE (control)            -2.42e-03    -1.65e-02         0.00    0.02 *  
ADE (treated)            -3.11e-03    -2.06e-02         0.00    0.02 *  
Total Effect              1.38e-03     1.27e-06         0.01    0.02 *  
Prop. Mediated (control)  1.85e+00     1.05e+00         4.93    0.02 *  
Prop. Mediated (treated)  1.51e+00     1.03e+00         4.23    0.02 *  
ACME (average)            4.15e-03     2.03e-06         0.02  <2e-16 ***
ADE (average)            -2.76e-03    -1.85e-02         0.00    0.02 *  
Prop. Mediated (average)  1.68e+00     1.04e+00         4.56    0.02 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Sample Size Used: 234 


Simulations: 100 

I'm also interested to know what estimates of "Prop. Mediated" is greater than 1. Should I be worried about that? I see some examples in the vignette where at least the upper confidence limit exceeds 1, but I'm not sure what this reveals about the modeling exercise or the underlying data.

Thank you!

Thank you for this great package.

The current mediate() function may be extended further by adding a function that computes odds ratios and/or risk ratios (with confidence intervals) for a dichotomous outcome. In the medical literature, we are often asked to report effect measures in these scales in addition the difference. Thank you very much for your consideration.

I am running into an issue when using mediations command to run the mediation analysis on multiple imputation data. I have weights in each of my datasets, but the packages seems to not be location them, as the error reads:
Error in eval(extras, data, env) : object 'weight' not found

The variable ipw_response_t are the truncated weights and are a numeric column in each of my imputed datasets. I specify it as below:
olsols <- mediations(datasets, treatment, mediators, outcome, covariates,
families=c("gaussian","gaussian")
, interaction=FALSE
, conf.level=.95
, sims=100,weights=c("ipw_response_t"))

Hello,

I am currently running a causal mediation analysis for a clinical trial where my mediator model uses a negative binomial regression model and my outcome model uses a logistic regression model.

model.x <- glm.nb(total.malaria.incident.events ~ Txarm + offset(years.in.study), data = delivery.analysis)
model.y = glm(compositeBO ~ Txarm + total.malaria.incident.events + offset(years.in.study), family = "binomial", data = delivery.analysis)
m.out = mediate(malaria.glm.nb, model.3.malaria.compositeBO, treat = "Txarm", mediator = "total.malaria.incident.events", boot=T)

However, I get the following error: Error in offset(years.in.study) : object 'years.in.study' not found
Does the mediate package allow for an "offset" function to be in the model because the time each participant contributes varies? If so, how can I correct this error? If not, is there another way to handle these models?

Thanks in advance!

At present, summary(mediation(...)) produces output like this:

Causal Mediation Analysis 

Quasi-Bayesian Confidence Intervals

               Estimate 95% CI Lower 95% CI Upper p-value    
ACME             0.1477       0.0237         0.28   0.018 *  
ADE              0.2998       0.1166         0.50  <2e-16 ***
Total Effect     0.4475       0.2293         0.66  <2e-16 ***
Prop. Mediated   0.3295       0.0741         0.60   0.018 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Sample Size Used: 200 


Simulations: 1000 

It would be nice if this could be supplemented by a tidy data frame, such that this would work:

as_data_frame(summary(mediation(...)))

To produce a data frame with these columns:

• term
• estimate
• p.value
• conf.low
• conf.high
• .prob (would be .95, .99 etc)
• replicates (sims)

(these chosen to match those typically produced by broom::tidy)

Hi,
Im trying to use the multimed function.
I can run the example code for multimed, but when I run it on my dataset, or this 'fake' dataset, it errors (see below).
treat is a 0/1 numeric variable.

What am I doing wrong? I assume its user error?
thanks!

library(palmerpenguins)
library(mediation)
#> Loading required package: MASS
#> Loading required package: Matrix
#> Loading required package: mvtnorm
#> Loading required package: sandwich
#> mediation: Causal Mediation Analysis
#> Version: 4.5.0
library(tidyverse)
library(reprex)

create a fake dataset to show example of error

glimpse(penguins)
#> Rows: 344
#> Columns: 8
#> $ species Adelie, Adelie, Adelie, Adelie, Adelie, Adelie, Adel~
#> $ island Torgersen, Torgersen, Torgersen, Torgersen, Torgerse~
#> $ bill_length_mm 39.1, 39.5, 40.3, NA, 36.7, 39.3, 38.9, 39.2, 34.1, ~
#> $ bill_depth_mm 18.7, 17.4, 18.0, NA, 19.3, 20.6, 17.8, 19.6, 18.1, ~
#> $ flipper_length_mm 181, 186, 195, NA, 193, 190, 181, 195, 193, 190, 186~
#> $ body_mass_g 3750, 3800, 3250, NA, 3450, 3650, 3625, 4675, 3475, ~
#> $ sex male, female, female, NA, female, male, female, male~
#> $ year 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007~

penguins %>%
filter(species!="Gentoo") %>%
mutate(treat = as.numeric(rbernoulli(n=220, 0.5))) %>%
drop_na()->df

summary(df)
#> species island bill_length_mm bill_depth_mm
#> Adelie :146 Biscoe : 44 Min. :32.1 Min. :15.50
#> Chinstrap: 68 Dream :123 1st Qu.:37.8 1st Qu.:17.50
#> Gentoo : 0 Torgersen: 47 Median :40.6 Median :18.40
#> Mean :42.0 Mean :18.37
#> 3rd Qu.:46.0 3rd Qu.:19.10
#> Max. :58.0 Max. :21.50
#> flipper_length_mm body_mass_g sex year treat
#> Min. :172.0 Min. :2700 female:107 Min. :2007 Min. :0.000
#> 1st Qu.:187.0 1st Qu.:3400 male :107 1st Qu.:2007 1st Qu.:0.000
#> Median :191.0 Median :3700 Median :2008 Median :1.000
#> Mean :191.9 Mean :3715 Mean :2008 Mean :0.514
#> 3rd Qu.:196.0 3rd Qu.:3994 3rd Qu.:2009 3rd Qu.:1.000
#> Max. :212.0 Max. :4800 Max. :2009 Max. :1.000

test<-multimed(outcome = "bill_length_mm",
med.main = "flipper_length_mm",
med.alt = "bill_depth_mm",
treat = "treat",
covariates = "year",
data=df, sims=100)
#> Warning in mean.default(data.1[, treat] * data.1[, med.main]^2): argument is not
#> numeric or logical: returning NA
#> Warning in ETM2 * sigma^2/VY: Recycling array of length 1 in vector-array arithmetic is deprecated.
#> Use c() or as.vector() instead.
#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA

#> Warning in mean.default(data.b[, treat]): argument is not numeric or logical:
#> returning NA
Created on 2022-05-18 by the reprex package (v2.0.1)

Could Package "mediate" support multi-categorical treatment?

I could find arguments 'control.value' & 'treat.value' in the function mediate, it seems that the function could support binary treatment. I am not sure whether mediate can support multi-categorical treatment (# groups >= 3) as well. Thank you for your help in advance!

David

I'm having trouble with mediation::mediate()

nb.terror.y.anx <- glm(data = terror_removed2,
                   certain~cond + event + WC + anx,
                   control = glm.control(maxit = 50000),
                   family = "poisson")

nb.terror.m.anx <- glm(data = terror_removed2,
                   anx~cond + event + WC + certain,
                   control = glm.control(maxit = 50000),
                   family = "poisson")

med.anx <- mediation::mediate(nb.terror.m.anx,nb.terror.y.anx, treat = "cond", 
                              mediator = "anx", sim = 2)

Works fine, but whenever sim>2, I get:

Error in if (xhat == 0) out <- 1 else { : 
  missing value where TRUE/FALSE needed

I saw there's still an open topic in a realted issue

Hope it's ok to flag it up again.

Thanks!

Hello,
thank you for making this great package!

I have a dataset with missing values so I am using the mediations() and amelidiate() functions. At the same time, I would also like to test whether the mediated effect differs by demographic so I would like to use the test.modmed() function. Can they be somehow used together?

This is most probably my own misunderstanding, but I'm wondering whether bootstrap CIs are correctly implemented. Two questions:

1. In line 1584 of mediate.py, z_0 is calculated by compare theta to mean theta.

            z.inv <- length(theta[theta < mean(theta)])/sims

However, eq. 2.8 in DiCiccio & Efron, 1996 is different: theta is compared to the original, non-bootstrapped statistic. The current implementation assumes that the mean of the bootstrap distribution equals to the original statistic. Is this necessarily true?

2. In line 1586, U is calculated, in a way that seem to follow eq. 6.7 in the paper:

U <- (sims - 1) * (mean(theta) - theta)

However, eq. 6.7 uses jackknife estimator, not the bootstrap estimator. Is the usage of the bootstrap estimator here justified?

Hi, does medsens function support the sensitivity analysis of causal mediation analysis with multilevel data? Thank you so much!

Due to the essence of multilevel data, I use lmer function to obtain med.fit and out.fit. Then, use these two outputs as inputs of the mediate function to get med.out. Finally, I try to use medsens function to perform the sensitivity analysis.

However, I get the error message below:

Error in medsens(med.out, rho.by = 0.1, effect.type = "indirect", sims = 100) :
mediate object fitted with non-supported model combinations

thx!

I am working on a package that takes user input and within its functions, uses lme4 to fit and return a "merMod" object.

I tried to pass these models to mediate and was getting a "mediator model is not yet implemented" error. Poking through the source code, I see why. There are a series of if statements like

if(isMer.m && getCall(model.m)[[1]]=="lmer")

The models I've passed it, though, have a different 1st element, such that getCall(model.m)[[1]] would be "lme4::lmer" rather than "lmer". Of course, mediate has the technical capacity to fit these models but just doesn't recognize them for what they are.

I think it may be as simple of a fix as changing the if statements to something like

if(isMer.m && class(model.m)[[1]]=="lmerMod")

I am using the mediations() and amelidiate() functions to estimate causal mediated effects with 5 imputed datasets. After running the amelidiate() function on the 5 mediation analyses from the mediations() function, the estimates are equivalent to the mediation based on the first imputed dataset. I did some digging in the source code and realized that the amelidiate() function uses the names from mediations() to organize the data but the names of the 5 mediation analyses are identical so it is only using the first dataset. I am going to rename the resulting list from the mediations() function so that I can use the amelidiate() function. The main application for these functions seems to be a scenario with multiple outcome or mediators and so that resulting names would be different in that case but multiple imputations is also mentioned as an application so I thought this might be a useful comment for future releases of the package or for anyone who is having a similar problem.

I was trying to fit a multilevel mediation analysis with crossed random factors using the mediate function from the mediation package in R. For the med.fit and out.fit models, I used lmer. However, I received the following error message:

"Error in mediate(med.fit, out.fit, treat = "condition_numfct", mediator = "attr_rating", :
mediate does not support more than two levels per model."

I am not sure, but I suspect that this error may be due to the mediate function not supporting crossed random factors. Can anyone confirm this? If so, does anyone know of any alternative packages or functions that can handle multilevel mediation analyses with crossed random factors?

Thank you for any help or insights you can provide!

Running the following model, with MapChange as a binary factor

medModel <- 
  lme4::glmer( MapChange ~ TallCondition + ( 1 | V1.ResponseID ), data=bin.tallSonaK, family=binomial(link="logit") )
outModel <- 
  lme4::lmer( SGAIN ~ MapChange + TallCondition + ( 1 | V1.ResponseID ), data=bin.tallSonaK )
med <-
  mediation::mediate(model.m=medModel, model.y=outModel,treat='TallCondition',mediator='MapChange',data=bin.tallSonaK)
summary(med)

yields the error message Error in if (xhat == 0) out <- 1 else { : missing value where TRUE/FALSE needed

I believe this is an error in function pval:

pval <- function(x, xhat){
  ## Compute p-values
  if (xhat == 0) out <- 1
  else {
    out <- 2 * min(sum(x > 0), sum(x < 0)) / length(x)
  }
  return(min(out, 1))
}

caused by xhat having value NA.

Interestingly I can remove this error by keeping MapChange as num instead of factor

Hi there,

thanks for the package. I managed to use it for a single mediator analysis, but I can't seem to get it working for two mediators, using the mediators function I constantly get this error - Error in rep(1, nrow(dataarg)) : invalid 'times' argument

This is my code :
mediations(datasets = positive, treatment = c("factor_scores"), mediators = c("downward_comparison", "upward_comparison"), outcome = c("self_esteem"), families=c("gaussian","gaussian"), LowerY = NULL, UpperY = NULL, interaction=FALSE, conf.level=.90, sims=100, boot = FALSE, weights=NULL)

There are no NA values in the columns. What am I missing?

Thanks,
Ivaylo

I reduced the code as far as I could and inserted a browser() statement into med.fun.

med.fun has 3 arguments: y.data, m.data and index.

boot requires the function to have an argument "data", but med.fun doesn't.

Therefore, y.data and m.data are not resampled.

Also, new.fit.M <- eval.parent(Call.M)
is problematic, as is the equivalent for Y.

In particular, Call.M and Call.Y reference the original data in the stacked environments. So both models are refit to the original data, not to y.data and m.data.

One approach to fixing this is to add
y.data <- y.data[index, ]
m.data <- m.data[index, ]
near the top of med.fun and to do

model.m <- update(model.m, data = m.data)
model.y <- update(model.y, data = y.data)

instead of re-evaluating Call.M and Call.Y.

Finally, because we're dealing with a treatment and control, I think the observed group sizes should be conditioned on, so the function should default to passing strata = y.data[, treat] into the call to boot (though you can argue I'm wrong about that).

Harry

The formula object in R allows variables to be specified flexibly e.g. y ~ factor(X) or y ~ log(X). Doing this when fitting models for mediate, however, results in error.

This example from the package documentation works as expected:

library("mediation")
data("jobs")
set.seed(4646)

model.m <- lm(job_seek ~ treat + depress1 + econ_hard + sex + age
              + occp + marital + nonwhite + educ + income, data = jobs)
model.y <- lm(depress2 ~ treat + job_seek + depress1 + econ_hard + sex + age
              + occp + marital + nonwhite + educ + income, data = jobs)

out.1 <- mediate(model.m, model.y, sims = 10, boot = TRUE, treat = "treat",
                 mediator = "job_seek")
out.2 <- mediate(model.m, model.y, sims = 10, treat = "treat",
                 mediator = "job_seek")

But would break if we replace sex with factor(sex):

> model.m <- lm(job_seek ~ treat + depress1 + econ_hard + factor(sex) + age
+               + occp + marital + nonwhite + educ + income, data = jobs)
> model.y <- lm(depress2 ~ treat + job_seek + depress1 + econ_hard + sex + age
+               + occp + marital + nonwhite + educ + income, data = jobs)
> 
> out.1 <- mediate(model.m, model.y, sims = 10, boot = TRUE, treat = "treat",
+                  mediator = "job_seek")
Running nonparametric bootstrap

Error in factor(sex) : object 'sex' not found

Or, if we try replacing age with log(age):

> model.m <- lm(job_seek ~ treat + depress1 + econ_hard + sex + log(age)
+               + occp + marital + nonwhite + educ + income, data = jobs)
> model.y <- lm(depress2 ~ treat + job_seek + depress1 + econ_hard + sex + age
+               + occp + marital + nonwhite + educ + income, data = jobs)
> 
> out.1 <- mediate(model.m, model.y, sims = 10, boot = TRUE, treat = "treat",
+                  mediator = "job_seek")
Running nonparametric bootstrap

Error in eval(predvars, data, env) : object 'age' not found

The root cause has to do with how model.frame() handles column names. For example:

> model.m <- lm(job_seek ~ treat + depress1 + econ_hard + factor(sex) + age
+               + occp + marital + nonwhite + educ + income, data = jobs)
> test <- predict(model.frame())
Error in terms.formula(formula, data = data) : 
  argument is not a valid model
> test <- predict(model.m)
> test <- predict(model.m, newdata = model.frame(model.m))
Error in factor(sex) : object 'sex' not found

This works fine with boot=F, but throws an error when boot=T. A real corner case, I admit:

library(mediation)
#> mediation: Causal Mediation Analysis
#> Version: 4.5.0
test <- function(f1, f2, boot){
  mtcars$mpg <- mtcars$mpg/max(mtcars$mpg)
  m_med <- glm(f1, family=binomial, mtcars)
  m_out <- glm(f2, family=binomial, mtcars)
  mediate(m_med, m_out, treat = "wt", mediator = 'am', boot = boot)
}

x <- test('am ~ wt', 'mpg ~ wt + am', F)
#> Warning in eval(family$initialize): non-integer #successes in a binomial glm!
x <- test('am ~ wt', 'mpg ~ wt + am', T)
#> Warning in eval(family$initialize): non-integer #successes in a binomial glm!
#> Running nonparametric bootstrap
#> Error in stats::model.frame(formula = f1, data = structure(list(am = c(1, : object 'f1' not found

^{Created on 2021-03-17 by the reprex package (v0.3.0)}

I have a question regarding the functionality of the medsens function.

I am conducting a causal mediation analysis with a binary outcome and a continous mediator. My data is subject to missing data and I have opted to use multiple imputation in my workflow (through MICE). The mediations package workflow allows me to conduct the mediation analysis on imputed data in a very straightforward manner (through the mediations and amelidiate functions). The summary and plot functions work appropiately; I have obtained the results of the main analysis.

I am now in the stage where I would like to pass the mediate object (that results from the amelidiate function, which I use to pool results from the imputed datasets) through the medsens function in order to conduct the sensitivity analysis. Unfortunately, this is where I run into issues. I receive an error message which mentions that the model combination used in my analysis is not supported. The model combination should however be supported if I am not mistaken (continous mediator, binary outcome. Families in the mediations call are correctly specified as "gaussian" and "binomial", respectively; the mediate object confirms this. No treatment-mediator interaction was specified).

The question: should the medsens function in theory support the workflow described in the vignet accompanying the 'mediation' package on multiple-imputed data? I have no issues running the workflow on the non-imputed data, nor when I run the analysis on the individual imputed datasets (I have not gone through pooling the results from the individual analyses so far).

Looking forward as to any reply. Thanks in advance!