--- title: "Stepwise Gaussian graphical model selection" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Stepwise Gaussian graphical model selection} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 6) library(psychnets) has_cograph <- requireNamespace("cograph", quietly = TRUE) ``` ## What stepwise model selection does A stepwise Gaussian graphical model estimates the conditional dependence structure among a set of continuous variables by combining model selection with maximum-likelihood estimation. As in all Gaussian graphical models, each edge represents a partial correlation between two variables after conditioning on all remaining variables in the network. An edge therefore indicates a conditional association that cannot be explained by the other measured variables, whereas the absence of an edge implies that the observed association is adequately accounted for by the remaining variables under the fitted model. The network structure is determined through an iterative stepwise search in which candidate graphs are compared using an information criterion, typically the Extended Bayesian Information Criterion (EBIC). Edges are sequentially added or removed until no further improvement in the selection criterion is obtained. Unlike penalized estimators such as the graphical lasso, sparsity is achieved through explicit model selection rather than shrinkage of the precision matrix. Following model selection, the precision matrix is re-estimated by maximum likelihood on the selected graph without regularization. Consequently, the retained edge weights are not biased toward zero by a penalty term and may be preferred when the objective is parameter estimation and interpretation of partial correlations within the selected network. This approach is particularly appropriate when sample size is sufficiently large relative to the number of variables, allowing reliable estimation without regularization. When the sample size is moderate to large and the underlying network is not excessively dense, stepwise Gaussian graphical models often recover graph structures similar to those obtained using the EBIC graphical lasso, particularly when both methods employ the same information criterion for model selection. However, because the graphical lasso performs simultaneous estimation and regularization, whereas the stepwise approach separates model selection from parameter estimation, the two methods may differ in finite samples. The graphical lasso generally exhibits greater numerical stability in high-dimensional settings ((p \geq n)) or when estimation is ill-conditioned, while the stepwise estimator avoids the shrinkage bias inherent to penalized estimation by refitting the selected model using unpenalized maximum likelihood. ## The data The worked example uses `SRL_GPT`, which contains 300 observations on five self-regulated learning constructs. Cognitive strategy use (`CSU`), intrinsic value (`IV`), self-efficacy (`SE`), self-regulation (`SR`), and test anxiety (`TA`) are composite scores recorded on a 1 to 7 scale. Rows are observations and columns are network nodes. ```{r data-preview} head(SRL_GPT) ``` The analysis treats these composite scores as continuous. Before applying the same workflow to another data set, the analyst should inspect missingness, distributions, unusual observations, and the measurement scale. The present data are complete, so every pairwise correlation uses all 300 observations. ## Fitting the network with `psychnet()` `psychnet()` estimates the stepwise Gaussian graphical model when `method = "ggm"`. It returns a fitted `psychnet` object containing the selected network and the quantities needed for subsequent analysis. The default $\gamma = 0$ uses the ordinary Bayesian information criterion, and the default stepwise search considers single-edge additions and removals. ```{r fit-network} stepwise_net <- psychnet(data = SRL_GPT, method = "ggm") stepwise_net ``` The fitted network contains 5 nodes and 9 edges. One of the 10 possible pairs is removed. The printed KKT residual is $1.33 \times 10^{-15}$, which will be interpreted in the numerical check below. ## Inspecting the selected edges with `summary()` `summary()` prints the fitted model and returns its edge table invisibly. The table contains `from`, `to`, and `weight`. Each weight is an unshrunk partial correlation estimated under the selected network structure. ```{r summarize-network} summary(stepwise_net) ``` The retained weights range from -0.370 to 0.422. The largest positive edge joins `CSU` and `IV` ($r = 0.422$). Cognitive strategy use and intrinsic value remain positively related above and beyond self-efficacy, self-regulation, and test anxiety. The `SR`-`TA` edge is negative ($r = -0.370$). Higher self-regulation is associated with lower test anxiety after the other learning constructs have been accounted for. The model removes the `CSU`-`TA` edge. This result means that the selected structure does not retain a unique conditional association between cognitive strategy use and test anxiety. It does not imply that their ordinary correlation must be zero. The `IV`-`TA` and `SE`-`TA` edges are positive but modest, at 0.157 and 0.137. Their magnitudes warrant cautious interpretation. Every edge is an adjusted association, so the network alone provides no direction of causation. ## Checking the fitted solution with `certificate()` `certificate()` reports whether the fitted precision matrix satisfies the conditions that define the maximum-likelihood estimate for the selected graph. It returns the columns `method`, `certificate`, `kind`, and `certified`. For this model, `certificate` is the largest numerical discrepancy across the defining conditions, and `certified` compares that value with a default tolerance of $10^{-6}$. ```{r fit-certificate} certificate(stepwise_net) ``` The residual is $1.33 \times 10^{-15}$ and `certified` is `TRUE`. This value is at the order of machine precision. It indicates that the final precision matrix satisfies the selected model's numerical conditions to numerical precision. The certificate checks computation. It does not evaluate sampling uncertainty, measurement quality, model assumptions, or causal interpretation. ## Comparing the selected model with the graphical lasso `psychnet()` estimates the graphical lasso when `method = "glasso"`. Fitting it to the same data provides a useful comparison because both methods estimate partial-correlation networks, but they select and estimate edges differently. ```{r fit-glasso} glasso_net <- psychnet(data = SRL_GPT, method = "glasso") glasso_net ``` ```{r summarize-glasso} summary(glasso_net) ``` The graphical lasso retains all 10 possible edges, including a weak `CSU`-`TA` edge of 0.049. The stepwise model removes that edge and retains 9. For shared edges, the stepwise estimates are often slightly larger in absolute magnitude because its final refit applies no shrinkage. The `CSU`-`IV` edge is 0.422 in the stepwise model and 0.412 in the graphical lasso. The `SR`-`TA` edge is -0.370 in the stepwise model and -0.350 in the graphical lasso. This comparison does not identify a universally preferable model. The stepwise model reports unshrunk partial correlations for a selected structure. The graphical lasso estimates all retained weights under an $\ell_1$ penalty. The research aim and the planned approach to validation should guide the choice. ## Describing node position with `net_centralities()` `net_centralities()` computes node-level summaries from the fitted graph. The default output contains `node`, `strength`, and `expected_influence`. Strength is the sum of absolute edge weights. Expected influence is the signed sum, so negative associations offset positive associations. ```{r node-centrality} net_centralities(stepwise_net) ``` Self-regulation has the largest strength (1.259), followed by cognitive strategy use (1.152). Self-regulation has expected influence 0.519 because its negative edge with test anxiety offsets its positive edges with the learning constructs. Cognitive strategy use has strength and expected influence of 1.152 because all of its retained edges are positive. Test anxiety has the smallest strength (0.663) and a slightly negative expected influence (-0.077). These values describe its position in the estimated network. They do not show that self-regulation is an intervention target or that changing one node would cause changes in its neighbours. ## Quantifying node predictability with `net_predict()` `net_predict()` computes how much of each node's variance is accounted for by the remaining nodes. For a Gaussian graphical model, it returns a tidy table with `node`, `type`, `metric`, `predictability`, and `accuracy`. Continuous nodes use the $R^2$ metric, while `accuracy` is `NA` because classification accuracy does not apply. ```{r node-predictability} net_predict(stepwise_net) ``` Cognitive strategy use has the highest predictability ($R^2 = 0.833$). Intrinsic value ($R^2 = 0.783$), self-regulation ($R^2 = 0.792$), and self-efficacy ($R^2 = 0.728$) are also strongly accounted for by the remaining constructs. Test anxiety has the lowest predictability ($R^2 = 0.166$), so most of its variance is not accounted for by this five-node model. These values describe in-sample explained variance. Out-of-sample prediction requires a separate validation sample or a resampling design that evaluates performance on observations excluded from model fitting. ## Visualizing the network with `cograph::splot()` `cograph::splot()` draws the fitted network when the optional `cograph` package is available. Psychological styling uses green for positive edges and red for negative edges. The predictability ring is requested explicitly for this model. ```{r plot-network, eval = has_cograph} cograph::splot(stepwise_net, psych_styling = TRUE, predictability = TRUE) ``` The figure summarizes the signs and relative magnitudes of the retained edges. Node placement comes from a layout algorithm and has no measurement scale. Interpretation should therefore use the printed edge and node tables as the primary numerical results. ## Sensitivity to rank correlations `psychnet()` accepts `cor_method = "spearman"` for a rank-based sensitivity analysis. This analysis asks whether the main structure persists when estimation depends on monotonic ranks and is less sensitive to extreme values. ```{r fit-spearman} spearman_net <- psychnet(data = SRL_GPT, method = "ggm", cor_method = "spearman") spearman_net ``` ```{r summarize-spearman} summary(spearman_net) ``` The Spearman model contains 7 edges, compared with 9 in the Pearson model. The positive `CSU`-`IV`, `CSU`-`SE`, `CSU`-`SR`, `IV`-`SR`, and `SE`-`SR` edges remain, as does the negative `SR`-`TA` edge. The `IV`-`SE` and `SE`-`TA` edges are removed. The main positive learning cluster and negative self-regulation to test-anxiety relation persist, while several weaker connections depend on the correlation estimator. ```{r certify-spearman} certificate(spearman_net) ``` The Spearman fit has a residual of $1.55 \times 10^{-15}$ and is certified. Its numerical solution therefore satisfies the selected graph's defining conditions to numerical precision. ## Reporting the analysis A report should state the variables, sample size, correlation estimator, missing-data procedure, selection criterion, value of $\gamma$, and whether the stepwise refinement was used. The edge table should be reported numerically. Interpretation should distinguish partial correlations from ordinary correlations and from causal effects. The present analysis used Pearson correlations, ordinary BIC ($\gamma = 0$), and the default single-edge stepwise refinement on 300 complete observations. The selected network contained 9 edges. The largest positive edge was `CSU`-`IV` ($r = 0.422$), and the negative `SR`-`TA` edge was -0.370. The maximum-likelihood refit had a KKT residual of $1.33 \times 10^{-15}$. A Spearman sensitivity model retained the main pattern but selected 7 edges. ## How stepwise selection works The fitting procedure begins by generating candidate network structures across a graphical-lasso penalty path. At this stage, the path is a practical way to produce networks ranging from sparse to dense. Each distinct edge pattern is a candidate structure. Every candidate structure is then refitted by unregularized maximum likelihood. The Bayesian information criterion scores the tradeoff between fit and the number of edges. The candidate with the smallest score becomes the starting network for a local stepwise refinement. The refinement tests single-edge additions and removals, retaining a change whenever it improves the criterion. It stops when no single-edge change improves the score. The resulting graph is a local optimum among the one-edge alternatives examined by the search. A different correlation estimator or model-selection setting can lead to a different structure, as the Spearman analysis showed. The selected network should therefore be treated as an estimated model whose stability needs evaluation in a full study. ## Mathematical foundations This section states the model and selection criterion precisely. It can be skipped when the practical workflow is the main interest. ### Gaussian graphical model Let $\mathbf{X}$ follow a multivariate normal distribution with covariance matrix $\boldsymbol{\Sigma}$ and precision matrix $\mathbf{K} = \boldsymbol{\Sigma}^{-1}$. For $i \ne j$, $$ K_{ij} = 0 \quad \Longleftrightarrow \quad X_i \perp\!\!\!\perp X_j \mid \mathbf{X}_{\setminus\{i,j\}}. $$ The reported partial correlation is $$ \rho_{ij \cdot \mathrm{rest}} = -\frac{K_{ij}}{\sqrt{K_{ii}K_{jj}}}. $$ ### Information criterion For a candidate graph with $E$ edges, sample size $n$, and $p$ nodes, the extended Bayesian information criterion is $$ \operatorname{EBIC}_{\gamma} = -2\ell(\widehat{\mathbf{K}}) + E\log(n) + 4E\gamma\log(p). $$ The default $\gamma = 0$ removes the final term and yields ordinary BIC. The selected graph is the candidate with the smallest criterion value, followed by the single-edge refinement. ### Constrained maximum-likelihood certificate Let $\mathbf{S}$ be the sample correlation matrix and $\mathbf{W} = \widehat{\mathbf{K}}^{-1}$ be the covariance matrix implied by the fitted precision matrix. For every retained edge and every diagonal entry, the constrained optimum satisfies the corresponding covariance equality between $\mathbf{W}$ and $\mathbf{S}$. Precision entries outside the selected graph are zero. `certificate()` reports the maximum absolute violation across these conditions. A value near machine precision indicates that the constrained fit satisfies them to numerical precision. ### Predictability For node $j$, the conditional residual variance is $1/K_{jj}$. Its explained variance is $$ R_j^2 = 1 - \frac{1}{K_{jj}S_{jj}}, $$ where $S_{jj}$ is the marginal variance. `net_predict()` evaluates this identity for each node in the selected Gaussian graphical model.