---
title: "Bivariate dyadic workflow"
bibliography: references.bib
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Bivariate dyadic workflow}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
```

## Overview

This vignette shows the bivariate workflow implemented in `dyadicMarkov`. In the bivariate setting, two categorical variables are observed repeatedly for the two members of a dyad. The bivariate method follows the global-and-local procedure described in @bollen2026.

The bivariate method uses matrix codes to identify the local dependence patterns. Partial bivariate patterns are denoted B1--B3, while complete bivariate patterns are denoted C, D1--D4, and E1--E4. When the global step identifies a univariate case, the A-family codes described in the univariate workflow apply. The pattern nomenclature is summarized in Table 2 of @bollen2026.

The current bivariate functions support `states = 2`. With two binary variables observed for two members, the previous state is described by four binary components. The empirical bivariate count matrix therefore has 16 rows and 2 columns.

Although the data are synthetic, the four columns can be read like real repeated observations from a dyad. For example, `V1` could represent one coded behavior or response, and `V2` could represent a second coded behavior or response observed at the same measurement occasions. The columns `FM_V1` and `SM_V1` then describe the two members on the main variable, while `FM_V2` and `SM_V2` describe the same two members on the second variable.

The 16 rows arise because the previous state combines four binary lagged components: first member on `V1`, second member on `V1`, first member on `V2`, and second member on `V2`. With two possible states for each component, this gives $2^4 = 16$ previous-state combinations. The 2 columns represent the possible next states of the first member on the current main variable.

## Data

The example data set `dyadic_bivariate_example` contains two categorical variables for the first member and the second member of a dyad. Each row corresponds to one measurement occasion.

```{r bivariate-data}
utils::data("dyadic_bivariate_example", package = "dyadicMarkov")

head(dyadic_bivariate_example)
dim(dyadic_bivariate_example)
```

The four chains are first-member and second-member sequences for the main variable (`V1`) and the second variable (`V2`).

```{r bivariate-states}
table(dyadic_bivariate_example$FM_V1)
table(dyadic_bivariate_example$SM_V1)
```

## Empirical bivariate transition counts

The first step is to construct the empirical bivariate transition count matrix with `countEmpBivariate()`. Rows represent the 16 possible previous-state combinations of the two dyadic variables. Columns represent the next state of the first member on the main variable.

```{r bivariate-counts}
emp_bi <- dyadicMarkov::countEmpBivariate(
  chainFM_V1 = dyadic_bivariate_example$FM_V1,
  chainSM_V1 = dyadic_bivariate_example$SM_V1,
  chainFM_V2 = dyadic_bivariate_example$FM_V2,
  chainSM_V2 = dyadic_bivariate_example$SM_V2,
  states = 2L
)

emp_bi
class(emp_bi)
dim(emp_bi)
```

## Global bivariate case

The function `bivariateCase()` performs the global step of the bivariate method. The global approach compares nested models within the likelihood-ratio test (LRT) framework. `bivariateCase()` implements the two chi-squared tests for the A1 and B1 comparisons used to identify the dependence case. `dyadicMarkov` evaluates these tests using Pearson's chi-squared statistic, $X^2 = \sum (O-E)^2/E$, to classify the analyzed sequence as a trivial, univariate, partial bivariate, or complete bivariate case.

```{r bivariate-case}
case_bi <- dyadicMarkov::bivariateCase(emp_bi, alpha = 0.05)

case_bi
case_bi$case
summary(case_bi)
```

This example is identified as a complete bivariate case. The appropriate local step is therefore to compare complete bivariate candidate patterns.

## Local pattern identification for a complete case

For a complete bivariate case, `completePattern()` first computes the G-squared deviance, $G^2 = 2\sum O\log(O/E)$, for each complete bivariate candidate structure. It then calculates $AIC = G^2 + 2k$ and returns the candidate with the smallest AIC.

```{r complete-pattern}
complete_bi <- dyadicMarkov::completePattern(emp_bi)

complete_bi
complete_bi$pattern
complete_bi$aic
```

In this example, the selected complete bivariate pattern is `D2`, labelled by the package as actor only on the main, actor-partner on the second.

## Analyzing each variable and member in turn

Each bivariate analysis is defined from a specific perspective. The sequence supplied as the first member is the sequence being analyzed, while the second member provides the partner sequence. Likewise, one variable is treated as the main variable and the other as the second variable. Swapping the two members therefore changes the member perspective of the analysis, while swapping the two variables changes which variable is treated as the main variable. To describe the dyad more completely, the workflow can be repeated for each combination of analyzed member and main variable. These are distinct analyses and may therefore lead to different global cases and local interaction patterns.

For compactness, the following vignette-local helper applies the exported functions in sequence; `analyze_bivariate()` is not part of the package API.

```{r bivariate-sequences-in-turn}
analyze_bivariate <- function(label, fm_v1, sm_v1, fm_v2, sm_v2) {
  emp <- dyadicMarkov::countEmpBivariate(
    chainFM_V1 = fm_v1,
    chainSM_V1 = sm_v1,
    chainFM_V2 = fm_v2,
    chainSM_V2 = sm_v2,
    states = 2L
  )

  case <- dyadicMarkov::bivariateCase(emp, alpha = 0.05)

  cat("\n", label, "\n", sep = "")
  print(case)

  if (identical(case$case, "complete")) {
    print(dyadicMarkov::completePattern(emp))
  }

  if (identical(case$case, "partial")) {
    print(dyadicMarkov::partialPattern(emp))
  }

  if (identical(case$case, "univariate")) {
    print(dyadicMarkov::univariatePattern(fm_v1, sm_v1, states = 2L, alpha = 0.05))
  }
}

d <- dyadic_bivariate_example

analyze_bivariate(
  "FM_V1 as analyzed sequence, V1 as main variable",
  d$FM_V1, d$SM_V1, d$FM_V2, d$SM_V2
)

analyze_bivariate(
  "SM_V1 as analyzed sequence, V1 as main variable",
  d$SM_V1, d$FM_V1, d$SM_V2, d$FM_V2
)

analyze_bivariate(
  "FM_V2 as analyzed sequence, V2 as main variable",
  d$FM_V2, d$SM_V2, d$FM_V1, d$SM_V1
)

analyze_bivariate(
  "SM_V2 as analyzed sequence, V2 as main variable",
  d$SM_V2, d$FM_V2, d$SM_V1, d$FM_V1
)
```

For this example, analyzing the four sequences in turn illustrates three branches of the procedure: complete bivariate cases when `FM_V1` and `SM_V1` are analyzed, a partial bivariate case for `FM_V2`, and a univariate case for `SM_V2`. The partial branch is therefore demonstrated with the same bivariate example data rather than with an artificial seeded example.

## Reading the global and local steps together

The global and local steps should be read together while keeping their statistics distinct. The global step compares nested models within the LRT framework through two chi-squared tests evaluated using Pearson's chi-squared statistic, $X^2$. The local partial and complete steps use the G-squared deviance, $G^2$, to calculate candidate AIC values. If `bivariateCase()` returns `trivial`, no subsequent local pattern is selected. If it returns `univariate`, the bivariate workflow returns to `univariatePattern()` using the first- and second-member sequences of the current main variable. If it returns `partial`, the local step is `partialPattern()`. If it returns `complete`, the local step is `completePattern()`.

## References