--- title: "Grouped Transition Models: One Network per Condition" 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{Grouped Transition Models: One Network per Condition} %\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) ``` # The question Does the same person regulate their effort differently on different kinds of day? The `esm_srl` study prompted 41 students up to a few times a day for momentary reports of self-regulation, motivation, and anxiety on 0-100 scales, and each report carries a `day_type` label: weekday or weekend. A single transition network estimated from all of it would answer how momentary `effort` moves between `low`, `mid`, and `high` *on average*, and in doing so would erase exactly the distinction the day-type column exists to draw. The `group` argument of `ts_tna()` and its siblings answers the conditional question instead: one network per day type, both cut from one shared state alphabet, packaged so that Nestimate's grouped verbs compare them directly. ```{r libs, eval=has_nestimate} library(tsn) library(Nestimate) data(esm_srl) ``` # One call, one network per condition `group` names a column of the data. Every level of that column becomes its own fully estimated `ts_tna` model. Effort has four missing values, so the incomplete rows are set aside first. ```{r grouped, eval=has_nestimate} effort_rows <- subset(esm_srl, !is.na(effort)) by_day <- ts_tna( effort_rows, value = "effort", id = "name", time = "occasion", group = "day_type", labels = c("low", "mid", "high") ) by_day ``` Two day types, two networks. `as.data.frame()` with `what = "groups"` is the tidy index of how much data stands behind each one, and it is where every grouped analysis should start. ```{r grouped-index, eval=has_nestimate} as.data.frame(by_day, what = "groups") ``` The weekday network rests on 2,033 observations in 238 sequences, the weekend network on 783 in 233. The `sequences` column is the one to read closely: a sequence is cut at every change of student *and* at every change of day type, so a working week of roughly two reports a day forms one sequence of eight to ten observations, while a weekend forms one of three to four. The two cut rules are not the same. A change of student separates measurements that were never consecutive; a change of day type falls between two genuinely adjacent measurements, and the transition across it, from Friday evening to Saturday morning, is excluded from both networks because it belongs to neither day type. The next section shows what goes wrong without that rule. # Why not subset the data and build each model by hand? Subsetting looks equivalent: keep the weekend rows, call `ts_tna()` on them, repeat for weekdays. It is not equivalent, and it fails in two distinct ways. ```{r weekend-alone, eval=has_nestimate} weekend_alone <- ts_tna( subset(effort_rows, day_type == "weekend"), value = "effort", id = "name", time = "occasion", labels = c("low", "mid", "high") ) head(summary(series_networks(weekend_alone)), 3) ``` ## First failure: the subset welds broken time back together The hand-built model holds each student's weekend reports as **one** sequence. The first student's 15 weekend reports above span many separate weekends, yet they form a single unbroken series. Across all 41 students the 783 weekend observations become 41 sequences carrying 742 transitions. The grouped index showed the truth: those observations form **233** separate weekend runs, which contain only 783 − 233 = 550 genuine transitions. The hand-built model therefore counts 192 transitions, about a quarter of its total, between measurements that were never adjacent in time. Each one joins the end of one weekend to the start of the next across an entire intervening week, and each is a transition that never occurred. The grouped model never counts them. A sequence is cut at every change of group, so a transition is only ever counted between genuinely consecutive measurements within one day type. ## Second failure: each subset learns its own thresholds The second failure is quieter. The default `discretization = "quantile"` cuts at the tertiles of whatever data it is given, so the hand-built weekend model learns cut points from weekend reports alone: ```{r weekend-alone-dist, eval=has_nestimate} summary( discretize( subset(effort_rows, day_type == "weekend"), value = "effort", id = "name", time = "occasion", labels = c("low", "mid", "high") ) ) ``` By construction, exactly a third of weekend reports land in each state (261 apiece), so the weekend composition can carry no information at all. The grouped model cuts the states from the **pooled** series before splitting, and its weekend column reports the real composition on the shared scale: ```{r grouped-dist, eval=has_nestimate} state_distribution(by_day) ``` On this panel the distortion happens to be mild, because weekday and weekend effort sit at nearly the same level (means of 58.5 and 59.2): the shared scale puts each weekend state within a percentage point or two of a third. That mildness is a property of these data, not of the method. For groups that differ in level, such as a clinical against a control cohort or a high- against a low-workload term, per-subset tertiles redefine the states inside each group, and the artifact grows with the very difference under study. The grouped model makes the safe construction the default; explicit `breaks` (as in `vignette("nestimate-workflow", package = "tsn")`) achieve the same discipline by hand. # The collection is data Every view of the grouped model is a data frame with a leading `group` column. The default view is the edge list: ```{r edges, eval=has_nestimate} head(as.data.frame(by_day), 6) ``` Because both groups share one alphabet, a single edge can be read across them: ```{r one-edge, eval=has_nestimate} subset(as.data.frame(by_day), from == "high" & to == "mid") ``` The probability that a high-effort report softens to `mid` is 0.254 on weekdays and 0.342 on weekends. This is the edge the inference section below puts to the test. `what = "series"` returns the per-observation table behind the networks, with the run-level sequence IDs, and `summary()` adds Nestimate's network-metric panel, one column per group: ```{r summary, eval=has_nestimate} summary(by_day) ``` Most rows are structural constants for a full three-state network, but the standard deviation of out-strength is not: 0.032 on weekdays against 0.102 on weekends. Weekday rows spread their probability almost evenly; weekend rows are the more polarized ones. # Composition and dynamics ## Composition The state distributions above showed both day types holding close to a third of reports in each state. The mosaic view draws the transition structure itself, one panel per group, with tile area encoding transition probability: ```{r mosaic, eval=has_nestimate, fig.width=8, fig.height=4.5, fig.alt="Mosaic panels of the weekday and weekend effort transition networks, in which tile area encodes transition probability; the weekend high-to-mid tile is visibly larger than its weekday counterpart."} mosaic_plot(by_day) ``` The two mosaics differ mainly at the `high` state: the weekend network routes more probability from `high` to `mid` (0.342 against 0.254) and retains less in `high` (0.534 against 0.631). This is the same softening the edge table reported, now shown as tile area rather than as a number. ## Mean first passage times `passage_time()` reports, per group, the expected number of steps to reach each state from each other state (Kemeny & Snell, 1976). ```{r passage, eval=has_nestimate} passage_time(by_day) ``` The travel times are broadly similar, since every state is reachable within three to seven steps on either day type, but the stationary distributions tilt: the weekday chain settles with 33.2% of its time in `high`, the weekend chain with 27.9%. Since the observed compositions are nearly identical, that tilt is carried by the dynamics, and it points at the same place the edge table did: weekends hold high effort less firmly. ## Entropy and centrality ```{r entropy, eval=has_nestimate} transition_entropy(by_day) ``` Both day types run near the `log2(3) = 1.585` ceiling, at 1.37 and 1.39 bits, so momentary effort is weakly predictable everywhere; neither day type is a regime of rigid habit. ```{r centrality, eval=has_nestimate} day_centrality <- net_centrality(by_day) day_centrality ``` ```{r centrality-plot, eval=has_nestimate, fig.width=7, fig.height=4.5, fig.alt="Faceted centrality chart comparing in-strength, betweenness, and diffusion of the three states across the weekday and weekend networks, with the mid state carrying the highest in-strength and all betweenness in both."} plot(day_centrality) ``` `mid` is the hub of both networks. It carries the highest in-strength (0.54 weekday, 0.62 weekend) and is the only state with nonzero betweenness, because movement between `low` and `high` routes through the middle of the scale. The message printed above the tables matters for any transition network: self-transitions are excluded by default, or the diagonal would drown the between-state structure. `net_edge_betweenness()` runs group-wise as well, with one caveat: its result drops back to a plain `netobject_group`, so apply tsn's tidy accessors before this step, not after it. ```{r betweenness, eval=has_nestimate} net_edge_betweenness(by_day) ``` # Drawing the networks There is no single picture of a grouped model; `plot()` draws one selected group with the full `plot.ts_tna()` surface. ```{r plot-weekday, eval=has_nestimate && has_cograph, fig.width=6, fig.height=5, fig.alt="Transition network of weekday effort: three nodes labelled low, mid, and high with strong self-loops, the high-to-high loop heaviest."} plot(by_day, group = "weekday", type = "network") ``` ```{r plot-weekend, eval=has_nestimate && has_cograph, fig.width=6, fig.height=5, fig.alt="Transition network of weekend effort: the same three nodes, with a lighter high-to-high self-loop and a heavier high-to-mid edge than the weekday network."} plot(by_day, group = "weekend", type = "network") ``` Because the two networks share one node set, the figures are directly comparable, and the difference to look for is concentrated at the `high` node: its self-loop is heavier on weekdays (0.631 against 0.534) and its outflow to `mid` heavier on weekends. # Are the differences real? ## The permutation test `permutation()` on a grouped model shuffles sequences between the groups to build a null distribution per edge. ```{r permutation, eval=has_nestimate} day_comparison <- permutation(by_day, iter = 1000, seed = 2026) summary(day_comparison) ``` One difference holds up under the null: `high -> mid` (0.254 weekday, 0.342 weekend, p = 0.023). Its mirror image, the fall in `high -> high` retention from 0.631 on weekdays to 0.534 on weekends, sits exactly on the significance boundary (p = 0.046) and moves around 0.05 under other seeds, so it should be read as suggestive rather than established. Together they tell one coherent story: on weekends, a high-effort spell is more likely to soften into the mid range and less likely to hold. The other seven edges do not differ, so weekday and weekend self-regulation share most of their structure. ## The metric panel `compare_model()` condenses the pair into matrix-comparison metrics. ```{r compare, eval=has_nestimate} compare_model(by_day, i = "weekday", j = "weekend") ``` A Pearson correlation of 0.966 between the two weight matrices, with a maximum absolute difference of 0.097 on the `high -> high` retention: close networks whose one real gap is how long high effort lasts. Always name `i` and `j`, because with more than two groups there is no meaningful default pair. ## Bootstrap and Bayesian certainty per group `bootstrap_network()` resamples sequences within each group; `certainty()` reaches a similar verdict from a Dirichlet posterior, without resampling. ```{r bootstrap, eval=has_nestimate} bootstrap_network(by_day, iter = 200, seed = 2026) ``` ```{r certainty, eval=has_nestimate} set.seed(2026) certainty(by_day) ``` The weekday network, with 2,033 observations behind it, keeps eight of nine bootstrap edges and all nine posterior-certain ones; the weekend network, with 783, keeps four and seven. Equal sequence counts (238 against 233) do not buy equal precision when the sequences differ in length: a weekend run is about a third as long as a weekday run, so the weekend estimates rest on less data. Reading the index table first is not a formality. ## Is first order enough, per group? `markov_order_test()` (after Anderson & Goodman, 1957) checks the Markov assumption group by group. ```{r markov, eval=has_nestimate} set.seed(2026) markov_order_test(by_day) ``` The weekday test prefers a higher order (BIC selects 2, AIC and the permutation LRT select 3): within a working week, yesterday's effort carries information beyond today's. The weekend test is more temperate, since BIC keeps order 1 and the LRT stops at 2, though its shorter runs give order-3 contexts little data in which to show up. Read both networks as one-step summaries, with the weekday summary the more visibly incomplete of the two. ## Stability under case dropping Two checks ask whether the results would survive a smaller sample: `centrality_stability()` for the centrality *orderings*, `casedrop_reliability()` for the edge weights themselves (both after Epskamp, Borsboom & Fried, 2018). ```{r stability, eval=has_nestimate} set.seed(2026) centrality_stability(by_day, iter = 100) ``` ```{r casedrop, eval=has_nestimate} day_reliability <- casedrop_reliability(by_day, iter = 100, seed = 2026) day_reliability ``` ```{r casedrop-plot, eval=has_nestimate, fig.width=7, fig.height=4, fig.alt="Case-dropping stability curves for the weekday and weekend networks, showing the correlation between full-sample and subsampled edge weights declining as larger fractions of sequences are dropped."} plot(day_reliability) ``` The verdicts are mixed, and the mixture is the point. The centrality orderings are serviceable (coefficients of 0.5 to 0.9), but the edge-weight CS-coefficients come back at 0.3 for weekdays and 0.5 for weekends, at or below the 0.5 bar for acceptable stability. Individual edge weights from this panel should be reported with their intervals, not as point values, and the permutation result above should carry the inferential weight. A stability check that returns an uncomfortable number has done its job. # Pruning the collection `net_prune()` applies group-wise and returns the same class, so a pruned collection flows back into every verb above. ```{r prune, eval=has_nestimate} net_prune(by_day, method = "threshold", threshold = 0.3) ``` Each day type keeps five of its nine edges at a 0.3 threshold. Prune after testing, not before: the `high -> mid` edge that carries the weekend difference sits at 0.254 on weekdays, so a premature prune would have deleted one side of the finding. # What to watch for Four failure modes, all of which announce themselves: * **A level with almost no data.** Every level of the grouping column gets a network, however little stands behind it. The `groups` index is the guard; read `sequences` and `observations` before edges. * **A group column that alternates too fast.** Every change of group cuts a sequence, so a column that flips at almost every observation leaves runs of length one and no transitions at all. In the extreme case, a group whose observations are *never* adjacent, the model warns that the group's network is empty. Day type is a well-behaved grouping variable precisely because it changes only twice a week. * **A group whose rows do not survive discretization.** Temporal discretizers such as `ordinal` drop tail rows that cannot fill an embedding window. If that removes a group entirely, the model warns that the group was lost to discretization. A group you excluded yourself via `series` is not warned about, because that exclusion was your decision. * **Comparing a grouped model with a hand-subset model.** The two are not estimates of the same thing, for both reasons in the section above. # When to use which * **`group`.** The condition already exists as a column, there may be more than two levels, and the comparison itself is the point. Pooled discretization and boundary-safe sequence cutting come built in. * **Two models plus `permutation(x, y)`.** The split is not a column but a derived rule (an occasion cutoff, an external roster), and there are exactly two sides. Fix the state alphabet yourself with explicit `breaks`, or the comparison inherits the threshold artifact. The companion vignette `vignette("nestimate-workflow", package = "tsn")` runs this design on the `srl` panel's course halves. * **`series_networks()`.** The unit of comparison is the individual student (one network per person), not a condition shared across students. * **A single pooled model.** No conditioning variable at all; the average dynamics are the question. # 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.