The sensitivity analysis reported in Böllenrücher et al. (in press) examines how sequence length affects pattern identification for dichotomous dyadic sequences. Five sequence lengths are considered: 30, 60, 90, 180, and 720 measurement points. The simulations focus on four patterns: A1, B2, D4, and E4. Following the notation in Table 2 of Böllenrücher et al. (in press), A1 denotes actor-partner, B2 partial actor only, D4 actor-partner on the main, actor only on the second, and E4 complete actor only.
This vignette has two purposes. First, it uses one simulated dyad to
illustrate how the univariate and bivariate dyadicMarkov
workflows behave at 30 and 90 measurement points. These worked examples
are illustrations. Second, the vignette reproduces the full
sequence-length analysis across all 1,000 simulated dyads and all five
sequence lengths.
The worked examples and the full sensitivity analysis use the
simulation data included with dyadicMarkov. Because the
bivariate method described in Böllenrücher et al.
(in press) is developed for states = 2, the
univariate worked example is also presented with states = 2
for consistency. The univariate method described in Bollenrücher et al. (2023) supports any integer
states ≥ 2, so the binary setting shown here is only one
possible application.
The package includes five simulation data sets for sequence lengths
30, 60, 90, 180, and 720. Each data set contains 1,000 simulated dyads.
For each dyad, four dichotomous sequences are available. To match the
package workflow notation, this vignette refers to them as
FM_V1, SM_V1, FM_V2, and
SM_V2.
utils::data("data_complete_30", package = "dyadicMarkov")
utils::data("data_complete_60", package = "dyadicMarkov")
utils::data("data_complete_90", package = "dyadicMarkov")
utils::data("data_complete_180", package = "dyadicMarkov")
utils::data("data_complete_720", package = "dyadicMarkov")The columns TM1, TM2, and so forth contain the ordered measurements.
The remaining columns identify the member, variable, simulated pattern,
dyad, and archived historical classification. More specifically,
caseSimulated records the interaction pattern used to
generate the sequence, whereas local stores the
classification produced by the historical analysis. The historical
terminology is retained in these fixed simulation datasets; current
dyadicMarkov output uses the current scientific
nomenclature. For consistency with the bivariate workflow, the two
archived simulation variables are relabelled V1 (main
variable) and V2 (second variable) below. The helper then
extracts one sequence for a given dyad, variable, and member.
simulation_variables <- setNames(
rev(sort(unique(as.character(data_complete_30$variable)))),
c("V1", "V2")
)
simulation_chain <- function(x, n, dyad, variable, member) {
measurement_columns <- paste0("TM", seq_len(n))
variable_code <- unname(simulation_variables[[variable]])
keep <- x$dyad == as.character(dyad) &
x$variable == variable_code &
x$members == as.character(member)
as.integer(unlist(
x[keep, measurement_columns, drop = FALSE],
use.names = FALSE
))
}
simulation_pattern <- function(x, dyad, variable, member) {
variable_code <- unname(simulation_variables[[variable]])
x[
x$dyad == as.character(dyad) &
x$variable == variable_code &
x$members == as.character(member),
"caseSimulated",
drop = TRUE
]
}The five data sets were simulated separately. Therefore, for example, a 30-measurement-point sequence is not the first 30 observations of a corresponding 90-measurement-point sequence.
For the univariate illustration, we select simulated dyad 241 and
analyze one of its variables at sequence lengths of 30 and 90
measurement points. As in the univariate workflow, the analyzed member
sequence is denoted FM and the partner sequence
SM. These two sequence lengths illustrate how the available
transition information and the identified pattern can change with
sequence length.
FM_30 <- simulation_chain(data_complete_30, 30L, 241L, "V2", 1L)
SM_30 <- simulation_chain(data_complete_30, 30L, 241L, "V2", 2L)
FM_90 <- simulation_chain(data_complete_90, 90L, 241L, "V2", 1L)
SM_90 <- simulation_chain(data_complete_90, 90L, 241L, "V2", 2L)The empirical transition counts show the amount of observed transition information available at each sequence length.
emp_uni_30 <- dyadicMarkov::countEmp(FM_30, SM_30, states = 2L)
emp_uni_90 <- dyadicMarkov::countEmp(FM_90, SM_90, states = 2L)
rbind(
`L = 30` = rowSums(emp_uni_30),
`L = 90` = rowSums(emp_uni_90)
)
#> FM1_SM1 FM1_SM2 FM2_SM1 FM2_SM2
#> L = 30 10 12 0 7
#> L = 90 19 29 20 21For the 30-measurement-point sequence, one previous dyadic state has no observed outgoing transition. For the 90-measurement-point sequence, all four previous dyadic states contribute observed transitions.
uni_30 <- dyadicMarkov::univariatePattern(
chainFM = FM_30,
chainSM = SM_30,
states = 2L,
alpha = 0.05
)
uni_90 <- dyadicMarkov::univariatePattern(
chainFM = FM_90,
chainSM = SM_90,
states = 2L,
alpha = 0.05
)
c(
`L = 30` = uni_30$pattern,
`L = 90` = uni_90$pattern
)
#> L = 30 L = 90
#> "IM (A0)" "APM (A1)"
c(
`L = 30` = simulation_pattern(
data_complete_30, 241L, "V2", 1L
),
`L = 90` = simulation_pattern(
data_complete_90, 241L, "V2", 1L
)
)
#> L = 30 L = 90
#> "A1" "A1"The 30-measurement-point sequence is identified as
IM (A0), whereas the 90-measurement-point sequence is
identified as APM (A1). For the 90-measurement-point
example, A1 is also the simulated pattern.
For the bivariate illustration, we use the same simulated dyad 241 at
sequence lengths of 30 and 90 measurement points. Following the
bivariate workflow notation, the analyzed sequence is
FM_V1, with SM_V1 on the main variable and
FM_V2 and SM_V2 on the second variable. As in
the univariate section, this is a worked example rather than the full
sensitivity analysis.
FM_V1_30 <- simulation_chain(data_complete_30, 30L, 241L, "V1", 1L)
SM_V1_30 <- simulation_chain(data_complete_30, 30L, 241L, "V1", 2L)
FM_V2_30 <- FM_30
SM_V2_30 <- SM_30
FM_V1_90 <- simulation_chain(data_complete_90, 90L, 241L, "V1", 1L)
SM_V1_90 <- simulation_chain(data_complete_90, 90L, 241L, "V1", 2L)
FM_V2_90 <- FM_90
SM_V2_90 <- SM_90
emp_bi_30 <- dyadicMarkov::countEmpBivariate(
chainFM_V1 = FM_V1_30,
chainSM_V1 = SM_V1_30,
chainFM_V2 = FM_V2_30,
chainSM_V2 = SM_V2_30,
states = 2L
)
emp_bi_90 <- dyadicMarkov::countEmpBivariate(
chainFM_V1 = FM_V1_90,
chainSM_V1 = SM_V1_90,
chainFM_V2 = FM_V2_90,
chainSM_V2 = SM_V2_90,
states = 2L
)The global bivariate analysis first identifies whether the sequence is trivial, univariate, partial bivariate, or complete bivariate.
case_bi_30 <- dyadicMarkov::bivariateCase(emp_bi_30, alpha = 0.05)
case_bi_90 <- dyadicMarkov::bivariateCase(emp_bi_90, alpha = 0.05)
c(
`L = 30` = case_bi_30$case,
`L = 90` = case_bi_90$case
)
#> L = 30 L = 90
#> "trivial" "complete"
c(
`L = 30` = simulation_pattern(
data_complete_30, 241L, "V1", 1L
),
`L = 90` = simulation_pattern(
data_complete_90, 241L, "V1", 1L
)
)
#> L = 30 L = 90
#> "D4" "D4"Both selected sequences were simulated under pattern D4. The 30-measurement-point sequence is classified as trivial. The 90-measurement-point sequence is classified as complete bivariate and therefore proceeds to the corresponding local pattern identification.
complete_bi_90 <- dyadicMarkov::completePattern(emp_bi_90)
complete_bi_90$pattern
#> [1] "actor-partner on the main, actor only on the second (D4)"The 90-measurement-point sequence selects pattern D4, recovering the pattern under which that sequence was simulated.
The worked examples above considered one simulated dyad. We now consider all 1,000 simulated dyads at each of the five sequence lengths used in Böllenrücher et al. (in press). For every dyad, the four member-variable combinations are classified, giving 20,000 classifications in total.
The code below is included for reproducibility rather than as part of
the usual dyadicMarkov workflow. It extracts each simulated
sequence, applies the same classification procedure illustrated above,
and compares the identified pattern with the pattern used to simulate
that sequence. Because the complete calculation is relatively expensive,
the chunk is not evaluated when the vignette is built. Readers
interested in reproducing the analysis can run the code locally.
The same complete reproduction is also implemented as an opt-in
extended test in the package source. From a source checkout, it can be
run by setting DYADICMARKOV_EXTENDED_TESTS=true and
executing devtools::test(filter = "extended-sensitivity").
The extended test also compares the reproduced sensitivity and
specificity values with those reported in Table 5 of Böllenrücher et al. (in press), using an
absolute tolerance of 0.01 to accommodate the small specificity
differences documented below.
# -------------------------------------------------------------------------
# Helper 1: recover the A-code from a univariate result.
#
# univariatePattern() returns labels such as:
# "IM (A0)"
# "APM (A1)"
# "AM (A2)"
# "PM (A3)"
# -------------------------------------------------------------------------
univariate_code <- function(x) {
code <- sub(
"^.*\\(([A-Z][0-9]+)\\)$",
"\\1",
x$pattern
)
if (identical(code, x$pattern)) {
stop("The univariate pattern code could not be extracted.")
}
code
}
# -------------------------------------------------------------------------
# Helper 2: recover the current pattern code selected by an AIC procedure.
#
# partialPattern() and completePattern() return:
# - $pattern: the selected pattern written out in full;
# - $aic: the candidate table, including the short matrix code.
#
# Matching these two fields gives codes such as B2, D4, or E4.
# -------------------------------------------------------------------------
selected_matrix_code <- function(x) {
i <- match(x$pattern, x$aic$pattern)
if (is.na(i)) {
stop("The selected pattern could not be matched to its AIC table.")
}
x$aic$matrix[i]
}
# -------------------------------------------------------------------------
# Helper 3: classify one selected member-variable combination.
#
# For each classification, the selected member is placed in the
# first-member position for the main variable, while the other member is
# placed in the second-member position. The second variable is supplied
# using the same member order.
#
# The global bivariate classification determines whether the result is:
# - trivial;
# - univariate;
# - partial bivariate;
# - complete bivariate.
#
# The corresponding local identification method is then applied.
# -------------------------------------------------------------------------
classify_local_pattern <- function(data, n, dyad, variable, member,
alpha = 0.05) {
other_member <- if (member == 1L) 2L else 1L
second_variable <- if (variable == "V1") "V2" else "V1"
# Main variable
main_first <- simulation_chain(
data, n, dyad, variable, member
)
main_second <- simulation_chain(
data, n, dyad, variable, other_member
)
# Second variable, using the same member order
second_first <- simulation_chain(
data, n, dyad, second_variable, member
)
second_second <- simulation_chain(
data, n, dyad, second_variable, other_member
)
# Empirical 16 x 2 transition-count matrix
empirical <- dyadicMarkov::countEmpBivariate(
chainFM_V1 = main_first,
chainSM_V1 = main_second,
chainFM_V2 = second_first,
chainSM_V2 = second_second,
states = 2L
)
# Global bivariate case
global <- dyadicMarkov::bivariateCase(
empirical,
alpha = alpha
)
# Proceed to the appropriate local identification
if (global$case == "trivial") {
return("trivial")
}
if (global$case == "univariate") {
local <- dyadicMarkov::univariatePattern(
main_first,
main_second,
states = 2L,
alpha = alpha
)
return(univariate_code(local))
}
if (global$case == "partial") {
local <- dyadicMarkov::partialPattern(empirical)
return(selected_matrix_code(local))
}
if (global$case == "complete") {
local <- dyadicMarkov::completePattern(empirical)
return(selected_matrix_code(local))
}
stop("Unknown bivariate case: ", global$case)
}
# -------------------------------------------------------------------------
# Helper 4: run the classification for all 1,000 dyads at one sequence
# length.
#
# Each dyad contributes four member-variable combinations:
#
# V1, first member
# V1, second member
# V2, first member
# V2, second member
#
# The simulation label is read directly from caseSimulated rather than
# entered manually.
# -------------------------------------------------------------------------
classify_one_length <- function(data, n, alpha = 0.05) {
dyads <- sort(unique(as.integer(data$dyad)))
orientations <- data.frame(
variable = c("V1", "V1", "V2", "V2"),
member = c(1L, 2L, 1L, 2L),
stringsAsFactors = FALSE
)
out <- vector(
"list",
length(dyads) * nrow(orientations)
)
k <- 1L
for (i in seq_along(dyads)) {
d <- dyads[i]
# Optional progress information when the full analysis is run locally
if (i %% 100L == 0L) {
message(
"Length ", n, ": completed ",
i, " of ", length(dyads), " dyads"
)
}
for (j in seq_len(nrow(orientations))) {
variable <- orientations$variable[j]
member <- orientations$member[j]
# Pattern used to simulate this member-variable combination
variable_code <- unname(simulation_variables[[variable]])
keep <- data$dyad == as.character(d) &
data$variable == variable_code &
data$members == as.character(member)
simulated <- unique(data$caseSimulated[keep])
if (length(simulated) != 1L) {
stop(
"Expected exactly one simulation label for dyad ",
d, ", variable ", variable,
", member ", member, "."
)
}
# Pattern identified by the current dyadicMarkov workflow
identified <- classify_local_pattern(
data = data,
n = n,
dyad = d,
variable = variable,
member = member,
alpha = alpha
)
out[[k]] <- data.frame(
length = n,
dyad = d,
variable = variable,
member = member,
simulated = simulated,
identified = identified,
stringsAsFactors = FALSE
)
k <- k + 1L
}
}
do.call(rbind, out)
}
# -------------------------------------------------------------------------
# Run all five sequence lengths.
#
# 5 lengths x 1,000 dyads x 4 member-variable combinations
# = 20,000 classifications.
# -------------------------------------------------------------------------
simulation_data <- list(
`30` = data_complete_30,
`60` = data_complete_60,
`90` = data_complete_90,
`180` = data_complete_180,
`720` = data_complete_720
)
classification_results <- do.call(
rbind,
lapply(names(simulation_data), function(n) {
classify_one_length(
data = simulation_data[[n]],
n = as.integer(n),
alpha = 0.05
)
})
)
# -------------------------------------------------------------------------
# Calculate sensitivity and specificity for each pattern and sequence length.
#
# Sensitivity = TP / (TP + FN)
# Specificity = TN / (TN + FP)
#
# The simulated pattern is treated as the reference classification.
# -------------------------------------------------------------------------
patterns <- c("A1", "B2", "D4", "E4")
lengths <- c(30L, 60L, 90L, 180L, 720L)
sensitivity_results <- data.frame(
Length = lengths
)
for (pattern in patterns) {
sensitivity <- numeric(length(lengths))
specificity <- numeric(length(lengths))
for (i in seq_along(lengths)) {
current <- classification_results[
classification_results$length == lengths[i],
]
truth <- current$simulated == pattern
identified <- current$identified == pattern
TP <- sum(truth & identified)
FN <- sum(truth & !identified)
TN <- sum(!truth & !identified)
FP <- sum(!truth & identified)
sensitivity[i] <- TP / (TP + FN)
specificity[i] <- TN / (TN + FP)
}
sensitivity_results[[paste(pattern, "Se")]] <-
round(sensitivity, 3)
sensitivity_results[[paste(pattern, "Sp")]] <-
round(specificity, 3)
}
sensitivity_resultsThe complete calculation above is not run during an ordinary vignette
build. The sensitivity_results values obtained from that
reproduction, rounded to three decimals, are stored below so that the
vignette can display the results without rerunning all 20,000
classifications. They reproduce the sensitivity and specificity results
reported in Table 5 of Böllenrücher et al. (in
press), with minor numerical differences in a few specificity
values.
# Results obtained from the full sensitivity-analysis reproduction above.
# Values are stored here so that the vignette does not rerun the 20,000
# classifications every time it is built.
sensitivity_results <- data.frame(
Length = c(30L, 60L, 90L, 180L, 720L),
`A1 Se` = c(0.125, 0.836, 0.936, 0.941, 0.946),
`A1 Sp` = c(0.992, 0.942, 0.951, 0.997, 1.000),
`B2 Se` = c(0.114, 0.680, 0.801, 0.835, 0.850),
`B2 Sp` = c(0.994, 0.940, 0.935, 0.991, 1.000),
`D4 Se` = c(0.000, 0.119, 0.426, 0.876, 0.951),
`D4 Sp` = c(0.999, 0.972, 0.955, 0.952, 0.960),
`E4 Se` = c(0.015, 0.358, 0.621, 0.765, 0.808),
`E4 Sp` = c(1.000, 0.997, 0.997, 0.999, 0.998),
check.names = FALSE
)
knitr::kable(
sensitivity_results,
digits = 3,
align = c("r", rep("r", 8)),
caption = paste(
"Sensitivity (Se) and specificity (Sp) for patterns A1, B2, D4,",
"and E4 across the five sequence lengths."
)
)| Length | A1 Se | A1 Sp | B2 Se | B2 Sp | D4 Se | D4 Sp | E4 Se | E4 Sp |
|---|---|---|---|---|---|---|---|---|
| 30 | 0.125 | 0.992 | 0.114 | 0.994 | 0.000 | 0.999 | 0.015 | 1.000 |
| 60 | 0.836 | 0.942 | 0.680 | 0.940 | 0.119 | 0.972 | 0.358 | 0.997 |
| 90 | 0.936 | 0.951 | 0.801 | 0.935 | 0.426 | 0.955 | 0.621 | 0.997 |
| 180 | 0.941 | 0.997 | 0.835 | 0.991 | 0.876 | 0.952 | 0.765 | 0.999 |
| 720 | 0.946 | 1.000 | 0.850 | 1.000 | 0.951 | 0.960 | 0.808 | 0.998 |
Sensitivity generally increases with sequence length, whereas specificity remains high across all five lengths. With 30 measurement points, sensitivity is low for all four patterns and is particularly low for the more complex bivariate patterns D4 and E4. Their identification improves progressively as more transition information becomes available. At 90 measurement points, sensitivity is already high for A1 and B2 and has also increased substantially for D4 and E4; further improvements are observed at 180 and 720 measurement points.
These results are consistent with Böllenrücher et al. (in press), where shorter sequences are more frequently associated with trivial classifications and the identification of more complex patterns requires longer sequences. From 90 measurement points onward, trivial classifications decrease markedly, leading Böllenrücher et al. (in press) to recommend approximately 90 measurement points as a practical minimum for the bivariate method.
This recommendation is methodological guidance from the sensitivity
study, rather than a requirement imposed by
dyadicMarkov.