SEMgraph (version 1.2.2)

SEMrun: Fit a graph as a Structural Equation Model (SEM)

Description

SEMrun() converts a (directed, undirected, or mixed) graph to a SEM and fits it. If a binary group variable (i.e., case/control) is present, node-level or edge-level perturbation is evaluated. This function can handle loop-containing models, although multiple links between the same two nodes (including self-loops and mutual interactions) and bows (i.e., a directed and a bidirected link between two nodes) are not allowed.

Usage

SEMrun(
  graph,
  data,
  group = NULL,
  fit = 0,
  algo = "lavaan",
  start = NULL,
  SE = "standard",
  n_rep = 1000,
  limit = 100,
  ...
)

Value

A list of 5 objects:

1. "fit", SEM fitted lavaan, ricf, or cggm object, depending on the MLE method specified by the algo argument;

2. "gest" or "dest", a data.frame of node-specific ("gest") or edge-specific ("dest") group effect estimates and P-values;

3. "model", SEM model as a string if algo = "lavaan", and NULL otherwise;

4. "graph", the induced subgraph of the input network mapped on data variables. Graph edges (i.e., direct effects) with P-value < 0.05 will be highlighted in red (beta > 0) or blue (beta < 0). If a group vector is given, nodes with significant group effect (P-value < 0.05) will be red-shaded (beta > 0) or lightblue-shaded (beta < 0);

5. "data", input data subset mapping graph nodes, plus group at the first column (if no group is specified, this column will take NA values).

Arguments

graph

An igraph object.

data

A matrix whith rows corresponding to subjects, and columns to graph nodes (variables).

group

A binary vector. This vector must be as long as the number of subjects. Each vector element must be 1 for cases and 0 for control subjects. If NULL (default), group influence will not be considered.

fit

A numeric value indicating the SEM fitting mode. If fit = 0 (default), no group effect is considered. If fit = 1, a "common" model is used to evaluate group effects on graph nodes. If fit = 2, a two-group model is used to evaluate group effects on graph edges.

algo

MLE method used for SEM fitting. If algo = "lavaan" (default), the SEM will be fitted using the NLMINB solver from lavaan R package, with standard errors derived from the expected Fisher information matrix. If algo = "ricf", the model is fitted via residual iterative conditional fitting (RICF; Drton et al. 2009), with standard error derived from randomization or bootstrap procedures. If algo = "cggm", model fitting is based on constrained Gaussian Graphical Modeling (CGGM), with DAG nodewise Lasso procedure and de-biasing asymptotic inference (Jankova & Van De Geer, 2019).

start

Starting value of SEM parameters for algo = "lavaan". If start is NULL (default), the algorithm will determine the starting values. If start is a numeric value, it will be used as a scaling factor for the edge weights in the graph object (graph attribute E(graph)$weight). For instance, a scaling factor is useful when weights have fixed values (e.g., 1 for activated, -1 for repressed, and 0 for unchanged interaction). Fixed values may compromise model fitting, and scaling them is a safe option to avoid this problem. As a rule of thumb, to our experience, start = 0.1 generally performs well with (-1, 0, 1) weights.

SE

If "standard" (default), with algo = "lavaan", conventional standard errors are computed based on inverting the observed information matrix. If "none", no standard errors are computed.

n_rep

Number of randomization replicates (default = 1000), for permutation flip or boostrap samples, if algo = "ricf".

limit

An integer value corresponding to the network size (i.e., number of nodes). Beyond this limit, the execution under algo = "lavaan" will run with SE = "none", if fit = 0, or will be ridirected to algo = "ricf", if fit = 1, or to algo = "cggm", if fit = 2. This redirection is necessary to reduce the computational demand of standard error estimation by lavaan. Increasing this number will enforce lavaan execution when algo = "lavaan".

...

Currently ignored.

Author

Mario Grassi

Details

SEMrun maps data onto the input graph and converts it into a SEM. Directed connections (X -> Y) are interpreted as direct causal effects, while undirected, mutual, and bidirected connections are converted into model covariances. SEMrun output contains different sets of parameter estimates. Beta coefficients (i.e., direct effects) are estimated from directed interactions and residual covariances (psi coefficients) from bidirected, undirected, or mutual interactions. If a group variable is given, exogenous group effects on nodes (gamma coefficients) or edges (delta coefficients) will be estimated. By default, maximum likelihood parameter estimates and P-values for parameter sets are computed by conventional z-test (= estimate/SE), and fits it through the lavaan function, via Maximum Likelihood Estimation (estimator = "ML", default estimator in lavOptions). In case of high dimensionality (n.variables >> n.subjects), the covariance matrix could not be semi-definite positive and thus parameter estimates could not be done. If this happens, covariance matrix regularization is enabled using the James-Stein-type shrinkage estimator implemented in the function pcor.shrink of corpcor R package. Argument fit determines how group influence is evaluated in the model, as absent (fit = 0), node perturbation (fit = 1), or edge perturbation (fit = 2). When fit = 1, the group is modeled as an exogenous variable, influencing all the other graph nodes. When fit = 2, SEMrun estimates the differences of the beta and/or psi coefficients (network edges) between groups. This is equivalent to fit a separate model for cases and controls, as opposed to one common model perturbed by the exogenous group effect. Once fitted, the two models are then compared to assess significant edge (i.e., direct effect) differences (d = beta1 - beta0). P-values for parameter sets are computed by z-test (= d/SE), through lavaan. As an alternative to standard P-value calculation, SEMrun may use either RICF (randomization or bootstrap P-values) or GGM (de-biased asymptotically normal P-values) methods. These algorithms are much faster than lavaan in case of large input graphs.

References

Pearl J (1998). Graphs, Causality, and Structural Equation Models. Sociological Methods & Research., 27(2):226-284. <https://doi.org/10.1177/0049124198027002004>

Yves Rosseel (2012). lavaan: An R Package for Structural Equation Modeling. Journal of Statistical Software, 48(2): 1-36. <https://www.jstatsoft.org/v48/i02/>

Pepe D, Grassi M (2014). Investigating perturbed pathway modules from gene expression data via Structural Equation Models. BMC Bioinformatics, 15: 132. <https://doi.org/10.1186/1471-2105-15-132>

Drton M, Eichler M, Richardson TS (2009). Computing Maximum Likelihood Estimated in Recursive Linear Models with Correlated Errors. Journal of Machine Learning Research, 10(Oct): 2329-2348. <https://www.jmlr.org/papers/volume10/drton09a/drton09a.pdf>

Jankova, J., & Van De Geer, S (2019). Inference in high-dimensional graphical models. In Handbook of Graphical Models (2019). Chapter 14 (sec. 14.2): 325-349. Chapman & Hall/CRC. ISBN: 9780429463976

Hastie T, Tibshirani R, Friedman J. (2009). The Elements of Statistical Learning (2nd ed.). Springer Verlag. ISBN: 978-0-387-84858-7

Grassi M, Palluzzi F, Tarantino B (2022). SEMgraph: An R Package for Causal Network Analysis of High-Throughput Data with Structural Equation Models. Bioinformatics, 38 (20), 4829–4830 <https://doi.org/10.1093/bioinformatics/btac567>

See Also

See fitAncestralGraph and fitConGraph for RICF algorithm and constrained GGM algorithm details, respectively.

Examples

#### Model fitting (no group effect)

sem0 <- SEMrun(graph = sachs$graph, data = log(sachs$pkc))
summary(sem0$fit)
head(parameterEstimates(sem0$fit))

# Graphs
gplot(sem0$graph, main = "significant edge weights")
plot(sem0$graph, layout = layout.circle, main = "significant edge weights")


#### Model fitting (common model, group effect on nodes)

sem1 <- SEMrun(graph = sachs$graph, data = log(sachs$pkc),
               group = sachs$group)

# Fitting summaries
summary(sem1$fit)
print(sem1$gest)
head(parameterEstimates(sem1$fit))

# Graphs
gplot(sem1$graph, main = "Between group node differences")
plot(sem1$graph, layout = layout.circle, main = "Between group node differences")


#### Two-group model fitting (group effect on edges)

sem2 <- SEMrun(graph = sachs$graph, data = log(sachs$pkc),
               group = sachs$group,
               fit = 2)

# Summaries
summary(sem2$fit)
print(sem2$dest)
head(parameterEstimates(sem2$fit))

# Graphs
gplot(sem2$graph, main = "Between group edge differences")
plot(sem2$graph, layout = layout.circle, main = "Between group edge differences")

# \donttest{
# Fitting and visualization of a large pathway:

g <- kegg.pathways[["Neurotrophin signaling pathway"]]
G <- properties(g)[[1]]
summary(G)

# Nonparanormal(npn) transformation
als.npn <- transformData(alsData$exprs)$data

g1 <- SEMrun(G, als.npn, alsData$group, algo = "cggm")$graph
g2 <- SEMrun(g1, als.npn, alsData$group, fit = 2, algo = "cggm")$graph

# extract the subgraph with node and edge differences
g2 <- g2 - E(g2)[-which(E(g2)$color != "gray50")]
g <- properties(g2)[[1]]

# plot graph
E(g)$color<- E(g2)$color[E(g2) %in% E(g)]
gplot(g, l="fdp", psize=40, main="node and edge group differences")
# }