--- title: "The Nestimate Compatibility Map" author: - name: "Mohammed Saqr" affiliation: "University of Eastern Finland" - name: "Sonsoles López-Pernas" affiliation: "University of Eastern Finland" - name: "Manuel J. Gómez" affiliation: "University of Murcia" output: rmarkdown::html_vignette: toc: true number_sections: true vignette: > %\VignetteIndexEntry{The Nestimate Compatibility Map} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.align = "center", fig.width = 7, fig.height = 4.5, dpi = 110, out.width = "100%" ) library(tsn) has_nestimate <- requireNamespace("Nestimate", quietly = TRUE) has_cograph <- requireNamespace("cograph", quietly = TRUE) && requireNamespace("igraph", quietly = TRUE) ``` # One object, two packages A `ts_tna()` model is not exported *to* Nestimate; it *is* a Nestimate object. This vignette maps what that buys, verb by verb: which parts of Nestimate's analysis surface accept a tsn transition model directly, what question each verb answers, what each result draws, and which parts do not apply and why. The companion vignette `vignette("nestimate-workflow", package = "tsn")` runs a single analysis end to end; this one is the reference you consult when deciding which verb to reach for. The running example is the model from that workflow: the 36 students of the `srl` daily panel, whose `effort` reports are cut at 40 and 70 on the 0-100 scale so the states mean something outside the sample. ```{r build, eval=has_nestimate} library(tsn) library(Nestimate) data(srl) students <- subset(srl, !is.na(effort)) effort_model <- ts_tna( students, value = "effort", id = "name", time = "day", discretization = "threshold", breaks = c(40, 70), labels = c("low", "moderate", "high") ) class(effort_model) ``` The class chain is the whole compatibility story. `ts_tna` is the tsn layer (it remembers the source series and discretization), `netobject` is what every Nestimate verb dispatches on, and `cograph_network` is what `cograph::splot()` renders. Nestimate's own validator agrees: ```{r validate, eval=has_nestimate} validate_netobject(effort_model) ``` And because the object is a `cograph_network`, the model itself is one `plot()` away: ```{r plot-model, eval=has_nestimate && has_cograph, fig.width=6.5, fig.height=5, fig.alt="Directed transition network with low, moderate, and high effort nodes; edge labels show transition probabilities, with the heaviest self-loop on the high state."} plot(effort_model, type = "network") ``` The strong diagonal is the first substantive fact: every state's most likely successor is itself, `high` most strongly of all at 0.575. # Reading the model as data Both packages provide accessors, and they answer slightly different questions. tsn's `as.data.frame()` returns every state pair including self-transitions; `as.matrix()` returns the row-stochastic weight matrix. ```{r tsn-accessors, eval=has_nestimate} head(as.data.frame(effort_model)) ``` Nestimate's extractors return the same estimates in its own shapes: ```{r extract, eval=has_nestimate} extract_transition_matrix(effort_model) extract_initial_probs(effort_model) extract_edges(effort_model) ``` One difference matters for transition networks: `extract_edges()` **excludes self-transitions**, returning the six between-state edges sorted by weight and omitting the diagonal. For this model the diagonal is where most of the probability lives (`high → high` is 0.575), so use `as.data.frame()` or `extract_transition_matrix()` whenever persistence is part of the question, and `extract_edges()` when only between-state movement is. The initial-state probabilities say the panel starts engaged: 41.7% of students' first recorded day is a `high`-effort day, 19.4% a `low` one. # Where the process spends its time `state_distribution()` and `state_frequencies()` count observations per state, the marginal composition before any dynamics. ```{r frequencies, eval=has_nestimate} state_frequencies(effort_model) ``` High-effort days make up 41.7% of the panel and low-effort days 24.4%. Three displays draw the composition from different angles. The mosaic is the proportion bar: ```{r mosaic, eval=has_nestimate, fig.width=6.5, fig.height=4, fig.alt="Mosaic plot of the effort model's state composition, with high and moderate days forming most of the area and low days a narrower band."} mosaic_plot(effort_model) ``` `sequence_plot()` is the sequence carpet, one row per student, one coloured cell per day: ```{r carpet, eval=has_nestimate, fig.width=7, fig.height=5, fig.alt="Sequence heatmap with one row per student and one coloured cell per day, showing same-state spells as horizontal streaks of colour."} sequence_plot(effort_model) ``` The streaks are persistence made visible: students hold a state for spells, not single days. `distribution_plot()` stacks the day-by-day state mix of the whole panel: ```{r composition-time, eval=has_nestimate, fig.width=7, fig.height=4, fig.alt="Stacked area chart of the panel's state proportions across the 156 course days, with bands for low, moderate, and high effort."} distribution_plot(effort_model) ``` This is composition over course time, the view in which a first-to-second-half shift would appear, and the workflow vignette's period comparison shows there is none. `plot_state_frequencies()` offers a fourth, marimekko-style variant of the same information. # Which states and edges carry the structure? `net_centrality()` ranks states. Self-transitions are excluded by default, because for a transition network the diagonal would otherwise dominate every strength measure. ```{r centrality, eval=has_nestimate} effort_centrality <- net_centrality(effort_model) effort_centrality ``` ```{r centrality-plot, eval=has_nestimate, fig.width=7, fig.height=4, fig.alt="Centrality chart of the three states showing high and moderate with similar in-strength and low well behind."} plot(effort_centrality) ``` `high` and `moderate` receive similar between-state traffic (in-strengths of 0.60 and 0.57) while `low` receives markedly less (0.38); the betweenness column is uniformly zero because in a fully connected three-state network no state is an obligatory bridge. `net_edge_betweenness()` asks the same question of edges: ```{r betweenness, eval=has_nestimate} effort_betweenness <- net_edge_betweenness(effort_model) effort_betweenness ``` ```{r betweenness-plot, eval=has_nestimate, fig.width=6.5, fig.height=4.5, fig.alt="Edge-betweenness network view in which every between-state edge carries the same unit load."} plot(effort_betweenness) ``` Every between-state edge carries the same load of 1: each lies on exactly one shortest path, and none is a privileged route. That uniformity is itself a description: movement between effort states does not funnel through any single corridor. # How the process moves ## Mean first passage times `passage_time()` converts the transition matrix into expected travel times between states (Kemeny & Snell, 1976). ```{r passage, eval=has_nestimate} effort_passage <- passage_time(effort_model) effort_passage ``` ```{r passage-plot, eval=has_nestimate, fig.width=6.5, fig.height=4.5, fig.alt="Heatmap of mean first passage times between the three states, with the longest times into the low state."} plot(effort_passage) ``` A `moderate` or `high` day is reached from anywhere within about two to four days on average; a `low` day takes four to five and a half. Low effort is the state this panel visits reluctantly, not a basin it falls into. The stationary distribution (24% low, 34% moderate, 42% high) is where the chain settles regardless of its start. ## Persistence, return, and sojourn `markov_stability()` summarizes each state's dynamics in one row: the self-transition probability, the stationary share, the mean recurrence time, and the mean sojourn (how long a visit lasts once entered). ```{r markov-stability, eval=has_nestimate} effort_spells <- markov_stability(effort_model) effort_spells ``` ```{r markov-stability-plot, eval=has_nestimate, fig.width=7, fig.height=4, fig.alt="Bar panels of persistence, stationary probability, return time, and sojourn time for the three states."} plot(effort_spells) ``` A high-effort spell lasts 2.4 days on average and recurs every 2.4; a low-effort spell lasts 1.7 days and recurs only every 4.1. Persistence and prevalence line up here, but they are logically distinct columns: a rare state can still be sticky. ## Is the chain well-behaved? `chain_structure()` reports the Markov-chain diagnostics that other verbs quietly assume: irreducibility (every state reachable from every other), aperiodicity, and reversibility. ```{r chain, eval=has_nestimate} summary(chain_structure(effort_model)) ``` The chain is irreducible, aperiodic, and regular, so the stationary distribution above is well-defined and unique. It is *not* reversible: the forward process is statistically distinguishable from its time reversal, which is itself a substantive fact about self-regulation: the routes up the effort scale and the routes down are not mirror images. # How predictable is the process? `transition_entropy()` gives the point estimate; `entropy_bayes()` adds a posterior credible interval by placing a Dirichlet prior on each row. ```{r entropy-bayes, eval=has_nestimate} set.seed(2026) effort_entropy <- entropy_bayes(effort_model) effort_entropy ``` ```{r entropy-bayes-plot, eval=has_nestimate, fig.width=7, fig.height=4, fig.alt="Posterior distributions of the entropy rate, stationary entropy, and redundancy, each a narrow density around its point estimate."} plot(effort_entropy) ``` The entropy rate is 1.478 bits against a ceiling of `log2(3) = 1.585`, or 0.93 on the normalised scale, with a credible interval three hundredths wide: with 5,600 daily reports, the *uncertainty* about how unpredictable the process is, is itself small. The redundancy of 0.073 bits is the gap between the marginal entropy and the conditional one: knowing today's effort removes under 5% of tomorrow's uncertainty. Daily effort is mostly noise around a weak first-order signal, and every edge below should be read with that in mind. `entropy_network()` maps the same quantity onto the edges, each entry an edge's contribution to the row's unpredictability, and returns a network of the same class, so it can be plotted or pruned like any other: ```{r entropy-network, eval=has_nestimate} entropy_network(effort_model) ``` # Does history beyond one step matter? The model assumes the next state depends only on the current one. `markov_order_test()` (Anderson & Goodman, 1957) tests that assumption. ```{r markov-order, eval=has_nestimate} set.seed(2026) effort_order <- markov_order_test(effort_model) effort_order ``` ```{r markov-order-plot, eval=has_nestimate, fig.width=7, fig.height=4, fig.alt="Model-selection chart of log-likelihood, AIC, and BIC across Markov orders zero to three, favouring orders above one."} plot(effort_order) ``` AIC and the permutation LRT select order 3 and BIC selects order 2: yesterday carries information about tomorrow beyond today. `path_dependence()` says how much: ```{r path-dependence, eval=has_nestimate} effort_history <- path_dependence(effort_model) effort_history ``` ```{r path-dependence-plot, eval=has_nestimate, fig.width=7, fig.height=4, fig.alt="Bar chart of the KL divergence contributed by each two-step context, all values below 0.12 bits."} plot(effort_history) ``` The weighted divergence is 0.038 bits: history shifts the next-step distribution by a few hundredths of a bit and flips the modal prediction in only two of nine contexts. This is the pair of numbers that reconciles a decisively significant order test with the continued use of a first-order model: detectable is not the same as large. The first-order network remains a fair one-step summary; it is not the full process. `path_counts()` is the raw material, every observed two-step path with its count: ```{r path-counts, eval=has_nestimate} path_counts(effort_model) ``` `high → high` alone is 24.0% of all transitions. And when the order test demands more than one step, `build_hon()` estimates a higher-order network under its own per-node promotion rule: ```{r hon, eval=has_nestimate} hon <- build_hon(effort_model, max_order = 2) pathways(hon) ``` HON promotes no node here, because no single context shifts the next-step distribution enough to justify splitting a state, which is why `pathways()` comes back empty. Two criteria disagreeing is a result, not a bug: the likelihood test sees aggregate improvement, HON asks a stricter per-node question. Report both. # Simplifying the model `net_prune()` removes weak edges and returns the same class; `net_pruning_details()` is the audit trail; `net_deprune()` undoes it. ```{r prune, eval=has_nestimate} pruned <- net_prune(effort_model, method = "threshold", threshold = 0.3) net_pruning_details(pruned) ``` Four edges fall, and they are a coherent four: every descending route (`moderate → low`, `high → low`, `high → moderate`) plus the direct `low → high` leap. What survives is the ascending ladder and the persistence loops. Because pruning is recorded rather than destructive, the original model is one call away: ```{r deprune, eval=has_nestimate} net_deprune(pruned) ``` # How sure are we about the estimates? Nestimate offers several answers resting on different machinery: sequence resampling, node resampling, case dropping, sample splitting, and a posterior. ## Sequence bootstrap `bootstrap_network()` resamples the 36 sequences and re-estimates. ```{r bootstrap, eval=has_nestimate} effort_boot <- bootstrap_network(effort_model, iter = 200, seed = 2026) effort_boot ``` `summary()` is the tidy one-row-per-edge table: ```{r bootstrap-summary, eval=has_nestimate} head(summary(effort_boot), 3) ``` Five of nine edges clear the stability criterion: the `moderate` and `high` persistence loops, the two ascending steps `low -> moderate` and `moderate -> high`, and the single short descent `high -> moderate`. The four that fail are the `low -> high` leap, the two steep descents into `low` (`moderate -> low`, `high -> low`), and the `low -> low` loop, which just misses (p = 0.060) because `low` is the rarest state and its outgoing row rests on the fewest transitions. This is not quite the set that pruning kept: pruning discarded `high -> moderate` as too small and held the large `low -> low` loop, so ranking edges by size and ranking them by stability trade one edge in each direction, the caution the pruning section named. ## Bayesian edge certainty `certainty()` reaches its verdict from a Dirichlet posterior on each row, with no resampling. ```{r certainty, eval=has_nestimate} set.seed(2026) certainty(effort_model) ``` All nine edges are posterior-certain, with intervals a few hundredths wide. The two verbs disagree because they ask different questions: `certainty()` asks whether an edge's posterior excludes instability given the pooled transitions, while the bootstrap asks whether the estimate survives swapping *students* in and out. An edge can be precisely estimated from the pooled data and still depend on who is in the panel, which is exactly what befalls the four edges the bootstrap could not confirm. ## Whole-network statistics: the vertex bootstrap `vertex_bootstrap()` implements the Snijders and Borgatti (1999) vertex bootstrap for network-level statistics. ```{r vertex-boot, eval=has_nestimate} vertex_bootstrap(effort_model, iter = 100, seed = 2026) ``` The mean edge weight carries an interval of [0.23, 0.30] because the vertex bootstrap resamples *nodes* and there are only three, a deliberately conservative design for small node sets. ## Case dropping, in two flavours `centrality_stability()` tracks the centrality *ordering* as sequences are dropped; `casedrop_reliability()` tracks the edge weights themselves. Both report the CS-coefficient of Epskamp, Borsboom & Fried (2018). ```{r stability, eval=has_nestimate} set.seed(2026) centrality_stability(effort_model, iter = 100) ``` ```{r casedrop, eval=has_nestimate} effort_reliability <- casedrop_reliability(effort_model, iter = 100, seed = 2026) effort_reliability ``` ```{r casedrop-plot, eval=has_nestimate, fig.width=7, fig.height=4, fig.alt="Case-dropping stability curve showing the correlation between full-sample and subsampled edge weights declining as larger fractions of the 36 sequences are dropped."} plot(effort_reliability) ``` The verdicts are moderate across the board: the in-strength ordering holds to a 0.60 drop, the out-strength ordering only to 0.40, and the edge weights sit exactly at the 0.50 bar. This is what 36 sequences buy, and it is the sample-size counterweight to the posterior certainty above: report orderings and intervals, not bare point values. ## Split-half reliability ```{r reliability, eval=has_nestimate} set.seed(2026) network_reliability(effort_model) ``` Randomly halving the panel 1,000 times yields two models that correlate at 0.85 on average with a mean absolute edge difference of 0.06: real agreement, with real room for movement. All four verbs in this section tell the same story at different grain: the backbone of the model is solid, the weaker edges are not yet pinned down. # Comparing two models Any two models on the same alphabet can be compared. The pair here splits the panel at occasion 78, the first and second halves of the course, with the same fixed `breaks`, so both models sit on one alphabet, and the same 36 students on both sides. ```{r halves, eval=has_nestimate} first_half <- ts_tna( subset(students, day <= 78), value = "effort", id = "name", time = "day", discretization = "threshold", breaks = c(40, 70), labels = c("low", "moderate", "high") ) second_half <- ts_tna( subset(students, day > 78), value = "effort", id = "name", time = "day", discretization = "threshold", breaks = c(40, 70), labels = c("low", "moderate", "high") ) ``` ## The difference itself `subtract_networks()` returns the edge-by-edge difference as an object cograph can draw: ```{r subtract, eval=has_nestimate} half_difference <- subtract_networks(first_half, second_half) half_difference ``` ```{r subtract-plot, eval=has_nestimate && has_cograph, fig.width=6.5, fig.height=5, fig.alt="Difference network of the two course halves, edges coloured by the sign of the change and all magnitudes small."} cograph::plot_difference(half_difference) ``` The largest gap on any edge is 0.033 (`low → moderate`). ## Four more verdicts on the same pair `compare_model()` condenses the pair into a metric panel; `bayes_compare()` asks whether each edge difference is credibly nonzero; `vertex_compare()` tests the network-level statistics; and the workflow vignette's `permutation()` run supplies the sequence-level test. ```{r compare, eval=has_nestimate} half_panel <- compare_model(first_half, second_half) half_panel ``` ```{r compare-plot, eval=has_nestimate, fig.width=7, fig.height=4.5, fig.alt="Comparison panel chart of the two half-course models across weight-deviation, correlation, and similarity metrics, all indicating near-identity."} plot(half_panel) ``` ```{r bayes-compare, eval=has_nestimate} set.seed(2026) bayes_compare(first_half, second_half) ``` ```{r vertex-compare, eval=has_nestimate} vertex_compare(first_half, second_half, iter = 200, seed = 2026) ``` Five verbs, one verdict. The matrices correlate at 0.990 with rank and sign agreement of 1; the posterior finds zero of nine edges credibly different; the vertex comparison finds no network-level statistic moved (all p ≥ 0.67); and the permutation test in the workflow vignette finds nothing at the sequence level either. The course's second half reproduces its first, in level and in dynamics. When the verbs *disagree*, a posterior flagging tiny-but-credible differences that a sequence-level permutation absorbs, report the magnitude next to whichever test is used; the grouped weekday/weekend comparison in `vignette("group-models", package = "tsn")` shows what a real, localized difference looks like on these same students. # Link prediction `predict_links()` scores unobserved dyads by six structural predictors: ```{r predict, eval=has_nestimate} predict_links(effort_model) ``` It runs, but a three-state transition network is fully connected between states, which leaves link prediction nothing to rank; it earns its place on larger, pruned alphabets, where a missing edge is a real hypothesis. # Grouped models When the comparison variable is a column of the data, `ts_tna(group = ...)` builds one network per level and returns a Nestimate `netobject_group`. The grouped verbs verified against it are `net_centrality()`, `state_distribution()`, `transition_entropy()`, `passage_time()`, `certainty()`, `net_prune()`, `compare_model()`, `permutation()`, `bootstrap_network()`, `markov_order_test()`, `centrality_stability()`, `casedrop_reliability()`, and `mosaic_plot()`; `summary()` and `as.data.frame()` return tidy per-group tables, and `plot()` draws one selected group. `net_edge_betweenness()` also runs, but its result drops back to a plain `netobject_group`. The dedicated vignette `vignette("group-models", package = "tsn")` works through all of it on the `esm_srl` weekday/weekend split. # What does not apply, and why Three families of Nestimate verbs reject a transition model, each for a principled reason rather than a missing feature: * **Sequence-data verbs.** `frequencies()`, `mark_first_state()`, `mark_terminal_state()`, `entropy_trajectory()`, `mosaic_analysis()`, `plot_mosaic()`, and `as_htna()` operate on a raw wide data frame of sequences, not on an estimated model. They belong *before* model construction in a workflow, applied to sequence data you hold yourself. * **Psychometric-network verbs.** `predictability()` and the `glasso`/`pcor`/`cor` estimator family need a precision or correlation structure estimated from multivariate data. A transition matrix is a conditional-probability object, not a conditional-independence one; the quantities these verbs compute are undefined for it. * **Verbs that need extra structure.** `cluster_network()` and `build_clusters()` require a chosen number of clusters and are aimed at multi-sequence typologies; `wtna()` needs one-hot action data; `sequence_compare()` wants sequence data plus a grouping vector. And `as_tna()` is unnecessary by design: a tsn model needs no conversion, which is the point of this vignette. # The compatibility map | Question | Verb | Tidy access | Plot | Notes | |---|---|---|---|---| | What are the estimates? | `as.data.frame()`, `as.matrix()`, `extract_transition_matrix()`, `extract_initial_probs()` | direct | via `plot()` on the model | tsn accessors include the diagonal | | Between-state edges only | `extract_edges()` | direct | - | **excludes self-transitions** | | Marginal composition | `state_distribution()`, `state_frequencies()` | data frame | `mosaic_plot()`, `plot_state_frequencies()` | | | Composition over time / carpet | `distribution_plot()`, `sequence_plot()` | - | draws directly | stacked area; sequence heatmap | | State/edge importance | `net_centrality()`, `net_edge_betweenness()` | data frame / network | `plot()` on both | loops excluded from centrality by default | | Travel times between states | `passage_time()` | print | `plot()` | Kemeny & Snell (1976) | | Persistence, recurrence, sojourn | `markov_stability()` | `$stability` table | `plot()` | | | Chain diagnostics | `chain_structure()` | `summary()` | `plot()` heatmap | irreducibility, reversibility | | Predictability | `transition_entropy()`, `entropy_bayes()`, `entropy_network()` | print / `summary()` | `plot()` on both entropy results | rate in bits; ceiling `log2(n_states)` | | History beyond one step | `markov_order_test()`, `path_dependence()`, `path_counts()`, `build_hon()` + `pathways()` | print / data frame | `plot()` on test and dependence | test the size, not only the p | | Simplification | `net_prune()`, `net_pruning_details()`, `net_deprune()` | same class back | via model `plot()` | prune after testing | | Edge uncertainty | `bootstrap_network()`, `certainty()` | `summary()` | none | resampling vs Bayesian; they can disagree | | Network-level uncertainty | `vertex_bootstrap()` | data frame | `plot()` | Snijders & Borgatti (1999) | | Sample-size reliability | `centrality_stability()`, `casedrop_reliability()`, `network_reliability()` | print | `plot()` on all three | CS-coefficient bar: 0.5 / 0.7 | | Two-model comparison | `permutation()`, `compare_model()`, `subtract_networks()`, `bayes_compare()`, `vertex_compare()` | `summary()` / data frame | `plot()`; `cograph::plot_difference()` | report magnitudes with tests | | Link prediction | `predict_links()` | print | none | needs a sparse alphabet to be informative | | Many-group comparison | `ts_tna(group=)` + grouped verbs | `as.data.frame()` | `plot(x, group=)`, `mosaic_plot()` | see `vignette("group-models")` | | Rendering | `plot()`, `cograph::splot()` | none | none | requires cograph | # References Anderson, T. W., & Goodman, L. A. (1957). Statistical inference about Markov chains. *The Annals of Mathematical Statistics*, 28(1), 89–110. Epskamp, S., Borsboom, D., & Fried, E. I. (2018). Estimating psychological networks and their accuracy: A tutorial paper. *Behavior Research Methods*, 50(1), 195–212. Kemeny, J. G., & Snell, J. L. (1976). *Finite Markov Chains*. Springer. Snijders, T. A. B., & Borgatti, S. P. (1999). Non-parametric standard errors and tests for network statistics. *Connections*, 22(2), 161–170. Xu, J., Wickramarathne, T. L., & Chawla, N. V. (2016). Representing higher-order dependencies in networks. *Science Advances*, 2(5), e1600028.