Identify groups/associations among longitudinal binary variables

Question

This question is an extension of a prior question. I have a longitudinal dataset with several binary variables, along with an id variable and a time variable. In my research project, these binary variables are observed behaviors during a videotaped task (0 = absent, 1 = present) and I am interested to examine how the behaviors correlate over time and which (if any) variables cluster/group or vary together across time. Conceptually, these observed behaviors are similar to items on a questionnaire and I want to see what kind of "subscales" appear in my data.

My prior question revealed that exploratory factor analysis (EFA) and principal components analysis (PCA) are likely not appropriate. EFA is controversial with binary data. PCA is conceptually inappropriate for my research question since I want to find groups/associations among variables, not reduce them.

Prior commenters mentioned multiple factor analysis (MFA) and multiple correspondence analysis (MCA). If I understand correctly, MFA requires that I already have identified groups among variables, which I do not have. MCA looks appropriate for binary data as its goal is to analyze relationships among several categorical variables, but all the examples I found of MCA involve reducing binary variables like present/absent into mere counts, which removes the nestedness/longitudinal component of my dataset. Dynamic factor analysis (DFA) has been described as "ordinary factor analysis, plus time dependence", but DFA examples I found consist of model building, which is not my current goal.

My main questions are:

1. Have I misunderstood or incorrectly described any of the above analyses (re: binary and longitudinal data)?
2. How should I achieve the goal of identifying groups/associations among longitudinal binary variables without reducing the dataframe?

Below creates an example dataframe:

library(tidyverse)
set.seed(123)
x_1 <- rbinom(n=100, size=1, prob=0.20)
x_2 <- rbinom(n=100, size=1, prob=0.25)
x_3 <- rbinom(n=100, size=1, prob=0.15)
y_4 <- rbinom(n=100, size=1, prob=0.80)
y_5 <- rbinom(n=100, size=1, prob=0.75)
y_6 <- rbinom(n=100, size=1, prob=0.75)
df <- data.frame(x_1,x_2,x_3,y_4,y_5,y_6) %>%
  mutate(time = c(rep(1:20), rep(1:25), rep(1:15), rep(1:15), rep(1:25))) %>%
  mutate(id = c(rep("001", 20), rep("002", 25), rep("003", 15), rep("004", 15), rep("005", 25)))