Skip to contents

Relying on a regression-based approach, the exactmed_cat() function calculates standard causal mediation effects when the outcome is binary and the mediator is categorical. More precisely, exactmed_cat() relies on binary and multinomial logistic regression models for the outcome and mediator, respectively, in order to compute exact conditional natural direct and indirect effects. The function returns point and interval estimates for the conditional natural effects without making any assumption regarding the rareness or commonness of the outcome (hence the term exact). For completeness, exactmed_cat() also calculates the conditional controlled direct effect at a specified level of the mediator. Natural and controlled effects estimates are reported using three different scales: odds ratio (OR), risk ratio (RR) and risk difference (RD). The interval estimates can be obtained either by the delta method or the bootstrap.

Usage

exactmed_cat(
  data,
  a,
  m,
  y,
  a1,
  a0,
  m_cov = NULL,
  y_cov = NULL,
  m_cov_cond = NULL,
  y_cov_cond = NULL,
  adjusted = TRUE,
  interaction = TRUE,
  Firth = FALSE,
  boot = FALSE,
  nboot = 1000,
  bootseed = 1991,
  confcoef = 0.95,
  blevel_m = NULL,
  hvalue_y = NULL,
  yprevalence = NULL,
  mf = NULL
)

Arguments

data

A named data frame that includes the exposure, mediator and outcome variables as well as the covariates to be adjusted for in the models. The exposure can be either binary or continuous. If a covariate is categorical, it has to be included in the data frame as a factor, character or logical variable.

a

The name of the binary or continuous exposure variable.

m

The name of the categorical mediator variable.

y

The name of the binary outcome variable.

a1

A value corresponding to the high level of the exposure.

a0

A value corresponding to the low level of the exposure.

m_cov

A vector containing the names of the adjustment variables (covariates) in the mediator model.

y_cov

A vector containing the names of the adjustment variables (covariates) in the outcome model.

m_cov_cond

A named vector (atomic vector or list) containing specific values for some or all of the adjustment covariates m_cov in the mediator model. Please consult the package vignette for details.

y_cov_cond

A named vector (atomic vector or list) containing specific values for some or all of the adjustment covariates y_cov in the outcome model. Please consult the package vignette for details.

adjusted

A logical variable specifying whether to obtain adjusted or unadjusted estimates. If adjusted = FALSE, vectors m_cov and y_cov are ignored by the procedure.

interaction

A logical variable specifying whether there are exposure-mediator interaction terms in the outcome model.

Firth

A logical variable specifying whether to compute conventional or penalized maximum likelihood estimates for the outcome logistic regression model.

boot

A logical value specifying whether the confidence intervals are obtained by the delta method or by percentile bootstrap.

nboot

The number of bootstrap replications used to obtain the confidence intervals if boot = TRUE. By default nboot = 1000.

bootseed

The value of the initial seed (positive integer) for random number generation if boot = TRUE.

confcoef

A number between 0 and 1 for the confidence coefficient (ex.: 0.95) of the interval estimates.

blevel_m

The reference level of the mediator. If it is not specified blevel_m is fixed to the first level of the mediator (default).

hvalue_y

The value corresponding to the high level of the outcome. If the outcome is already coded as a numerical binary variable taking 0 or 1 values, then by default hvalue_y = 1.

yprevalence

The prevalence of the outcome in the population (a number between 0 and 1). Option used when case-control data are used. The low level of the outcome is treated as the control level.

mf

The level of the mediator at which the conditional controlled direct effect is computed. If it is not specified, mf is fixed at the reference level (blevel_m) of the mediator (default).

Value

An object of class results_cat is returned:

ne.or

Natural effects estimates on OR scale.

ne.rr

Natural effects estimates on RR scale.

ne.rd

Natural effects estimates on RD scale.

cde

Controlled direct effect estimates.

med.reg

Summary of the mediator regression.

out.reg

Summary of the outcome regression.

If boot==TRUE, the returned object also contains:

boot.ne.or

Bootstrap replications of natural effects on OR scale.

boot.ne.rr

Bootstrap replications of natural effects on RR scale.

boot.ne.rd

Bootstrap replications of natural effects on RD scale.

boot.cde.or

Bootstrap replications of controlled direct effect on OR scale.

boot.cde.rr

Bootstrap replications of controlled direct effect on RR scale.

boot.cde.rd

Bootstrap replications of controlled direct effect on RD scale.

boot.ind

Indices of the observations sampled in each bootstrap replication (one replication per column).

Details

By default, exactmed_cat() reports mediation effects evaluated at the sample-specific mean values of the numerical covariates (including the dummy variables created internally by the function to represent the non-reference levels of the categorical covariates). In order to estimate mediation effects at specific values of some covariates (that is, stratum-specific effects), the user needs to provide named vectors m_cov_cond and/or y_cov_cond containing those values or levels. The adjustment covariates appearing in both m_cov and y_cov (common adjustment covariates) must have the same values; otherwise, exactmed_cat()'s execution is aborted and an error message is displayed in the R console.

The Firth parameter allows to reduce the bias of the outcome logistic regression coefficients estimators when facing a problem of separation or quasi-separation. The bias reduction is achieved by the brglmFit fitting method of the brglm2 package. More precisely, estimates are obtained via penalized maximum likelihood with a Jeffreys prior penalty, which is equivalent to the mean bias-reducing adjusted score equation approach in Firth (1993).

When the data come from a case-control study, the yprevalence parameter should be used and its value ideally correspond to the true outcome prevalence. exactmed_cat() accounts for the ascertainment in the sample by employing weighted regression techniques that use inverse probability weighting (IPW) with robust standard errors. These errors are obtained via the vcovHC function of the R package sandwich. Specifically, we use the HC3 type covariance matrix estimator (default type of the vcovHC function).

For the mediation effects expressed on the multiplicative scales (odds ratio, OR; risk ratio, RR), the exactmed_cat() function returns delta method confidence intervals by exponentiating the lower and upper limits of the normal confidence intervals obtained for the logarithmic transformations of the effects. The exactmed_cat() function also provides the estimated standard errors of natural and controlled direct effects estimators that are not log-transformed, where those are derived using a first order Taylor expansion (e.g., \(\hat{SE}(\hat{OR})=\hat{OR} \times \hat{SE}(\log(\hat{OR}))\)). The function performs Z-tests (null hypothesis: there is no effect) computing the corresponding two-tailed p-values. Note that for the multiplicative scales, the standard scores (test statistics) are obtained by dividing the logarithm of an effect estimator by the estimator of the corresponding standard error (e.g., \(\log(\hat{OR}) / \hat{SE}(\log(\hat{OR}))\)). No log-transformation is applied when working on the risk difference scale.

Note

The exactmed_cat() function only works for complete data. Users can apply multiple imputation techniques (e.g., R package mice) or remove observations of variables used in mediation analysis that have missing values (NA).

References

Samoilenko M, Blais L, Lefebvre G. Comparing logistic and log-binomial models for causal mediation analyses of binary mediators and rare binary outcomes: evidence to support cross-checking of mediation results in practice. Observational Studies.2018;4(1):193-216.

Samoilenko M, Lefebvre G. Parametric-regression-based causal mediation analysis of binary outcomes and binary mediators: moving beyond the rareness or commonness of the outcome. American Journal of Epidemiology.2021;190(9):1846-1858. doi:10.1093/aje/kwab055 .

Samoilenko M, Lefebvre G. An exact regression-based approach for the estimation of the natural direct and indirect effects with a binary outcome and a continuous mediator. Statistics in Medicine.2023; 42(3): 353–387. doi:10.1002/sim.9621 .

Firth D. Bias reduction of maximum likelihood estimates. Biometrika.1993;80:27-38. doi:10.2307/2336755 .

Examples

exactmed_cat(
  data = datamed_cat, a = "X", m = "M", y = "Y", a1 = 1, a0 = 0,
  m_cov = c("C1", "C2"), y_cov = c("C1", "C2")
)
#> 
#> Natural effects on OR scale: 
#> 
#>                 Estimate Std.error    2.5%   97.5%      P.val
#> Direct effect    4.41076   0.77459 3.12629 6.22295 < 2.22e-16
#> Indirect effect  0.75882   0.06912 0.63476 0.90713  0.0024451
#> Total effect     3.34697   0.49821 2.50004 4.48081 4.4409e-16
#> 
#> Natural effects on RR scale: 
#> 
#>                 Estimate Std.error    2.5%   97.5%      P.val
#> Direct effect    1.85088   0.13312 1.60752 2.13108 < 2.22e-16
#> Indirect effect  0.92653   0.02233 0.88379 0.97135  0.0015453
#> Total effect     1.71491   0.11852 1.49765 1.96367 5.9952e-15
#> 
#> Natural effects on RD scale: 
#> 
#>                 Estimate Std.error     2.5%    97.5%      P.val
#> Direct effect    0.34503   0.03670  0.27309  0.41697 < 2.22e-16
#> Indirect effect -0.05514   0.01754 -0.08951 -0.02076  0.0016682
#> Total effect     0.28989   0.03365  0.22393  0.35585 < 2.22e-16
#> 
#> Controlled direct effect (m='a'): 
#> 
#>          Estimate Std.error     2.5%     97.5%      P.val
#> OR scale 45.28588  23.96135 16.05391 127.74523 5.7465e-13
#> RR scale  7.19900   1.56042  4.70728  11.00966 < 2.22e-16
#> RD scale  0.74056   0.06235  0.61836   0.86275 < 2.22e-16
#> 
#> 
#> Mediator model: 
#> 
#> 
#> Call:
#> mlogit(formula = M ~ 1 | X + C1 + C2, data = data_long, reflevel = ref_level, 
#>     method = "nr")
#> 
#> Frequencies of alternatives:choice
#>   a   b   c   d   e 
#> 0.2 0.2 0.2 0.2 0.2 
#> 
#> nr method
#> 6 iterations, 0h:0m:0s 
#> g'(-H)^-1g = 2.06E-05 
#> successive function values within tolerance limits 
#> 
#> Coefficients :
#>                Estimate Std. Error z-value  Pr(>|z|)    
#> (Intercept):b  0.539165   0.202979  2.6563  0.007901 ** 
#> (Intercept):c  0.622236   0.209098  2.9758  0.002922 ** 
#> (Intercept):d  0.435369   0.217929  1.9978  0.045743 *  
#> (Intercept):e -0.089258   0.244928 -0.3644  0.715541    
#> X:b            1.024819   0.228423  4.4865 7.241e-06 ***
#> X:c            1.688852   0.247460  6.8248 8.808e-12 ***
#> X:d            2.222046   0.259408  8.5658 < 2.2e-16 ***
#> X:e            2.750817   0.286159  9.6129 < 2.2e-16 ***
#> C1:b          -0.654903   0.220269 -2.9732  0.002947 ** 
#> C1:c          -1.167718   0.237979 -4.9068 9.257e-07 ***
#> C1:d          -1.537355   0.248169 -6.1948 5.836e-10 ***
#> C1:e          -2.271721   0.276774 -8.2079 2.220e-16 ***
#> C2:b           0.828863   0.135509  6.1167 9.555e-10 ***
#> C2:c           1.662424   0.156101 10.6497 < 2.2e-16 ***
#> C2:d           1.997098   0.164942 12.1079 < 2.2e-16 ***
#> C2:e           2.934705   0.188921 15.5341 < 2.2e-16 ***
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> Log-Likelihood: -1325.6
#> McFadden R^2:  0.17635 
#> Likelihood ratio test : chisq = 567.66 (p.value = < 2.22e-16)
#> 
#> Outcome model: 
#> 
#> 
#> Call:
#> glm(formula = Yform, family = binomial(link = "logit"), data = data)
#> 
#> Coefficients:
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept) -2.67124    0.27264  -9.798  < 2e-16 ***
#> X            3.81300    0.52911   7.206 5.75e-13 ***
#> Mb           1.06895    0.31187   3.428 0.000609 ***
#> Mc           2.19709    0.33371   6.584 4.59e-11 ***
#> Md           2.78578    0.36477   7.637 2.22e-14 ***
#> Me           3.37526    0.40984   8.236  < 2e-16 ***
#> X:Mb        -1.69550    0.63866  -2.655 0.007936 ** 
#> X:Mc        -3.00850    0.61553  -4.888 1.02e-06 ***
#> X:Md        -4.11908    0.61667  -6.680 2.40e-11 ***
#> X:Me        -4.49990    0.62393  -7.212 5.50e-13 ***
#> C1           1.34154    0.16052   8.358  < 2e-16 ***
#> C2          -0.48389    0.09375  -5.162 2.45e-07 ***
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> (Dispersion parameter for binomial family taken to be 1)
#> 
#>     Null deviance: 1377.4  on 999  degrees of freedom
#> Residual deviance: 1112.0  on 988  degrees of freedom
#> AIC: 1136
#> 
#> Number of Fisher Scoring iterations: 5
#> 

exactmed_cat(
  data = datamed_cat, a = "X", m = "M", y = "Y", a1 = 1, a0 = 0,
  m_cov = c("C1", "C2"), y_cov = c("C1", "C2"), yprevalence = 0.1
)
#> 
#> Natural effects on OR scale: 
#> 
#>                 Estimate Std.error    2.5%   97.5%      P.val
#> Direct effect    5.10315   1.60423 2.75584 9.44980 2.1639e-07
#> Indirect effect  0.66606   0.15266 0.42503 1.04379    0.07623
#> Total effect     3.39902   0.63694 2.35422 4.90750 6.6158e-11
#> 
#> Natural effects on RR scale: 
#> 
#>                 Estimate Std.error    2.5%   97.5%      P.val
#> Direct effect    4.07813   1.02846 2.48770 6.68536 2.4930e-08
#> Indirect effect  0.72668   0.12473 0.51910 1.01729   0.062874
#> Total effect     2.96352   0.49660 2.13389 4.11570 8.9856e-11
#> 
#> Natural effects on RD scale: 
#> 
#>                 Estimate Std.error     2.5%   97.5%      P.val
#> Direct effect    0.18856   0.05096  0.08868 0.28843 0.00021528
#> Indirect effect -0.06828   0.04347 -0.15348 0.01692 0.11626210
#> Total effect     0.12028   0.01930  0.08246 0.15810 4.5611e-10
#> 
#> Controlled direct effect (m='a'): 
#> 
#>          Estimate Std.error     2.5%     97.5%      P.val
#> OR scale 50.94107  32.31878 14.69049 176.64440 5.8083e-10
#> RR scale 34.03192  15.81876 13.68453  84.63366 3.2419e-14
#> RD scale  0.32863   0.12981  0.07422   0.58305    0.01135
#> 
#> 
#> Mediator model: 
#> 
#> 
#> t test of coefficients:
#> 
#>               Estimate Std. Error t value  Pr(>|t|)    
#> (Intercept):b  0.46021    0.28632  1.6073  0.108304    
#> (Intercept):c  0.23238    0.30606  0.7593  0.447878    
#> (Intercept):d -0.17027    0.33552 -0.5075  0.611933    
#> (Intercept):e -0.94990    0.37577 -2.5279  0.011630 *  
#> X:b            1.36477    0.28420  4.8022 1.813e-06 ***
#> X:c            2.55073    0.33390  7.6392 5.172e-14 ***
#> X:d            3.57849    0.38840  9.2135 < 2.2e-16 ***
#> X:e            4.14570    0.43632  9.5015 < 2.2e-16 ***
#> C1:b          -0.87770    0.30533 -2.8746  0.004132 ** 
#> C1:c          -1.57312    0.33565 -4.6867 3.167e-06 ***
#> C1:d          -2.13493    0.34608 -6.1690 1.003e-09 ***
#> C1:e          -2.80812    0.39506 -7.1080 2.260e-12 ***
#> C2:b           0.74211    0.18065  4.1080 4.323e-05 ***
#> C2:c           1.83417    0.20581  8.9118 < 2.2e-16 ***
#> C2:d           2.14591    0.22479  9.5461 < 2.2e-16 ***
#> C2:e           3.14727    0.25797 12.2000 < 2.2e-16 ***
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> 
#> Outcome model: 
#> 
#> 
#> z test of coefficients:
#> 
#>             Estimate Std. Error  z value  Pr(>|z|)    
#> (Intercept) -5.30443    0.32627 -16.2579 < 2.2e-16 ***
#> X            3.93067    0.63443   6.1955 5.808e-10 ***
#> Mb           1.13266    0.33207   3.4109 0.0006475 ***
#> Mc           2.53737    0.38554   6.5813 4.664e-11 ***
#> Md           3.14575    0.42668   7.3726 1.673e-13 ***
#> Me           3.72830    0.50238   7.4212 1.160e-13 ***
#> X:Mb        -1.70977    0.74079  -2.3080 0.0209973 *  
#> X:Mc        -3.24768    0.72475  -4.4811 7.426e-06 ***
#> X:Md        -4.26365    0.72327  -5.8950 3.747e-09 ***
#> X:Me        -4.67521    0.75025  -6.2315 4.619e-10 ***
#> C1           1.39184    0.18345   7.5868 3.278e-14 ***
#> C2          -0.61007    0.11767  -5.1848 2.162e-07 ***
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 

m_cov_cond <- c(C1 = 0.1, C2 = 0.4)
y_cov_cond <- c(C1 = 0.1, C2 = 0.4)

exactmed_cat(
  data = datamed_cat, a = "X", m = "M", y = "Y", a1 = 1, a0 = 0,
  m_cov = c("C1", "C2"), y_cov = c("C1", "C2"),
  m_cov_cond = m_cov_cond, y_cov_cond = y_cov_cond
)
#> 
#> Natural effects on OR scale: 
#> 
#>                 Estimate Std.error    2.5%   97.5%      P.val
#> Direct effect    1.80534   0.30422 1.29755 2.51186 0.00045536
#> Indirect effect  0.84292   0.05871 0.73536 0.96621 0.01414837
#> Total effect     1.52175   0.24283 1.11305 2.08052 0.00850960
#> 
#> Natural effects on RR scale: 
#> 
#>                 Estimate Std.error    2.5%   97.5%      P.val
#> Direct effect    1.36805   0.12324 1.14663 1.63222 0.00050347
#> Indirect effect  0.92152   0.03049 0.86365 0.98326 0.01350957
#> Total effect     1.26068   0.11211 1.05904 1.50071 0.00918811
#> 
#> Natural effects on RD scale: 
#> 
#>                 Estimate Std.error     2.5%    97.5%     P.val
#> Direct effect    0.14608   0.04110  0.06553  0.22663 0.0003785
#> Indirect effect -0.04262   0.01732 -0.07657 -0.00866 0.0138989
#> Total effect     0.10347   0.03901  0.02701  0.17992 0.0079922
#> 
#> Controlled direct effect (m='a'): 
#> 
#>          Estimate Std.error     2.5%     97.5%      P.val
#> OR scale 45.28588  23.96135 16.05391 127.74523 5.7465e-13
#> RR scale 12.20703   3.32817  7.15378  20.82978 < 2.22e-16
#> RD scale  0.68575   0.09315  0.50317   0.86833 1.8185e-13
#> 
#> 
#> Mediator model: 
#> 
#> 
#> Call:
#> mlogit(formula = M ~ 1 | X + C1 + C2, data = data_long, reflevel = ref_level, 
#>     method = "nr")
#> 
#> Frequencies of alternatives:choice
#>   a   b   c   d   e 
#> 0.2 0.2 0.2 0.2 0.2 
#> 
#> nr method
#> 6 iterations, 0h:0m:0s 
#> g'(-H)^-1g = 2.06E-05 
#> successive function values within tolerance limits 
#> 
#> Coefficients :
#>                Estimate Std. Error z-value  Pr(>|z|)    
#> (Intercept):b  0.539165   0.202979  2.6563  0.007901 ** 
#> (Intercept):c  0.622236   0.209098  2.9758  0.002922 ** 
#> (Intercept):d  0.435369   0.217929  1.9978  0.045743 *  
#> (Intercept):e -0.089258   0.244928 -0.3644  0.715541    
#> X:b            1.024819   0.228423  4.4865 7.241e-06 ***
#> X:c            1.688852   0.247460  6.8248 8.808e-12 ***
#> X:d            2.222046   0.259408  8.5658 < 2.2e-16 ***
#> X:e            2.750817   0.286159  9.6129 < 2.2e-16 ***
#> C1:b          -0.654903   0.220269 -2.9732  0.002947 ** 
#> C1:c          -1.167718   0.237979 -4.9068 9.257e-07 ***
#> C1:d          -1.537355   0.248169 -6.1948 5.836e-10 ***
#> C1:e          -2.271721   0.276774 -8.2079 2.220e-16 ***
#> C2:b           0.828863   0.135509  6.1167 9.555e-10 ***
#> C2:c           1.662424   0.156101 10.6497 < 2.2e-16 ***
#> C2:d           1.997098   0.164942 12.1079 < 2.2e-16 ***
#> C2:e           2.934705   0.188921 15.5341 < 2.2e-16 ***
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> Log-Likelihood: -1325.6
#> McFadden R^2:  0.17635 
#> Likelihood ratio test : chisq = 567.66 (p.value = < 2.22e-16)
#> 
#> Outcome model: 
#> 
#> 
#> Call:
#> glm(formula = Yform, family = binomial(link = "logit"), data = data)
#> 
#> Coefficients:
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept) -2.67124    0.27264  -9.798  < 2e-16 ***
#> X            3.81300    0.52911   7.206 5.75e-13 ***
#> Mb           1.06895    0.31187   3.428 0.000609 ***
#> Mc           2.19709    0.33371   6.584 4.59e-11 ***
#> Md           2.78578    0.36477   7.637 2.22e-14 ***
#> Me           3.37526    0.40984   8.236  < 2e-16 ***
#> X:Mb        -1.69550    0.63866  -2.655 0.007936 ** 
#> X:Mc        -3.00850    0.61553  -4.888 1.02e-06 ***
#> X:Md        -4.11908    0.61667  -6.680 2.40e-11 ***
#> X:Me        -4.49990    0.62393  -7.212 5.50e-13 ***
#> C1           1.34154    0.16052   8.358  < 2e-16 ***
#> C2          -0.48389    0.09375  -5.162 2.45e-07 ***
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> (Dispersion parameter for binomial family taken to be 1)
#> 
#>     Null deviance: 1377.4  on 999  degrees of freedom
#> Residual deviance: 1112.0  on 988  degrees of freedom
#> AIC: 1136
#> 
#> Number of Fisher Scoring iterations: 5
#> 

C1b <- factor(sample(c("a", "b", "c"), nrow(datamed_cat), replace = TRUE))
datamed_cat$C1 <- C1b

m_cov_cond <- list(C1 = "c", C2 = 0.4)
y_cov_cond <- list(C1 = "c", C2 = 0.4)

exactmed_cat(
  data = datamed_cat, a = "X", m = "M", y = "Y", a1 = 1, a0 = 0,
  m_cov = c("C1", "C2"), y_cov = c("C1", "C2"),
  m_cov_cond = m_cov_cond, y_cov_cond = y_cov_cond
)
#> 
#> Natural effects on OR scale: 
#> 
#>                 Estimate Std.error    2.5%   97.5%      P.val
#> Direct effect    2.79709   0.47413 2.00642 3.89934 1.2945e-09
#> Indirect effect  0.73321   0.05085 0.64003 0.83996 7.6501e-06
#> Total effect     2.05087   0.31705 1.51477 2.77669 3.3819e-06
#> 
#> Natural effects on RR scale: 
#> 
#>                 Estimate Std.error    2.5%   97.5%      P.val
#> Direct effect    1.51882   0.10932 1.31897 1.74894 6.3874e-09
#> Indirect effect  0.90494   0.02139 0.86398 0.94784 2.3732e-05
#> Total effect     1.37444   0.09639 1.19793 1.57696 5.7588e-06
#> 
#> Natural effects on RD scale: 
#> 
#>                 Estimate Std.error     2.5%    97.5%      P.val
#> Direct effect    0.24298   0.03824  0.16802  0.31793 2.1079e-10
#> Indirect effect -0.06762   0.01517 -0.09735 -0.03788 8.3164e-06
#> Total effect     0.17536   0.03686  0.10311  0.24761 1.9648e-06
#> 
#> Controlled direct effect (m='a'): 
#> 
#>          Estimate Std.error     2.5%     97.5%      P.val
#> OR scale 38.82017  19.92499 14.19601 106.15697 1.0127e-12
#> RR scale  5.27108   1.09818  3.50397   7.92937 1.5543e-15
#> RD scale  0.71878   0.05547  0.61005   0.82751 < 2.22e-16
#> 
#> 
#> Mediator model: 
#> 
#> 
#> Call:
#> mlogit(formula = M ~ 1 | X + C1 + C2, data = data_long, reflevel = ref_level, 
#>     method = "nr")
#> 
#> Frequencies of alternatives:choice
#>   a   b   c   d   e 
#> 0.2 0.2 0.2 0.2 0.2 
#> 
#> nr method
#> 6 iterations, 0h:0m:0s 
#> g'(-H)^-1g = 1.46E-06 
#> successive function values within tolerance limits 
#> 
#> Coefficients :
#>                 Estimate Std. Error z-value  Pr(>|z|)    
#> (Intercept):b -0.0929394  0.2221877 -0.4183 0.6757335    
#> (Intercept):c  0.0148042  0.2261272  0.0655 0.9478011    
#> (Intercept):d -0.3404610  0.2403031 -1.4168 0.1565419    
#> (Intercept):e -1.0307691  0.2735875 -3.7676 0.0001648 ***
#> X:b            0.9707211  0.2258020  4.2990 1.716e-05 ***
#> X:c            1.5464564  0.2411365  6.4132 1.425e-10 ***
#> X:d            2.0510189  0.2509649  8.1725 2.220e-16 ***
#> X:e            2.5224051  0.2740241  9.2050 < 2.2e-16 ***
#> C1b:b          0.3801689  0.2599386  1.4625 0.1435951    
#> C1b:c          0.0011054  0.2718014  0.0041 0.9967551    
#> C1b:d         -0.0588294  0.2827345 -0.2081 0.8351720    
#> C1b:e         -0.0565024  0.3032325 -0.1863 0.8521832    
#> C1c:b          0.2308921  0.2619274  0.8815 0.3780410    
#> C1c:c         -0.1498283  0.2752639 -0.5443 0.5862297    
#> C1c:d         -0.0186318  0.2816653 -0.0661 0.9472593    
#> C1c:e         -0.1318734  0.3061992 -0.4307 0.6667020    
#> C2:b           0.7564177  0.1303864  5.8014 6.578e-09 ***
#> C2:c           1.5371809  0.1486267 10.3426 < 2.2e-16 ***
#> C2:d           1.8367214  0.1563525 11.7473 < 2.2e-16 ***
#> C2:e           2.7075166  0.1772995 15.2709 < 2.2e-16 ***
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> Log-Likelihood: -1363.8
#> McFadden R^2:  0.15261 
#> Likelihood ratio test : chisq = 491.25 (p.value = < 2.22e-16)
#> 
#> Outcome model: 
#> 
#> 
#> Call:
#> glm(formula = Yform, family = binomial(link = "logit"), data = data)
#> 
#> Coefficients:
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept) -1.67108    0.24653  -6.778 1.22e-11 ***
#> X            3.65894    0.51326   7.129 1.01e-12 ***
#> Mb           0.72431    0.29721   2.437 0.014808 *  
#> Mc           1.73050    0.31306   5.528 3.25e-08 ***
#> Md           2.19788    0.33835   6.496 8.26e-11 ***
#> Me           2.41026    0.37002   6.514 7.32e-11 ***
#> X:Mb        -1.42592    0.61862  -2.305 0.021167 *  
#> X:Mc        -2.88105    0.59566  -4.837 1.32e-06 ***
#> X:Md        -3.95284    0.59535  -6.639 3.15e-11 ***
#> X:Me        -4.11068    0.59951  -6.857 7.05e-12 ***
#> C1b          0.03451    0.17089   0.202 0.839959    
#> C1c          0.20335    0.17342   1.173 0.240982    
#> C2          -0.32515    0.08795  -3.697 0.000218 ***
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> (Dispersion parameter for binomial family taken to be 1)
#> 
#>     Null deviance: 1377.4  on 999  degrees of freedom
#> Residual deviance: 1187.0  on 987  degrees of freedom
#> AIC: 1213
#> 
#> Number of Fisher Scoring iterations: 4
#>