This is a step-by-step guide to illustrate the contents of the paper “Quantifying the therapeutic alliance: A methodological framework for conducting observational moment-by-moment research”.

Setup

Custom functions

The analyses in this tutorial rely on four custom R functions (permanent_group, sequence_descriptives, sequence_heatmap, and sequence_trajectory), all defined in dyadObs_functions.r. The chunk below downloads this file to your R working directory and loads it. To find or edit the file afterwards, run getwd() in the R console — it will be saved there as dyadObs_functions.r.

#download.file(
#  "https://raw.githubusercontent.com/dfernandezregueras/dyadObs/main/dyadObs_functions.r",
#  destfile = "dyadObs_functions.r"
#)
source("dyadObs_functions.r")

Data files

This tutorial requires four data files. Download each one using the links below and save all of them to the folder shown by getwd() in R.

1. Preprocessing (observational data)

1.1 Load the observational data and merge the different files

Sometimes, data from different observations are stored in separate files. To prepare the data, we will merge them into a single file. The Observer XT provides .xlsx data files that contain several variables. For the data analysis, only a few of them are relevant:

  • Behavior: Contains the names of the coded events in a longitudinal format.
  • Observation: Contains the name assigned to each observation session. When coding, researchers must ensure that each session has a unique name. A common practice is to include an ID for the case and the session number (e.g., case1_session1). This variable represents repeated measures in a long format, as it includes all coded events within a session.
  • Subject: Identifies the individual who emits the behavior (e.g., “Therapist” or “Client”).

The researcher may choose to retain the remaining variables. However, for simplicity, we will remove them since they are not needed for our analyses. If the researcher chooses to retain the additional variables and the data frames contain different variables (as in our case), the merge function should be used.

library(readxl)

# Load the data files

data_1 <- read_excel("data_1.xlsx")[, c("Observation", "Subject", "Behavior")]
data_2 <- read_excel("data_2.xlsx")[, c("Observation", "Subject", "Behavior")]
data_3 <- read_excel("dyad 31.xlsx")[, c("Observation", "Subject", "Behavior")]
data_4 <- read_excel("dyad 32.xlsx")[, c("Observation", "Subject", "Behavior")]



# Merge the data (only three variables)

data <- rbind(data_1, data_2, data_3, data_4)


library(writexl)
write_xlsx(data, "Previous data.xlsx")

# Maintain the remaining variables

data <- Reduce(function(x, y) merge(x, y, by = c("Observation", "Subject", "Behavior"), all = T),
                      list(data_1, data_2, data_3, data_4))

1.2 Assign the sessions to the remaining variables: Phase of therapy, ID, Group

library(readxl)

data <- read_excel("Previous data.xlsx")



# ASSIGN ID & Treatment phase

# Visualize the list of Observation sessions

unique_sessions <- unique(data$Observation)

unique_sessions <- as.data.frame(unique_sessions)


# Assign ID and Treatment phases based on the list of Observation sessions

session_lookup <- data.frame(
  Observation = unique(data$Observation),
  ID = rep(1:32, times = c(5, 5, 5, 5, 5, 5, 3, 5, 5, 5, 4, 5, 5, 5, 5, 5, 5, 5, 3, 5, 5, 5, 5, 5, 4, 5, 5, 5, 5, 5, 5, 5)),
  Treatment_Phase = c("Middle of intervention", "End of intervention", "Assessment", "Treatment planning", "Beginning of intervention",
                      "End of intervention","Assessment","Treatment planning","Beginning of intervention","Middle of intervention",
                      "Middle of intervention","End of intervention","Assessment","Treatment planning","Beginning of intervention",
                      "Beginning of intervention","Middle of intervention", "Assessment", "End of intervention", "Treatment planning",
                      "Middle of intervention","End of intervention","Assessment","Treatment planning","Beginning of intervention",
                      "Beginning of intervention","Middle of intervention","End of intervention","Assessment","Treatment planning",
                      "Assessment","Treatment planning","Beginning of intervention",
                      "Middle of intervention","End of intervention","Assessment","Treatment planning","Beginning of intervention",
                      "Middle of intervention","End of intervention","Assessment","Treatment planning","Beginning of intervention",
                      "Middle of intervention","End of intervention","Assessment","Treatment planning","Beginning of intervention",
                      "Assessment","Treatment planning","Beginning of intervention", "End of intervention",
                      "Middle of intervention","End of intervention","Assessment","Treatment planning","Beginning of intervention",
                      "Beginning of intervention","Middle of intervention","End of intervention","Assessment","Treatment planning",
                      "Middle of intervention","Assessment", "End of intervention", "Treatment planning","Beginning of intervention",
                      "Beginning of intervention","Middle of intervention","End of intervention","Assessment","Treatment planning",
                      "Middle of intervention","End of intervention","Assessment","Treatment planning","Beginning of intervention",
                      "Middle of intervention","End of intervention","Assessment","Beginning of intervention","Treatment planning",
                      "Beginning of intervention","Middle of intervention","End of intervention","Assessment","Treatment planning",
                      "Assessment","Treatment planning","Beginning of intervention",
                      "End of intervention","Assessment","Treatment planning","Beginning of intervention","Middle of intervention",
                      "Beginning of intervention","Middle of intervention","Assessment","End of intervention","Treatment planning",
                      "Beginning of intervention","Middle of intervention","Assessment","End of intervention","Treatment planning",
                      "Middle of intervention","End of intervention","Assessment","Treatment planning","Beginning of intervention",
                      "Middle of intervention","End of intervention","Assessment","Treatment planning","Beginning of intervention",
                      "Assessment","Treatment planning","Beginning of intervention","End of intervention",
                      "Middle of intervention","End of intervention","Assessment","Treatment planning","Beginning of intervention",
                      "End of intervention","Assessment","Treatment planning","Beginning of intervention","Middle of intervention",
                      "Beginning of intervention","Middle of intervention","End of intervention","Assessment","Treatment planning",
                      "Middle of intervention","End of intervention","Assessment","Treatment planning","Beginning of intervention",
                      "End of intervention","Assessment","Treatment planning","Beginning of intervention","Middle of intervention",
                      "Middle of intervention","End of intervention","Assessment","Treatment planning","Beginning of intervention",
                      "Middle of intervention","End of intervention","Assessment","Treatment planning","Beginning of intervention"
                      )
)


# Merge the assigned IDs and Treatment phases with the rest of data

library(dplyr)

data_merged <- left_join(data, session_lookup, by = "Observation")


# Check phases' names for potential typos

#unique(data_merged$Treatment_Phase)


##################################################################################

# ASSIGN GROUPS

# Visualize the list of IDs

unique_ID <- unique(data_merged$ID)

unique_ID <- as.data.frame(unique_ID)

# Assign Groups based on the list of IDs

ID_lookup <- data.frame(
  ID = unique(data_merged$ID),
  Group = c("Poor", "Poor", "Good", "Poor", "Poor", "Good", "Poor", "Good", "Good", "Poor", "Good", "Good", "Poor", "Good", "Good", "Poor", "Poor", "Good", "Poor", "Good", "Average", "Average", "Average", "Average", "Average", "Average", "Average", "Average", "Average", "Average", "Average", "Average")
)


# Merge the assigned Grups with the rest of data

data_complete <- left_join(data_merged, ID_lookup, by = "ID")

# Check Groups' names for potential typos

#unique(data_complete$Group)

head(data_complete)
library(writexl)
write_xlsx(data_complete, "Preprocessed data - phase names.xlsx")

1.3 Permanent grouping of behaviors

Some behaviors, either for theoretical reasons or due to their low frequency, need to be grouped. In future sections, we will present a more flexible approach for grouping behaviors for specific analyses without the need to permanently modify the dataset. However, some researchers may prefer to group behaviors permanently beforehand if they are certain that these behaviors should remain unified under the same label. In such cases, researchers may find the permanent_group function useful.

  • permanent_group allows researchers to group behaviors under a single label before conducting analyses. Its simple usage requires researchers to input the names of the behaviors to be grouped, separated by a + symbol, as well as the name of the variable that store the behaviors and the name of the new grouped behavior created.
library(readxl)
data <- read_excel("Preprocessed data - phase names.xlsx")

As an example, the client behaviors Partial agreement and Disagreement, which are theoretically related, will be unified into the new behavior Group_1.

data <- permanent_group(data, "Behavior", "Partial agreement C + Disagreement C", "Group_1")

After preprocessing the data and obtaining the final dataset, it is always good practice to store it safely to avoid losing information and repeating the entire process.

library(writexl)
write_xlsx(data, "Preprocessed data.xlsx")

2. Descriptive statistics (observational data)

Descriptive data will be useful for future decisions regarding the analyses. To provide a comprehensive description of the data, the function sequence_descriptives is able to provide descriptive statistics for both individual behaviors and sequences of an indefinite number of behaviors. The contents of this section will, therefore, apply to both individual behaviors and sequences.

2.1 Individual behaviors and sequences: descriptive statistics, normality tests and density distribution plots

The sequence_descriptives function is designed to analyze the relative frequencies of behavior sequences across multiple phases of treatment. It provides descriptive statistics, tests for normality (Shapiro-Wilk test), and generates histograms with density distributions to visualize the distribution of behaviors/sequences. The function is highly flexible, allowing researchers to group behaviors using the + symbol without permanently altering the dataset.

Key Features

  1. Flexible Behavior Grouping:

    • Researchers can group behaviors using the + symbol (e.g., "BehaviorA + BehaviorB") to treat them as the same during analysis.
    • This grouping is applied dynamically and does not modify the original dataset.

  2. Descriptive Statistics:

    • Computes the following statistics for each treatment phase:
      • Sample size: Number of observation sessions per phase.
      • Mean: Mean relative frequency of behaviors/sequences.
      • Median: Median relative frequency of behaviors/sequences.
      • Standard deviation: Variability of relative frequencies.
      • Percentile 5: 5th percentile of relative frequencies.
      • Percentile 95: 95th percentile of relative frequencies.
      • Skewness: Measure of distribution asymmetry.
      • Standard error of skewness: Error associated with skewness.
      • Kurtosis: Measure of distribution tailedness.
    • Includes the Shapiro-Wilk test for normality, with a p-value and a variable indicating whether the distribution is Normal or Not normal.

  3. Visualizations:

    • Generates boxplots for each treatment phase.
    • Allows combining plots for all phases into a single display and group comparisons based on a variable (e.g. Intervention A vs Intervention B). Boxplots comparing groups could be displayed in the same plot or separated in different plots.

  4. Group Analysis:

    • Provides options to filter and compare specific groups or phases.

  5. Relative or Absolute Frequencies:

    • Allows absolute or relative frequencies analyses of behaviors/sequences per observation session.

  6. Requested Libraries Handling:

    • Checks for required packages (dplyr, ggplot2, gridExtra, e1071 and stats [included in the in base R]) and provides a warning if they are not installed.

Function Arguments

Argument Description
data The dataset containing the observations. It should be a data.frame where each row represents a subject’s observation within a phase of the treatment.
treatment_phase_var The name of the variable that defines the treatment phase in the dataset (e.g., "Treatment_Phase"). This variable is used to group observations by phase.
phase_order A vector specifying the order of the treatment phases to be analyzed. Example: c("Phase1", "Phase2", "Phase3"). If NULL, phases are ordered as they appear in the data.
behavior_var The name of the variable containing the observed behaviors for each subject (e.g., "Behavior"). This variable should include behavior sequences to analyze.
sequence A vector or list of strings specifying the behavior sequence to analyze. Use + to group behaviors (e.g., c("BehaviorA + BehaviorB", "BehaviorC").
combine_plots A logical value (TRUE or FALSE). If TRUE, boxes for all groups are combined into a single display. If FALSE, plots are displayed separately for each group.
group_var The name of the grouping variable used to divide the data into different groups (e.g., "Group"). If NULL, the data is analyzed as a whole.
selected_groups A vector of group names to filter and analyze. Example: c("Group1", "Group2"). If NULL, all groups in group_var are analyzed.
use_relative_frequency A logical value (TRUE or FALSE). If TRUE, relative frequencies will be computed. If FALSE, absolute frequencies will be computed. 
library(readxl)
data <- read_excel("Preprocessed data - phase names.xlsx")

Output and usage example

results <- sequence_descriptives(
  data = data,
  treatment_phase_var = "Treatment_Phase",
  phase_order = c("Assessment", "Treatment planning", "Beginning of intervention",
                  "Middle of intervention", "End of intervention"),
  behavior_var = "Behavior",
  sequence = c("Empathy T", "Total agreement C"),
  combine_plots = T,  # FALSE for separate plots per group
  group_var = "Group",
  selected_groups = c("Good", "Poor"),
  use_relative_frequency = T
)

library(DT)

stats_df <- results$statistics


numeric_cols <- names(stats_df)[sapply(stats_df, is.numeric)]

datatable(stats_df,
          options = list(
            pageLength = 5,
            scrollX = TRUE,
            dom = 'Bfrtip'
          )) %>%
  formatRound(columns = setdiff(numeric_cols, "Sample_Size"), digits = 3) %>%
  formatRound(columns = "Sample_Size", digits = 0)

3. Data analyses: answering the proposed research questions

3.1 What behaviors/sequences are more frequent?

The sequence_heatmap function is designed to analyze and visualize individual or sequences of behaviors. It generates a heatmap that displays the absolute or relative frequency of specific behaviors or sequences of behaviors across different moments (e.g., treatment phases). The function is highly customizable, allowing users to define behaviors/sequences and treatment phases, as well as specify the type of data desired (absolute or relative frequency). The researcher will also be able to flexibly group behaviors that they want to be treated as the same using a +, without permanently changing the data frame. It also allows for group comparison.

Key Features

  1. Behavior and Sequence Analysis:
    • Handles grouped behaviors flexibly (e.g., "BehaviorA + BehaviorB"), treating them as one single behavior.
  2. Customizable Inputs:
    • Users can specify the columns for behaviors, treatment phases, and sessions in the dataset.
    • Allows customization of treatment phase levels and the order of behaviors/sequences in the heatmap.
  3. Frequency Calculation:
    • Computes absolute frequencies (counts) or relative frequencies (proportions) of behaviors or sequences.
    • Calculates mean frequencies per phase by normalizing the sum of frequencies by the total number of sessions in each phase.
  4. Heatmap Visualization:
    • Generates a heatmap using ggplot2, with behaviors/sequences on the y-axis and treatment phases on the x-axis.
    • The fill color of the heatmap represents the absolute or relative frequency of behaviors/sequences.
  5. Custom Ordering:
    • Allows users to order behaviors/sequences in the heatmap either:
      • Personalized: In the same order as provided in custom_sequences.
      • By Frequency: Based on the total frequency of behaviors/sequences (highest to lowest).
  6. Group Separation:
    • Supports separation of data by groups (e.g., "Good", "Poor") if a group_var is provided.
    • Generates separate heatmaps for each group, ensuring that the analysis is tailored to each subset of the data.
  7. Requested Libraries Handling:
    • Checks for and loads required packages (dplyr, ggplot2, tidyr, purrr, stringr) and installs them if they are not installed.

Function Arguments

Argument Description
data The input dataset containing behavioral data. This should be a data frame or with columns for behaviors, treatment phases, session identifiers and, optionally, groups.
behavior_col The column name in data that contains behaviors (default: "Behavior").
phase_col The column name in data that contains treatment phases (default: "Treatment_Phase"). This column should categorize observations into distinct phases.
session_col The column name in data that contains session identifiers (default: "Observation"). This column should uniquely identify each session or observation.
custom_sequences A vector of behaviors or sequences to analyze. Provide a list of behaviors and/or sequences using the c() command. Example: list(c("BehaviorA + BehaviorB", "BehaviorC"), c("BehaviorD", "BehaviorE"), c("Behavior F + Behavior G")).
phase_levels A vector specifying the order of treatment phases (default: NULL). Introduce a c() with the ordered names of the phases. Example: c("Phase1", "Phase2", "Phase3"). If NULL, phases are ordered as they appear in the data.
freq_type The type of frequency to compute: "absolute" (default) or "relative". Use "absolute" for raw counts or "relative" for proportions relative to total behaviors/sequences in each session. It computes the mean value across all sessions in each treatment phase.
custom_order The ordering of behaviors/sequences in the heatmap: "personalized" (default) or "frequency". Use "personalized" to maintain the order in which behaviors/sequences were introduced in custom_sequences. Use "frequency" to order behaviors/sequences from higher to lower total frequency across all phases.
group_var The column name in data that contains group identifiers (default: NULL). If provided, the function will generate separate heatmaps for each group. Example: "Group".
selected_groups A vector of groups to include in the analysis (default: NULL). If provided, only the specified groups will be analyzed. Example: c("Good", "Poor").

Output: example usage

The function returns a heatmap or various heatmaps depending if the group argument includes a group variable. The following callings show different possible outputs.

2 groups compared - Relative frequency - Personalized order

sequence_heatmap(
  data = data,
  behavior_col = "Behavior",
  phase_col = "Treatment_Phase",
  session_col = "Observation",
  custom_sequences = list(
    c("Task with justification T", "Total agreement C"),
    c("Task with justification T", "Disagreement C + Partial agreement C"),
    c("Task with justification T", "Partial agreement C"),
    c("Task without justification T", "Total agreement C"),
    c("Feedback task T", "Total agreement C"),
    c("Empathy T", "Total agreement C"),
    c("Present success T", "Total agreement C")
  ),
  phase_levels = c("Assessment", "Treatment planning", "Beginning of intervention", "Middle of intervention", "End of intervention"),
  freq_type = "relative",
  custom_order = "personalized",
  group_var = "Group",
  selected_groups = c("Good", "Poor")
)

No group comparison - Absolute frequency - Order by total frequency

sequence_heatmap(
  data = data,
  behavior_col = "Behavior",
  phase_col = "Treatment_Phase",
  session_col = "Observation",
  custom_sequences = list(
    c("Total agreement T"),
    c("Empathy T"),
    c("Partial agreement C"),
    c("Task without justification T", "Total agreement C"),
    c("Partial agreement C + Disagreement C")),
  phase_levels = c("Assessment", "Treatment planning", "Beginning of intervention", "Middle of intervention", "End of intervention"),
  freq_type = "absolute",
  custom_order = "frequency",
  group_var = NULL,
  selected_groups = c("Good", "Poor") # Only works if group_var is not NULL
)

3.2 How do behaviors/sequences evolve over time?

The sequence_trajectory function models how specific behavior sequences evolve across treatment phases using a standardized metric to allow comparison between different behaviors/sequences and different groups. It enables researchers to explore patterns of change (e.g., increasing, decreasing, stable) in the relative frequency of sequences, optional comparisons between groups. The function supports multiple polynomial shapes (linear, quadratic, cubic) to fit and visualize behavioral trends.

Key Features

  1. Trajectory Modeling:
    • Computes the relative frequency in a standardized metric of specific behaviors or behavior sequences across treatment phases.
    • Aggregates and averages frequencies to identify trends across sessions.
  2. Group-Based Analysis:
    • If a grouping variable is provided, trajectories are computed separately for each group (e.g., "Good" vs "Poor"), enabling comparative analysis.
  3. Behavior Grouping with Flexibility:
    • Supports combined behavior analysis using the + operator (e.g., "BehaviorA + BehaviorB") for evaluating equivalent/similar behaviors.
  4. Clean Visualization:
    • Plots behavioral trajectories using ggplot2 with clear phase separation and optional faceting by group.
  5. Phase Control:
    • Users can control the ordering of treatment phases.
  6. Flexible Confidence Intervals:
    • Allows users to choose between t-based (default) and bootstrap confidence intervals in the plots.
    • This is especially useful when working with small samples or data that may depart from normality.

Argument Descriptions

Argument Description
data The input data frame containing behavioral observation data.
sequence A character vector (or list of vectors) representing the target behavior and/or sequence(s) to analyze. Each sequence is matched as an exact temporal pattern.
behavior_var The column name containing behavior codes (default: "Behavior").
phase_var The column name indicating treatment phases (default: "Treatment_Phase").
phase_order Optional. A vector specifying the desired order of phases. If NULL, the function uses the factor levels in the data.
observation_var Column identifying individual observations/sessions (default: "Observation").
id_var Column identifying unique subjects or participants (default: "ID").
group_variable Optional. The column name for a grouping variable (e.g., "Group"). Required if comparing across groups.
groups_compared Optional. A vector of group labels to include in the analysis (e.g., c("Treatment", "Control")). Required if group_variable is specified.
shapes A character vector indicating the polynomial shapes to fit: "linear", "quadratic", and/or "cubic" (default includes all three).
ci_method The method used to compute confidence intervals in the plots. Options are "t" (default, using the t-distribution) and "bootstrap" (using bootstrap percentile confidence intervals).
boot_seed Optional integer used to make bootstrap confidence intervals reproducible. Only relevant when ci_method = "bootstrap". If NULL (default), the random-number stream is not seeded and error bars will vary slightly between runs.

Output: example usage

The function returns a list with plots and model information:

  • For single-group analysis:
    • A list of results for each sequence.
    • A combined vertical plot (all_sequences_plot) showing the relative frequency trajectories across phases.
  • For group comparisons:
    • Results broken down by group and sequence.
    • A comparison_plot showing trajectory plots side-by-side across groups for all sequences.

Each plot includes: - Individual lines for each subject. - Group average trajectory with confidence intervals. - Best-fit polynomial model (in red) with AIC info.

Two single behaviors compared

The output shown answers the research question How do behaviors/sequences evolve over time?, so the group_variable argument is NULL.

# Two single behaviors compared

result <- sequence_trajectory(
  data = data,
  sequence = list(c("Task with justification T"), c("Empathy T")),
  behavior_var = "Behavior",
  phase_var = "Treatment_Phase",
  phase_order = c("Assessment", "Treatment planning", "Beginning of intervention", "Middle of intervention", "End of intervention"),
  observation_var = "Observation",
  id_var = "ID",
  group_variable = NULL,
  groups_compared = NULL,
  shapes = c("linear", "quadratic", "cubic")
)

result[["all_sequences_plot"]]

3.3 Do changes in behaviors/sequences differ across groups?

The sequence_trajectory allows for comparisons between groups using the argument group_variable, which needs to be filled with the name of the variable that contains the groups’ names. In the argument groups_compared, we need to indicate which groups from the group_variable we want to compare.

Output

One sequence - Two groups compared

result <- sequence_trajectory(
  data = data,
  sequence = list(c("Task with justification T", "Total agreement C")),
  behavior_var = "Behavior",
  phase_var = "Treatment_Phase",
  phase_order = c("Assessment", "Treatment planning", "Beginning of intervention", "Middle of intervention", "End of intervention"),
  observation_var = "Observation",
  id_var = "ID",
  group_variable = "Group",
  groups_compared = c("Good", "Poor"),
  shapes = c("linear", "quadratic", "cubic")
)

result[["comparison_plot"]]

Two sequences - Two groups compared

result <- sequence_trajectory(
  data = data,
  sequence = list(c("Summarizing T + Clarifying T + Finishing sentence T", "Total agreement C"), c("Empathy T", "Total agreement C")),
  behavior_var = "Behavior",
  phase_var = "Treatment_Phase",
  phase_order = c("Assessment", "Treatment planning", "Beginning of intervention", "Middle of intervention", "End of intervention"),
  observation_var = "Observation",
  id_var = "ID",
  group_variable = "Group",
  groups_compared = c("Good", "Poor"),
  shapes = c("linear", "quadratic", "cubic")
)

result[["comparison_plot"]]