regmed: Mediation model for 1 exposure, multiple mediators, and 1 outcome variable
This section focuses on the original version of regmed that handles a single exposure and a single outcome. The user-level functions for single outcome/exposure are:
regmed.prefilter()pre-filters mediators if number of mediators is greater than the number of subjectsregmed.grid()penalized model for grid of lambdasregmed.grid.bestfit()get best fit model from grid of fitsregmed.fit()fit penalized model with specified lambdaplot.regmed.grid()plot results of regmed.gridregmed.edges()create edges from a fit of a model for use in plotting and model fit by lavaanplot.regmed.edges()plot directed graph based on edges: exposure -> mediator -> outcomeregmed.lavaan.model()setup model to estimate unpenalized parameters by lavaann sem functionregmed.lavaan.dat()setup data to input to lavaann sem functionsummary.lavaan()summarize lavaan sem model fit
Pre-filter mediators
If the number of mediators exceeds the sample size, model fitting can become unstable. To reduce the number of mediators to fit, the regmed.prefilter uses sure independence screening (Fan & Lv, 2008) to reduce the number of potential mediators. This is based on ranking marginal correlations and then selecting the highest ranked values such that the number of parameters is less than the sample size. Because mediation depends on the two correlations, \(cor(x,med)\) and \(cor(med,y)\) we rank the absolute values of their products, \(|cor(x, med) * cor(med, y)|\), and choose the highest k ranked values to determine which potential mediators to include in penalized mediation models. If k is not specified, the default value of k is n/2, where n is the sample size, because each mediator results in two parameters alpha and beta.
To speed calculations solely for demonstration purposes, we choose k=10 mediators. Also, we select the first exposure variable and the first outcome variable to demonstrate methods for regmed. The function regmed.prefilter returns a list of x, mediator, and y, with each of these variables centered and scaled by default, and accounts for missing data by subsetting to subjects that are not missing any of the variables.
[1] "med.1" "med.2" "med.74" "med.88" "med.90" "med.99" "med.112"
[8] "med.129" "med.133" "med.178"Grid of lambda penalties
To fit a series of models over a grid of penalty “lambda” values, it is necessary to define a vector of lambda values. It is best to arrange this vector from largest to smallest lambda values to assure that the largest lambda (first value in the vector) results in no selected parameters (alpha, beta, delta), and the smallest lambda (last value in the vector) selects multiple parameters, and that a lambda between the largest and smallest results in a minimum Bayesian Information Criterion (BIC). Starting with a large lambda and then decreasing lambda values is best, because this provides a “warm start” for each subsequent fit, with final fitted parameter values for a specific lambda used as initial values for the next lambda value. The demonstration below extracts variables that have been processed by regmed.prefilter.
Fit regmed for grid of lambda penalties
lambda.grid <- seq(from = 0.4, to = 0.01, by = -0.05)
x1 <- dat.filter$x
y1 <- dat.filter$y
med <- dat.filter$mediator
fit.grid <- regmed.grid(x1, med, y1, lambda.grid, frac.lasso = 0.8)The plot method for regmed.grid gives two figures that represent the size of the coefficients over the grid of lambdas, and the BIC by by the grid of lambdas.
Call:
regmed.grid(x = x1, mediator = med, y = y1, lambda.vec = lambda.grid,
frac.lasso = 0.8)
lambda converge iter df bic
1 0.40 TRUE 289 9 1048.155
2 0.35 TRUE 8 9 1041.148
3 0.30 TRUE 5 12 1047.839
4 0.25 TRUE 14 15 1053.422
5 0.20 TRUE 6 17 1053.775
6 0.15 TRUE 14 19 1055.445
7 0.10 TRUE 15 21 1059.531
8 0.05 TRUE 9 21 1059.378
Select best model from grid by minimum BIC
Call:
NULL
Coefficients:
alpha beta alpha*beta
med.1 0.37711665 0.1143248 0.04311379
med.2 0.32632288 0.0000000 0.00000000
med.74 0.00000000 0.0000000 0.00000000
med.88 -0.10086322 0.0000000 0.00000000
med.90 0.00000000 0.0000000 0.00000000
med.99 0.00000000 0.0000000 0.00000000
med.112 -0.14386830 0.0000000 0.00000000
med.129 0.04906179 0.0000000 0.00000000
med.133 0.05666569 0.0000000 0.00000000
med.178 0.00000000 0.0000000 0.00000000
sum of alpha*beta = 0.04311379
delta = 0.5050343
sum of delta + alpha*beta = 0.5481481
var(x) = 1
var(y) = 1.101601 Fit single regmed model with user-specified lambda
If user prefers to specify a lambda penalty and not search through a grid, the function regmed.fit is available.
Call:
regmed.fit(x = x1, mediator = med, y = y1, lambda = 0.3, frac.lasso = 0.8)
Coefficients:
alpha beta alpha*beta
med.1 0.40109014 0.126420020 5.070582e-02
med.2 0.35582882 0.000000000 0.000000e+00
med.74 0.00000000 0.000000000 0.000000e+00
med.88 -0.13906482 0.000000000 0.000000e+00
med.90 0.00000000 0.000000000 0.000000e+00
med.99 -0.01045770 -0.001413514 1.478211e-05
med.112 -0.17587944 0.000000000 0.000000e+00
med.129 0.07206755 0.000000000 0.000000e+00
med.133 0.07530145 0.000000000 0.000000e+00
med.178 -0.01914775 0.000000000 0.000000e+00
sum of alpha*beta = 0.05072061
delta = 0.5386689
sum of delta + alpha*beta = 0.5893895
var(x) = 1
var(y) = 1.116486 Plot directed graph of fit model
The demonstration below shows how to determine edges, where an edge is defined by vertex-1 pointing to vertex-2, and the vertices are the variables selected in a model. For example, if alpha[1] and beta[1] are both estimated to be non-zero, then x -> med[1] and med[1] -> y. The function regmed.edges has an option to choose how vertices and edges are selected. By type=“mediators”, vertices and edges are selected only if the product \(alpha * beta\) is non-zero, because this is required for the definition of a mediator. By specifying type = “any”, all vertices and edges that represent non-zero parameters are selected. For example, if alpha[1] is non-zero, and beta[1] is zero, the edge x -> med[1] is selected, but the edge med[1] -> y is not selected. To determine if terms are zero, a threshold variable eps is used (see help for regmed.edges).
For fit.best, the summary above shows that only med.1 is selected as a mediator, because \(|alpha * beta| > eps\). In this case, choosing type=“mediators” selects only med.1 to be included in edges.
In contrast, choosing type=“any” for fit.best includes more mediator variables with edges, even though they do not all meet the definition of a mediator.
Lavaan sem: Estimate parameters of selected model without imposing penalty
Because penalized models can overly shrink parameter estimates towards zero, it can be beneficial to refit the selected model without penalties. The model without penalties is a standard structural equation model, which can be fit with the function sem from the package lavaan. To facilitate this step, the function regmed.lavaan.model uses the edges and fit of a model to create a text string that represents the model syntax for sem. This function abides by our assumption that the residual covariances between x and med, between x and y, and between med and y, are all zero. It also assumes that the covariances of the mediators, a matrix in the fit of a regmed object, are fixed in the sem model fitting.
The demonstration below shows how to set up the data set and models for sem. The function regmed.lavaan.dat assures that x, med, y, are centered and scaled and subset to subjects without any missing data.
dat <- regmed.lavaan.dat(x1, med, y1)
mod.best <- regmed.lavaan.model(edges.med, fit.best)
mod.any <- regmed.lavaan.model(edges.any, fit.best)The demonstration below shows how to use lavaan:::sem() to fit specified models and view results
lhs op rhs est se z pvalue ci.lower ci.upper
1 med.1 ~ x1 0.551 0.084 6.577 0.000 0.387 0.716
2 y1 ~ med.1 0.293 0.093 3.147 0.002 0.110 0.475
3 y1 ~ x1 0.415 0.093 4.445 0.000 0.232 0.598
4 y1 ~~ y1 0.602 0.085 7.071 0.000 0.435 0.769 lhs op rhs est se z pvalue ci.lower ci.upper
1 med.1 ~ x1 0.551 0.084 6.577 0.000 0.387 0.716
2 med.2 ~ x1 0.489 0.088 5.576 0.000 0.317 0.661
3 med.88 ~ x1 -0.197 0.099 -2.002 0.045 -0.390 -0.004
4 med.112 ~ x1 -0.264 0.097 -2.723 0.006 -0.454 -0.074
5 med.129 ~ x1 0.232 0.098 2.369 0.018 0.040 0.423
6 med.133 ~ x1 0.200 0.098 2.032 0.042 0.007 0.393
7 y1 ~ med.1 0.293 0.093 3.147 0.002 0.110 0.475
8 y1 ~ x1 0.415 0.093 4.445 0.000 0.232 0.598
9 y1 ~~ y1 0.602 0.085 7.071 0.000 0.435 0.769
