Running a Functional Mediation Model
Now we call the functional mediation function. Only 10 bootstrap samples are used below for this quick illustration. At least a few hundred bootstraps are recommended in practice in order to increase precision and power.
model1 <- funmediation(data=the_data,
treatment=X,
mediator=M,
outcome=Y,
id=subject_id,
time=t,
nboot=10)
#> Warning in funmediation(data = the_data, treatment = X, mediator = M, outcome =
#> Y, : The number of bootstrap samples was very small and results will be
#> dubious.
#> Warning in norm.inter(t, (1 + c(conf, -conf))/2): extreme order statistics used
#> as endpoints
#> Warning in norm.inter(t, alpha): extreme order statistics used as endpoints
#> Warning in norm.inter(t, (1 + c(conf, -conf))/2): extreme order statistics used
#> as endpoints
#> Warning in norm.inter(t, alpha): extreme order statistics used as endpointsThe warning message, ``extreme order statistics used as endpoints,’’ occurs because too few bootstrap samples are used for the bootstrap confidence intervals to be valid; however, we ignore this here in order for the vignette example to run quickly.
The print function gives an initial overview of the
results.
print(model1)
#> =======================================================
#> Functional Regression Mediation Function Output
#> =======================================================
#> TREATMENT: X
#> MEDIATOR: M (Assumed Continuous)
#> OUTCOME: Y (Assumed Continuous)
#> =======================================================
#> DIRECT EFFECT OF X ...
#> Direct effect estimate:
#> -0.07192656
#> Direct effect standard error:
#> 0.1047166
#> Direct effect bootstrap estimate:
#> -0.1102083
#> Direct effect bootstrap confidence interval:
#> ... by normal method:
#> -0.2939 , 0.0735
#> =======================================================
#> INDIRECT EFFECT OF X ...
#> Indirect effect parametric estimate:
#> -0.2522696
#> Indirect effect bootstrap estimate:
#> -0.2460247
#> Indirect effect bootstrap confidence interval:
#> ... by normal method:
#> -0.3517 , -0.1404
#> =======================================================
#> Computation time:
#> Time difference of 18.58889 secs
#> =======================================================
#> Time-varying Effects Model Predicting MEDIATOR from TREATMENT:
#> Response variable: MEDIATOR
#> Response outcome distribution type: gaussian
#> Time interval: 0.01 to 1
#> Number of subjects: 500
#> Effects specified as time-varying: (Intercept), TREATMENT
#> You can use the plot function to view the plots.
#>
#> Back-end model fitted in mgcv::bam function:
#> Method fREML
#> Formula:
#> MEDIATOR ~ TREATMENT + s(time, bs = "ps", by = NA, pc = 0, k = 7,
#> fx = FALSE) + s(time, bs = "ps", by = TREATMENT, pc = 0,
#> m = c(2, 1), k = 7, fx = FALSE)
#> Pseudolikelihood AIC: 84897.63
#> Pseudolikelihood BIC: 84937.51
#>
#> =======================================================
#> Functional Regression Model Predicting OUTCOME from MEDIATOR
#> and TREATMENT:
#>
#> Family: gaussian
#> Link function: identity
#>
#> Formula:
#> OUTCOME ~ s(x = MEDIATOR.tmat, by = L.MEDIATOR) + TREATMENT
#>
#> Estimated degrees of freedom:
#> 2 total = 4
#>
#> REML score: 742.2459
#> Scalar terms in functional regression model:
#> (Intercept) TREATMENT
#> 0.16499534 -0.07192656
#> =======================================================
#> Parametric model for Predicting OUTCOME from TREATMENT
#> ignoring MEDIATOR:
#>
#> Call:
#> glm(formula = glm_formula)
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 0.47332 0.06671 7.096 4.45e-12 ***
#> TREATMENT -0.30617 0.09730 -3.147 0.00175 **
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> (Dispersion parameter for gaussian family taken to be 1.179171)
#>
#> Null deviance: 598.90 on 499 degrees of freedom
#> Residual deviance: 587.23 on 498 degrees of freedom
#> AIC: 1505.3
#>
#> Number of Fisher Scoring iterations: 2
#>
#> =======================================================The first part of the output above summarizes the estimated indirect
effect and its bootstrap confidence interval. The estimate for the
indirect effect \(\int
\alpha_X(t)\beta_M(t)\) is given as about -0248. Choosing the
wider confidence interval in order to be conservative, a 95% confidence
interval extends from about -0.35 to -0.14. (This interval would not
actually be reliable in practice because of the low number of bootstrap
samples, which could be remedied by choosing a higher
nboot). Thus, there appears to be a significant negative
indirect effect of \(X\) on \(Y\) mediated through \(M(t)\).
The second part of the printed output summarizes the TVEM used to
predict \(M(t)\) from \(X\). It only provides summary information
on the model that was fit. The time-varying coefficients are not
actually printed, because they are functions rather than single numbers.
They can be plotted using the plot function, as described
later. They are also stored in the model1 output object, as
described later.
The third part summarizes the scalar-on-function functional
regression model used to predict \(Y\)
from \(X\) and \(M(t)\). This model involves scalar
(non-time-varying) coefficients for the intercept and the direct effect
of \(X\), and a functional
(time-varying) coefficient for the effect of \(M(t)\). The scalar coefficients are
printed, but as before, the functional coefficient needs to be plotted
using the plot function.
The fourth section of the printed output simply presents an estimate
of the total effect of \(X\) on \(Y\) ignoring \(M(t)\), obtained using the glm
function. \(X\) is shown to have a
total negative effect on \(Y\), with an
estimated coefficient of -0.30617 and \(p\)-value of 0.00175. This
suggests that the \(X=1\) group has
statistically significantly lower average \(Y\) than the \(X=0\) group, irrespective of \(M\).
The plot function plots the results. Three kinds of
plots are available, represented by the options tvem,
pfr, and coef. The tvem plot
shows the estimates for \(\alpha_0(t)\)
and \(\alpha_1(t)\).
The intercept plot appears to be significantly above zero and significantly increasing. It can be interpreted as the mean \(M(t)\) for the control (\(X=0\)) group because \(E(M(t)|X=0)=\alpha_0(t) + \alpha_X(t) \times 0=\alpha_0(t)\). The treatment effect plot is generally nonzero and decreasing. It can be interpreted as \(E(M(t)|X=1)-E(M(t)|X=0)\), a time-specific treatment effect on the mediator, and it suggests that as time goes on the \(X=1\) group tends to become lower in average \(M(t)\) than the \(X=0\) group.
The pfr plot shows the estimate for \(\beta_M(t)\).
The function tends to be nonzero at least for higher values of time, and generally seems to be steadily increasing. One possible interpretation is that the mediator becomes more important at later times. Alternatively, perhaps individuals who increase more (or decrease less) on average over time on the mediator tend to have a higher value on the outcome (see Dziak et al., 2019).
The pfrgam plot does the same thing as the
pfr plot, but with a slightly different implementation
based on the plot.gam method from the mgcv
library.
Last, the coef plot summarizes all of the most important
coefficients estimated in the model. The upper left and upper right
panes show the estimates of \(\alpha_0(t)\) and \(\alpha_X(t)\) from the model of the effect
of \(X\) on \(M(t)\). The lower left shows the estimate
of \(\beta_M(t)\) from the model of the
effect of \(M(t)\) on \(Y\) given \(X\). In the lower right panel, the estimate
of the indirect effect is printed (not plotted, because it is a single
number).
Because we have already looked at the relevant plots before using the
coef option, there is nothing new to interpret in the
coef plot; it is just intended as a convenient summary.
Exploring Information in the Output Objects
The output object created, called model1, contains
useful information both on the model fit to the original data and on the
bootstrapping results.
model1$original_results$time_gridis the grid of time points on which the coefficient functions were estimated.model1$original_results$alpha_int_estimate,model1$original_results$alpha_int_se,model1$original_results$alpha_X_estimate, andmodel1$original_results$alpha_X_seare the functional coefficients and pointwise standard errors for the intercept and treatment in the TVEM, at each point of the time grid.model1$original_results$beta_int_estimate,model1$original_results$beta_int_se,model1$original_results$beta_X_estimateandmodel1$original_results$beta_X_seare the estimates and standard errors for the scalar coefficients representing the intercept and the direct treatment effect in the scalar-outcome functional linear model.model1$original_results$beta_M_estimateandmodel1$original_results$beta_M_seare the estimates and pointwise standard errors representing the effect of the time-varying mediator at each point of the time grid.beta_M_pvalueis the \(p\)-value for an overall test of whether the mediator is significantly associated with the outcome, i.e., whether \(\beta_M(t)\) can be shown not to be all zeroes; see the documentation for thepfrfunction in therefundpackage for more information.model1$original_results$tau_int_estimate, together with$tau_int_se,$tau_X_estimate, and$tau_X_se, represent the estimates of the intercept and coefficient in the model for the total effect of \(X\) on \(Y\). Thus,model1$original_results$tau_X_estimateis the estimated total effect andmodel1$original_results$tau_X_pvalueis the \(p\)-value for testing whether this effect is zero.model1$original_results$indirect_effect_estimategives the estimated value of \(\int \alpha_X(t)\beta_M(t)\) obtained without bootstrapping.model1$original_results$tvem_XM_details,$funreg_MY_details, and$total_effect_detailsgive the fitted model objects fromtvem,pfrandglm. They might be useful for plots or diagnostics but might contain more technical detail than some users are interested in.
Including a Covariate
It is possible to run the model with one or more additional
covariates. To demonstrate this, the
simulate_mediation_example has an option
make_covariate_S to create an extra subject-level
covariate. The data is generated in such a way that \(S\) is not actually related to \(X\), \(M\)
or \(Y\), but the analyst is assumed
not to know this – that is, the simulated true value of its coefficient
is zero.
simulation_with_covariate <- simulate_funmediation_example(nsub=500,
make_covariate_S=TRUE);
data_with_covariate <- simulation_with_covariate$dataset; To include the covariate \(S_i\)
with an assumption of a possibly time-varying effect \(\alpha_S(t)\) on \(M(t)\) (recall that the covariate can have
a time-varying effect even though it may not have time-varying values),
use the covariates_on_outcome and
tve_covariates_on_mediator arguments:
model_with_tve_covariate <- funmediation(data=data_with_covariate,
treatment=X,
mediator=M,
outcome=Y,
tve_covariates_on_mediator = ~S,
covariates_on_outcome = ~S,
id=subject_id,
time=t,
nboot=10)
#> Warning in funmediation(data = data_with_covariate, treatment = X, mediator =
#> M, : The number of bootstrap samples was very small and results will be
#> dubious.
#> Warning in norm.inter(t, (1 + c(conf, -conf))/2): extreme order statistics used
#> as endpoints
#> Warning in norm.inter(t, alpha): extreme order statistics used as endpoints
#> Warning in norm.inter(t, (1 + c(conf, -conf))/2): extreme order statistics used
#> as endpoints
#> Warning in norm.inter(t, alpha): extreme order statistics used as endpointsNotice that the tilde (~) character needs to be included for the
covariates. This is done so that R will allow multiple covariates to be
listed in the style of the right side of a formula, e.g.,
~S1 + S2.
The covariates in tve_covariates_on_mediator and
covariates_on_outcome need to be subject-level and should
not have time-varying values; the function does not currently support
models with time-varying values other than the mediator. However, the
relationship of the covariates to the mediator can be specified as
either time-varying or time-invariant. The relationship of the covariate
to the outcome, in contrast, can only be specified as time- invariant in
the current version of the package. That is, \(M(t)\) is the only time-varying predictor
in the functional regression predicting the outcome.
The results of the model fit can be viewed using the
print and plot functions:
print(model_with_tve_covariate)
#> =======================================================
#> Functional Regression Mediation Function Output
#> =======================================================
#> TREATMENT: X
#> MEDIATOR: M (Assumed Continuous)
#> OUTCOME: Y (Assumed Continuous)
#> =======================================================
#> DIRECT EFFECT OF X ...
#> Direct effect estimate:
#> 0.1125763
#> Direct effect standard error:
#> 0.1009101
#> Direct effect bootstrap estimate:
#> 0.06677104
#> Direct effect bootstrap confidence interval:
#> ... by normal method:
#> -0.124 , 0.2575
#> =======================================================
#> INDIRECT EFFECT OF X ...
#> Indirect effect parametric estimate:
#> -0.3018684
#> Indirect effect bootstrap estimate:
#> -0.2930111
#> Indirect effect bootstrap confidence interval:
#> ... by normal method:
#> -0.3663 , -0.2197
#> =======================================================
#> Computation time:
#> Time difference of 27.21542 secs
#> =======================================================
#> Time-varying Effects Model Predicting MEDIATOR from TREATMENT:
#> Response variable: MEDIATOR
#> Response outcome distribution type: gaussian
#> Time interval: 0.01 to 1
#> Number of subjects: 500
#> Effects specified as time-varying: (Intercept), TREATMENT, S
#> You can use the plot function to view the plots.
#>
#> Back-end model fitted in mgcv::bam function:
#> Method fREML
#> Formula:
#> MEDIATOR ~ TREATMENT + S + s(time, bs = "ps", by = NA, pc = 0,
#> k = 7, fx = FALSE) + s(time, bs = "ps", by = TREATMENT, pc = 0,
#> m = c(2, 1), k = 7, fx = FALSE) + s(time, bs = "ps", by = S,
#> pc = 0, m = c(2, 1), k = 7, fx = FALSE)
#> Pseudolikelihood AIC: 83986.37
#> Pseudolikelihood BIC: 84047.99
#>
#> =======================================================
#> Functional Regression Model Predicting OUTCOME from MEDIATOR
#> and TREATMENT:
#>
#> Family: gaussian
#> Link function: identity
#>
#> Formula:
#> OUTCOME ~ s(x = MEDIATOR.tmat, by = L.MEDIATOR) + TREATMENT +
#> S
#>
#> Estimated degrees of freedom:
#> 2 total = 5
#>
#> REML score: 729.0772
#> Scalar terms in functional regression model:
#> (Intercept) TREATMENT S
#> -0.025787118 0.112576286 -0.001643264
#> =======================================================
#> Parametric model for Predicting OUTCOME from TREATMENT
#> ignoring MEDIATOR:
#>
#> Call:
#> glm(formula = glm_formula)
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 0.36962 0.08331 4.437 1.12e-05 ***
#> TREATMENT -0.16672 0.09601 -1.737 0.0831 .
#> S -0.02827 0.09593 -0.295 0.7684
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> (Dispersion parameter for gaussian family taken to be 1.149205)
#>
#> Null deviance: 574.74 on 499 degrees of freedom
#> Residual deviance: 571.16 on 497 degrees of freedom
#> AIC: 1493.5
#>
#> Number of Fisher Scoring iterations: 2
#>
#> =======================================================
plot(model_with_tve_covariate, what_plot="tvem")
Notice that we now get a time-varying effects plot which includes a coefficient for the effect of \(S\), in addition to the coefficients for the effect of \(X\). In this example \(\alpha_S(t)\) is not significantly nonzero at any time and does not change over time, so there is no evidence that covariate \(S\) is important in predicting the outcome \(Y\). As before, the \(X=1\) group has lower values on the mediator than the \(X=0\) group (i.e., their contrast is negative) at least at later time points.
Alternatively, by specifying tie_covariates_on_mediator
instead of tve_covariates_on_mediator, the covariate can be
assumed to have a time-invariant effect rather than a time-varying
effect on \(M\).
model_with_tie_covariate <- funmediation(data=data_with_covariate,
treatment=X,
mediator=M,
outcome=Y,
tie_covariates_on_mediator = ~S,
covariates_on_outcome = ~S,
id=subject_id,
time=t,
nboot=10)
#> Warning in funmediation(data = data_with_covariate, treatment = X, mediator =
#> M, : The number of bootstrap samples was very small and results will be
#> dubious.
#> Warning in norm.inter(t, (1 + c(conf, -conf))/2): extreme order statistics used
#> as endpoints
#> Warning in norm.inter(t, alpha): extreme order statistics used as endpoints
#> Warning in norm.inter(t, (1 + c(conf, -conf))/2): extreme order statistics used
#> as endpoints
#> Warning in norm.inter(t, alpha): extreme order statistics used as endpoints
print(model_with_tie_covariate)
#> =======================================================
#> Functional Regression Mediation Function Output
#> =======================================================
#> TREATMENT: X
#> MEDIATOR: M (Assumed Continuous)
#> OUTCOME: Y (Assumed Continuous)
#> =======================================================
#> DIRECT EFFECT OF X ...
#> Direct effect estimate:
#> 0.1125763
#> Direct effect standard error:
#> 0.1009101
#> Direct effect bootstrap estimate:
#> 0.11074
#> Direct effect bootstrap confidence interval:
#> ... by normal method:
#> -0.0443 , 0.2658
#> =======================================================
#> INDIRECT EFFECT OF X ...
#> Indirect effect parametric estimate:
#> -0.302025
#> Indirect effect bootstrap estimate:
#> -0.2680548
#> Indirect effect bootstrap confidence interval:
#> ... by normal method:
#> -0.3417 , -0.1944
#> =======================================================
#> Computation time:
#> Time difference of 21.21128 secs
#> =======================================================
#> Time-varying Effects Model Predicting MEDIATOR from TREATMENT:
#> Response variable: MEDIATOR
#> Response outcome distribution type: gaussian
#> Time interval: 0.01 to 1
#> Number of subjects: 500
#> Effects specified as time-varying: (Intercept), TREATMENT
#> You can use the plot function to view the plots.
#>
#> Effects specified as non-time-varying:
#> estimate standard_error
#> S -0.03008131 0.05790221
#>
#> Back-end model fitted in mgcv::bam function:
#> Method fREML
#> Formula:
#> MEDIATOR ~ TREATMENT + S + s(time, bs = "ps", by = NA, pc = 0,
#> k = 7, fx = FALSE) + s(time, bs = "ps", by = TREATMENT, pc = 0,
#> m = c(2, 1), k = 7, fx = FALSE)
#> Pseudolikelihood AIC: 83993.2
#> Pseudolikelihood BIC: 84042.46
#>
#> =======================================================
#> Functional Regression Model Predicting OUTCOME from MEDIATOR
#> and TREATMENT:
#>
#> Family: gaussian
#> Link function: identity
#>
#> Formula:
#> OUTCOME ~ s(x = MEDIATOR.tmat, by = L.MEDIATOR) + TREATMENT +
#> S
#>
#> Estimated degrees of freedom:
#> 2 total = 5
#>
#> REML score: 729.0772
#> Scalar terms in functional regression model:
#> (Intercept) TREATMENT S
#> -0.025787118 0.112576286 -0.001643264
#> =======================================================
#> Parametric model for Predicting OUTCOME from TREATMENT
#> ignoring MEDIATOR:
#>
#> Call:
#> glm(formula = glm_formula)
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 0.36962 0.08331 4.437 1.12e-05 ***
#> TREATMENT -0.16672 0.09601 -1.737 0.0831 .
#> S -0.02827 0.09593 -0.295 0.7684
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> (Dispersion parameter for gaussian family taken to be 1.149205)
#>
#> Null deviance: 574.74 on 499 degrees of freedom
#> Residual deviance: 571.16 on 497 degrees of freedom
#> AIC: 1493.5
#>
#> Number of Fisher Scoring iterations: 2
#>
#> =======================================================
plot(model_with_tie_covariate, what_plot="pfr")
The functional effect plot here is very similar to what it was before. Now there are scalar estimates for \(\alpha_S=-0.03008\) and \(\beta_S=-0.00164\) but no plot for \(\alpha_S(t)\) because \(\alpha_S\) is now assumed not to vary over time. Again, \(S\) is not statistically significant either as a predictor of the mediator or of the outcome, because the coefficients are not large relative to their standard errors.
Using a Binary Mediator
The mediator may be binary rather than continuous, in which case the the TVEM uses the logistic link function. For purposes of this vignette demonstration, we dichotomize \(M\) in the simulated dataset and use that as our binary mediator (this usually would not be advisable for real data because it involves loss of information).
data_with_covariate$binary_M <- 1*(data_with_covariate$M > mean(data_with_covariate$M))
head(data_with_covariate)
#> subject_id t X M Y S binary_M
#> 1 1 0.01 1 -2.9541303 -0.44206 0 0
#> 2 1 0.05 1 -0.5834873 -0.44206 0 0
#> 3 1 0.06 1 0.6919501 -0.44206 0 1
#> 4 1 0.08 1 0.4411048 -0.44206 0 1
#> 5 1 0.15 1 -1.2007641 -0.44206 0 0
#> 6 1 0.17 1 -1.0369571 -0.44206 0 0Now we call the funmediation function and include the
binary_mediator=TRUE argument.
model_with_binary_M <- funmediation(data=data_with_covariate,
treatment=X,
mediator=binary_M,
outcome=Y,
id=subject_id,
time=t,
binary_mediator=TRUE,
nboot=10)
#> Warning in funmediation(data = data_with_covariate, treatment = X, mediator =
#> binary_M, : The number of bootstrap samples was very small and results will be
#> dubious.
#> Warning in norm.inter(t, (1 + c(conf, -conf))/2): extreme order statistics used
#> as endpoints
#> Warning in norm.inter(t, alpha): extreme order statistics used as endpoints
#> Warning in norm.inter(t, (1 + c(conf, -conf))/2): extreme order statistics used
#> as endpoints
#> Warning in norm.inter(t, alpha): extreme order statistics used as endpointsprint(model_with_binary_M)
#> =======================================================
#> Functional Regression Mediation Function Output
#> =======================================================
#> TREATMENT: X
#> MEDIATOR: binary_M (Assumed Binary)
#> OUTCOME: Y (Assumed Continuous)
#> =======================================================
#> DIRECT EFFECT OF X ...
#> Direct effect estimate:
#> 0.03910025
#> Direct effect standard error:
#> 0.1004427
#> Direct effect bootstrap estimate:
#> 0.07949018
#> Direct effect bootstrap confidence interval:
#> ... by normal method:
#> -0.1906 , 0.3497
#> =======================================================
#> INDIRECT EFFECT OF X ...
#> Indirect effect parametric estimate:
#> -0.8550336
#> Indirect effect bootstrap estimate:
#> -0.2207669
#> Indirect effect bootstrap confidence interval:
#> ... by normal method:
#> -0.3184 , -0.1231
#> =======================================================
#> Computation time:
#> Time difference of 30.91431 secs
#> =======================================================
#> Time-varying Effects Model Predicting MEDIATOR from TREATMENT:
#> Response variable: MEDIATOR
#> Response outcome distribution type: binomial
#> Time interval: 0.01 to 1
#> Number of subjects: 500
#> Effects specified as time-varying: (Intercept), TREATMENT
#> You can use the plot function to view the plots.
#>
#> Back-end model fitted in mgcv::bam function:
#> Method fREML
#> Formula:
#> MEDIATOR ~ TREATMENT + s(time, bs = "ps", by = NA, pc = 0, k = 7,
#> fx = FALSE) + s(time, bs = "ps", by = TREATMENT, pc = 0,
#> m = c(2, 1), k = 7, fx = FALSE)
#> Pseudolikelihood AIC: 27100.7
#> Pseudolikelihood BIC: 27134.79
#>
#> =======================================================
#> Functional Regression Model Predicting OUTCOME from MEDIATOR
#> and TREATMENT:
#>
#> Family: gaussian
#> Link function: identity
#>
#> Formula:
#> OUTCOME ~ s(x = MEDIATOR.tmat, by = L.MEDIATOR) + TREATMENT
#>
#> Estimated degrees of freedom:
#> 2 total = 4
#>
#> REML score: 731.6047
#> Scalar terms in functional regression model:
#> (Intercept) TREATMENT
#> -0.62912300 0.03910025
#> =======================================================
#> Parametric model for Predicting OUTCOME from TREATMENT
#> ignoring MEDIATOR:
#>
#> Call:
#> glm(formula = glm_formula)
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 0.35608 0.06942 5.129 4.18e-07 ***
#> TREATMENT -0.16720 0.09591 -1.743 0.0819 .
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> (Dispersion parameter for gaussian family taken to be 1.147098)
#>
#> Null deviance: 574.74 on 499 degrees of freedom
#> Residual deviance: 571.25 on 498 degrees of freedom
#> AIC: 1491.6
#>
#> Number of Fisher Scoring iterations: 2
#>
#> =======================================================
plot(model_with_binary_M, what_plot="coef")
The plots are similar to those shown in the previous example with numerical mediator, except that the TVEM coefficient functions for the intercept and treatment effects on the outcome are interpreted on the logit scale as in logistic regression.
Note that if the binary_mediator=TRUE argument is not
specified, the TVEM will be fit using an identity link (hence a linear
model at any given time \(t\)) as for a
normally distributed outcome. The function would not detect that
binary_M consists of zeroes and ones and no other numbers.
The identity link is probably not the best approach but might be
interesting as a comparison. If the mediator is binary, it is probably
better to use binary_mediator=TRUE in order to specify that
the \(\alpha\) coefficients are to be
interpreted on the logit (logistic regression) scale. A future version
of this package might allow log link functions as an alternative to
logit link functions, but this is currently not available.
