Two DEMATEL models can have the same number of factors, the same rating scale and the same panel size, and still require opposite interventions. Which one you are holding is a property of the matrix’s spectrum rather than of its subject matter, so you cannot look up your intervention logic, you have to diagnose it.
This vignette works through that diagnosis on the two published matrices the package ships.
library(spectralDEMATEL)
fefo <- worked_matrices$fefo_stock_control # 13 factors
resl <- worked_matrices$resilience_capabilities # 15 factorsBefore any number means anything, the matrix has to satisfy the conditions the diagnostics rest on. Every one of them comes back as data — nothing here stops, warns or prints.
ck <- assumption_checks(fefo)
ck[ck$verdict != "pass", c("check", "verdict", "value")]
#> check verdict value
#> 6 zero_diagonal warn 13One warning, and it is worth having. This matrix has a non-zero diagonal: thirteen factors carry a self-influence. That is admissible and everything is computed as usual, but self-influence inflates coupling, so the value is not comparable with a study that set its diagonal to zero.
The second matrix is a harder case.
ck2 <- assumption_checks(resl)
ck2[ck2$verdict != "pass", c("check", "verdict", "value")]
#> check verdict value
#> 7 strong_connectivity fail 4
#> 8 connectivity_margin skipped NAStrong connectivity fails. A binary verdict would stop there; the point of returning the checks as data is that they can say which factors:
The chain of cause is visible in the matrix itself. Row 14 is blank, so factor 14 dispatches nothing. Column 7 is blank, so factor 7 receives nothing. Factor 8’s only source is 7, and factor 4’s only sources are 7 and 8 — so neither can be reached from the rest of the system either.
which(rowSums(resl) == 0) # dispatches nothing
#> [1] 14
which(colSums(resl) == 0) # receives nothing
#> [1] 7That is usually a coding slip in one row rather than a property of the system, which is exactly why naming the factors matters: a user is told where to look.
The engine computes the diagnostics anyway. This is deliberate — the published study reports them, and an engine that refused would make the check unimplementable. But the numbers carry the verdict with them, and here they carry a visible mark of it:
d_resl <- spectral_diagnostics(resl)
round(d_resl$entry_points[stranded], 4) # the top-ranked entry point is here
#> [1] 0.2851 0.2813 0.5391 0.0000
round(d_resl$accumulation[stranded], 4) # yet nothing accumulates at them
#> [1] 0.0000 0.0000 0.0000 0.3913The highest-ranked entry point in the whole system is one of the four factors cut off from it, and the left eigenvector puts exactly zero weight on three of them. For a strongly connected matrix the entry profile is unique and strictly positive; when A2 fails it is neither, and a ranking read off it is describing something other than the system as a whole.
Here is the mechanism on a matrix small enough to read. Five factors, densely and unevenly connected:
core <- matrix(c(0, 4, 3, 1, 1,
4, 0, 3, 1, 1,
2, 2, 0, 1, 1,
1, 1, 1, 0, 1,
1, 1, 1, 1, 0), 5, 5, byrow = TRUE)
storage.mode(core) <- "double"
d_core <- spectral_diagnostics(core)
c(irreducible = is_irreducible(core),
coupling = round(d_core$mu_max, 3), hierarchy = round(d_core$hierarchy_sd, 3))
#> irreducible coupling hierarchy
#> 1.000 0.778 0.162Coupling above 0.5 and hierarchy above 0.10: strongly coupled, with influence entering at a few factors.
Now add a sixth factor that dispatches heavily to all five and receives nothing — an observer, a reporting layer, a context variable someone added late.
observed <- rbind(cbind(core, 0), 0)
observed[6, 1:5] <- 4
d_obs <- spectral_diagnostics(observed)
c(irreducible = is_irreducible(observed),
coupling = round(d_obs$mu_max, 3), hierarchy = round(d_obs$hierarchy_sd, 3))
#> irreducible coupling hierarchy
#> 0.000 0.350 0.224Coupling has fallen from 0.78 to 0.35, crossing the 0.5 threshold, so
the system reads as dampened where it was amplified. Nothing about the
five original factors changed. What changed is the normalising constant:
dematel() divides by the larger of the maximum row and
column total, and factor 6’s row is now the largest, so the whole
normalised matrix shrank.
The check catches this. The verdict is not bureaucratic.
d_fefo <- spectral_diagnostics(fefo)
data.frame(
row.names = c("FEFO stock control", "Resilience capabilities"),
factors = c(d_fefo$n, d_resl$n),
coupling = c(d_fefo$mu_max, d_resl$mu_max),
multiplier = c(d_fefo$multiplier, d_resl$multiplier),
indirect = c(d_fefo$indirect_dominant, d_resl$indirect_dominant),
dominance = c(d_fefo$dominance, d_resl$dominance)
)
#> factors coupling multiplier indirect dominance
#> FEFO stock control 13 0.8278 5.808 TRUE 0.02789
#> Resilience capabilities 15 0.3512 1.541 FALSE 0.19937The two are opposites. FEFO stock control is strongly coupled: total effects run 5.8 times the direct ones, so indirect propagation dominates and an intervention has to be appraised on total effects rather than on the ratings as entered. Its mode dominance of 0.028 says one propagation pattern accounts for essentially all of the behaviour.
Resilience capabilities barely amplifies at all: a multiplier of 1.54 means indirect effects are a third of the direct ones. Direct effects are most of the story, and the ratings can be read close to face value.
Coupling, the dominant eigenvalue and the multiplier are one quantity on three scales, which the package asserts rather than assumes:
Hierarchy has three readings and they do not run the same way.
data.frame(
row.names = c("FEFO stock control", "Resilience capabilities"),
sd = c(d_fefo$hierarchy_sd, d_resl$hierarchy_sd),
gini = c(d_fefo$hierarchy_gini, d_resl$hierarchy_gini),
pr = c(d_fefo$hierarchy_pr, d_resl$hierarchy_pr)
)
#> sd gini pr
#> FEFO stock control 0.03776 0.07195 0.9420
#> Resilience capabilities 0.14776 0.37057 0.4522hierarchy_sd and hierarchy_gini are
high when influence enters at a few factors.
hierarchy_pr is low. Resilience
capabilities is the more concentrated of the two on all three readings,
which you can only see if you know that the third column runs
backwards.
Anything that displays one of these must say which reading it used
and which way it runs. This project has inverted the direction once
already, in code that had been reviewed, and the mistake produced
entirely plausible output. hierarchy_sd is the reading
behind every published corpus-level result.
The standard DEMATEL deliverable is prominence, the sum of a factor’s row and column totals in the total-relation matrix. It adds two different roles together: how much a factor dispatches, and how much it absorbs. The left and right eigenvectors separate them.
T_fefo <- dematel(fefo)$T
prominence <- rowSums(T_fefo) + colSums(T_fefo)
ranks <- data.frame(
factor = seq_len(d_fefo$n),
prominence = rank(-prominence),
entry = rank(-d_fefo$entry_points),
accumulation = rank(-d_fefo$accumulation)
)
head(ranks[order(ranks$entry), ], 5)
#> factor prominence entry accumulation
#> 5 5 1 1 5
#> 4 4 7 2 10
#> 3 3 5 3 9
#> 13 13 2 4 1
#> 12 12 3 5 3The three rankings disagree, and the disagreements are the interesting part.
Factor 4 is the second strongest entry point and only seventh by prominence. Factor 3 is third by entry and fifth by prominence. Both inject more influence than the standard deliverable credits them with.
Factor 13 runs the other way: second by prominence, first by accumulation, only fourth by entry. Prominence adds dispatch to absorption, so a factor that absorbs heavily while dispatching moderately can outrank the factor actually driving the system. Pressure applied at 13 lands where effects already accumulate, which is close to the opposite of intervening.
Only factor 5 tops all three, which is why a single blended score is persuasive right up to the point where it is wrong.
sensitivity_matrix() gives the derivative of the
dominant eigenvalue in every off-diagonal entry, so the question moves
from which factor matters to which relationship, if
strengthened or weakened, moves the system most.
s <- sensitivity_matrix(fefo)
strongest <- which(s$total == max(s$total, na.rm = TRUE), arr.ind = TRUE)
strongest
#> row col
#> [1,] 13 5Never read that ranking without the condition number.
ev_condition is at least 1 by construction and bounds
how far a first-order estimate can be trusted. At 1.01 these derivatives
can be read at face value; a value in the tens means the derivative is
locally uninformative, and assumption_checks() raises a
warning verdict past 5.
The derivative is signed, and some entries are negative:
Strengthening one of those lowers coupling — the lever for a system someone wants to calm rather than drive, which is the usual situation in risk work. The reason is not mysterious, and the function separates it out:
neg <- which(s$total < 0, arr.ind = TRUE)[1, ]
c(local = s$local[neg[1], neg[2]],
normalization = s$normalization[neg[1], neg[2]],
total = s$total[neg[1], neg[2]])
#> local normalization total
#> 0.08575 -0.87277 -0.78702The local effect is always an amplification. The normalisation effect is never one, and it is non-zero only in the row or column that currently sets the normalising constant. Where the second outweighs the first, raising an influence rating lowers the dominant eigenvalue. Shown as a single total, that reads as a bug; shown as two components, it is a mechanism.
A diagnosis is worth less if the same numbers would arise from the
rating distribution alone. surrogate_ensemble() shuffles
the off-diagonal entries, holding the number of factors, the density and
the exact multiset of values fixed, and recomputes.
ens <- surrogate_ensemble(fefo, B = 200, seed = 42)
c(observed = d_fefo$hierarchy_sd,
surrogate_median = median(ens$hierarchy_sd),
share_at_least_observed = mean(ens$hierarchy_sd >= d_fefo$hierarchy_sd))
#> observed surrogate_median share_at_least_observed
#> 0.03776 0.04834 0.86000Read that carefully, because it does not say what a null test usually says. The observed hierarchy is lower than the surrogate median: shuffling this matrix’s own ratings at random produces a more concentrated entry profile than the matrix actually has, about 86% of the time. The diffuseness of this system is not an artefact of its rating distribution — it is structure the analyst built, and it survives the comparison.
The same ensemble answers the question for coupling and dominance, and there the verdict is sharper still:
data.frame(
row.names = c("coupling", "dominance", "hierarchy_sd"),
observed = c(d_fefo$mu_max, d_fefo$dominance, d_fefo$hierarchy_sd),
surrogate_min = c(min(ens$mu_max), min(ens$dominance), min(ens$hierarchy_sd)),
surrogate_max = c(max(ens$mu_max), max(ens$dominance), max(ens$hierarchy_sd))
)
#> observed surrogate_min surrogate_max
#> coupling 0.82783 0.62943 0.82046
#> dominance 0.02789 0.03104 0.07792
#> hierarchy_sd 0.03776 0.01912 0.07584Coupling sits above every one of the 200 surrogates and mode dominance below every one of them. Shuffling this matrix’s ratings cannot reproduce either: the system is more strongly coupled, and more completely governed by a single propagation mode, than its own entry distribution can account for. Only hierarchy falls inside the ensemble, and on the low side.
That is the comparison the source paper could not run across its corpus, because the precomputed spectral file holds no raw matrices. A user of this package has their own matrix, so they can.
The seed is exposed so a figure can be reproduced. This is the only
function here that is not instant — each draw is an eigendecomposition
plus a connectivity check — so run it on demand rather than eagerly, and
lower B if a user is waiting.
A shuffle that disconnects the graph is rejected and redrawn. On a sparse matrix almost every shuffle disconnects it, so the redrawing can go on a long time. A directed cycle is the extreme case: its edges are exactly as many as strong connectivity needs, and essentially no rearrangement of them survives.
surrogate_ensemble() therefore stops after
max_attempts draws, 50 * B by default, and
returns what it has. Reaching that bound is not an error and does not
throw. It is recorded in the result.
attributes(ens)[c("requested", "attempts", "complete")]
#> $requested
#> [1] 200
#>
#> $attempts
#> [1] 200
#>
#> $complete
#> [1] TRUEcomplete is the one to read. A share computed over two
hundred draws and a share computed over seven are not the same claim,
and the row count alone does not say which you have.
if (isTRUE(attr(ens, "complete"))) {
mean(ens$mu_max >= d_fefo$mu_max)
} else {
"ensemble truncated; too few draws to quote a share"
}
#> [1] 0surrogate_position() passes the same three values
through, so the check is available whichever entry point you use. A
matrix that admits no valid shuffle at all returns NULL, as
spectral_diagnostics() does on inadmissible input.
Every diagnosis records the version of the engine that produced it.
A metric definition has moved in this project’s history before. When one moves again, a diagnosis recorded with its version stays interpretable rather than silently meaning something new. The checks travel in the same object for the same reason: numbers whose validity depends on something the caller was told separately are numbers whose validity will eventually be lost.