Simulation Workshop: Evaluating Graph Recovery And Edge Uncertainty

Why Evaluate Graph Estimators by Simulation?

A sparse graph estimator is useful only to the extent that it recovers relevant edges and estimates their strengths under conditions resembling the intended application. A single fitted graph cannot reveal how often a method misses true edges, adds false ones, or produces intervals that fail to cover. A simulation makes these questions measurable because the generating precision matrix, and therefore its graph, are known.

A sound comparison changes one or more design factors systematically while holding the rest fixed. Here the central factors are sample size and edge signal strength. The outcomes include graph support precision and recall, matrix estimation error, and coverage of edgewise confidence intervals. These quantities answer different questions; none alone is a complete ranking of methods.

Build a Controlled Data-Generating Model

The graph is encoded by the off-diagonal zeros of \(\Theta\), so a simulation can control graph structure directly by specifying \(\Theta\). It must also ensure that the matrix is positive definite. One convenient sufficient condition is strict diagonal dominance for a symmetric matrix with positive diagonal entries: if each diagonal entry exceeds the sum of the absolute off-diagonal entries in its row, the matrix is positive definite.

For a simple reproducible design, use a chain graph on \(p=50\) vertices. Set all diagonal entries to 2, and for each \(j=1,\ldots,49\), set \(\Theta_{j,j+1}=\Theta_{j+1,j}\) to \(0.4\) when \(j\) is odd and \(0.3\) when \(j\) is even. Set all other off-diagonal entries to zero. Multiply every nonzero off-diagonal entry by a signal-strength factor \(\gamma\in\{0.5,1,2\}\), leaving the diagonal fixed.

The largest absolute off-diagonal row sum before scaling is \(0.4+0.3=0.70\): an interior vertex between edges of those two sizes has that sum, and all other rows have a sum no larger than \(0.70\). At the largest signal level, the maximum is \(2(0.70)=1.40<2\). Thus the matrix is strictly diagonally dominant at every chosen signal level and is positive definite. This check matters: simulating from a purported precision matrix that is not positive definite does not define the intended Gaussian model.

For each design point, generate independent data sets from \(N_{50}(0,\Theta^{-1})\), using \(n\in\{50,100,200\}\). This gives nine combinations of sample size and signal strength. A larger study could add graph structures, such as hubs, but should report each structure separately rather than averaging unlike generating models into an opaque overall score.

Specify Estimation and Scoring Before Simulating

Compare at least two procedures, for example graphical lasso and nodewise lasso with a prespecified symmetrization rule. The comparison is only interpretable if tuning and edge inclusion are fully specified. One workable plan is to select each method's penalty using a prespecified extended BIC rule on the same candidate penalty grid, then apply a fixed coefficient threshold \(\tau\) to define \(\hat{E}\). For nodewise regressions, state whether an undirected edge is selected by the OR rule (either endpoint selects the other) or the AND rule (both do). Do not choose the tuning rule after viewing which method performs best at a particular design point.

For each replicate, record true positives \(TP=|E\cap\hat{E}|\), false positives \(FP=|\hat{E}\setminus E|\), and false negatives \(FN=|E\setminus\hat{E}|\). Then compute support precision and recall. Precision is the fraction of estimated edges that are true; recall is the fraction of true edges that are recovered. If no edges are estimated, define precision as zero and flag that convention in the report. If the true graph has no edges, recall is undefined; this chain design avoids that case.

For matrix error, one useful scale-adjusted measure is normalized Frobenius error, \(\|\hat{\Theta}-\Theta\|_F/\sqrt{p}\). It measures errors in all entries, not just mistakes in edge support; a separate off-diagonal error can be useful when diagonal accuracy is not the focus. Report which version is used. A method can have good support recall but poor coefficient estimates, or small matrix error while missing weak edges.

Evaluate Edgewise Interval Uncertainty

Graph selection and interval estimation are related but distinct tasks. For a concrete marginal interval assessment, one can prespecify a debiased precision estimator. Let \(S=n^{-1}\sum_{i=1}^n(X_i-\bar{X})(X_i-\bar{X})^\mathsf{T}\) be the sample covariance matrix, and let \(\hat{\Theta}\) be a symmetric initial sparse estimate. A commonly used debiasing step is

Under suitable high-dimensional regularity conditions, the off-diagonal entry \(\tilde{\Theta}_{jk}\) is approximately normal around \(\Theta_{jk}\). For Gaussian observations, the asymptotic variance for an off-diagonal precision entry has leading form \(\Theta_{jj}\Theta_{kk}+\Theta_{jk}^2\). A plug-in Wald interval is therefore

This is a marginal, nominal 95% interval, not a simultaneous confidence band for every edge. Its theoretical justification depends on the regime and regularity conditions; simulation coverage is an empirical diagnostic, not a proof that those conditions hold. To avoid selection-induced distortion, evaluate intervals for prespecified pairs, or for all pairs with true edges and null pairs separately, rather than reporting coverage only among edges selected from the same data.

For a set \(A\) of pairs, empirical coverage is the proportion of simulation replicates in which the interval contains the true coefficient. Report coverage for true edges and null edges as distinct summaries. For null edges the target is zero; for true edges the target is their nonzero coefficient. Also report average interval length, because very wide intervals can attain high coverage while conveying little precision.

Aggregate Results With Monte Carlo Uncertainty

Use \(R\) independent replicates at each combination of \(n\) and \(\gamma\), and retain replicate-level outcomes. For any performance statistic \(T_r\) from replicate \(r\), the average \(\bar{T}=R^{-1}\sum_{r=1}^R T_r\) estimates the expected performance under that fixed design. The following result explains why replication averages are appropriate.

Proof. By linearity of expectation,

where identical distribution gives \(E[T_r]=E[T_1]\). Independence makes all cross-covariances zero, so

The sample standard deviation of the \(T_r\)'s divided by \(\sqrt{R}\) estimates the Monte Carlo standard error of the mean. For coverage, each indicator is Bernoulli, giving the estimated standard error \(\sqrt{\hat{q}(1-\hat{q})/R}\), as used above. More replicates reduce simulation noise; they do not repair a poorly chosen data-generating model or an unfair tuning comparison. \(\square\)

For each method, show results by sample size and signal strength, preferably in a table or plot with Monte Carlo error bars. Do not pool all nine conditions into one average: a method that excels at strong signals and large samples may fail badly for weak signals and small samples. Include the number of replicates, the tuning and symmetrization rules, the edge threshold, the interval procedure, and the definition of every reported metric.

Common Pitfalls and Interpretation

A common mistake is to treat a high recall value as proof of a good graph. A procedure that selects nearly every possible pair can recover almost all true edges, hence high recall, while having very low precision. Conversely, a method that selects only a few highly convincing edges can have high precision but low recall. Report both, alongside estimation error and uncertainty results.

A second mistake is to report confidence-interval coverage only for edges selected by the same data. Selection makes the reported set data-dependent, so ordinary marginal coverage need not persist after conditioning on selection. The correction is to define the target set in advance, report true-edge and null-pair coverage separately, and state clearly whether any selection adjustment or independent sample split is used. Likewise, marginal 95% coverage does not mean that all intervals jointly cover their targets with probability 0.95.

Finally, verify the generating model numerically as well as conceptually. Check symmetry, positive definiteness, and the intended edge support before drawing samples. For the chain design above, the largest base row sum is \(0.70\), and at \(\gamma=2\) it is \(1.40\), still less than the diagonal value 2. Rechecking such simple quantities prevents a simulation from answering a different question than the one intended.