Multivariate Dynamics & Networks

Intensive Longitudinal Data

Overview

This tutorial covers two types of models for examining multivariate dynamics and networks: An N = 1 idiographic approach, and an N = all nomothetic approach.

For the N = 1 model we use a unified structural equation model (uSEM) that can be implemented in R using the pompom package.

For the N = all model we use a multilevel vector autoregression model (mlVAR) that can be implemented in R using the mlVAR package.

Outline

This script covers the following steps:

A. N = 1 model - unified Structural Equation Model (uSEM)

  • Set up the uSEM model

  • Check model summary (model fit statistics)

  • Data visualization and interpretation of results

B. N = all Multilevel VAR Model

  • Set up the mlVAR model

  • Check model summary(model fit statistics)

  • Data visualization and interpretation of results

C. Selected Readings

Preliminaries

Loading Libraries

# Check to see if necessary packages are installed, and install if not
packages <- c("psych", "pompom", "mlVAR", "ggplot2", "dplyr")
if (length(setdiff(packages, rownames(installed.packages()))) > 0) {
  install.packages(setdiff(packages, rownames(installed.packages())))  
}

# Load packages
library(psych)    #data description
library(ggplot2)  #data visualization
library(pompom)   #uSEM
library(mlVAR)    #mlVAR models
library(dplyr)    #data manipulation

Loading Data

We make use of the extended AMIB Data, a multiple time-scale data set that has been shared for teaching purposes. The data include person-level dispositions, daily diary assessments, and ecological momentary assessments that were obtained after everyday social interactions (event-contingent sampling).

The data are shared with intent to facilitate learning about analysis of intensive longitudinal data (from experience sampling, daily diary, or EMA designs). The data are posted only for teaching purposes and should not be used for research.

Use requires citation and acknowledgement of both this website and the following paper:

Ram, N., Conroy, D. E., Pincus, A. L., Hyde, A. L., & Molloy, L. E. (2012). Tethering theory to method: Using measures of intraindividual variability to operationalize individuals’ dynamic characteristics. In G. Hancock & J. Harring (Eds.), Advances in longitudinal methods in the social and behavioral sciences (pp. 81-110). New York: Information Age.

Specifically, the Phase 2 interaction-level data. This group of 30 participants provided up to 8 reports per day for 21 days. Thus, the time-series are relatively long (> 100 observations per person, M = 140, range = 64-168).

Loading interaction-level file (T = many, 43 assessments per person, on average).

#set filepath for data file
filepath <- "https://raw.githubusercontent.com/The-Change-Lab/collaborations/refs/heads/main/AMIB/AMIB_interactionP2.csv"

#read in the .csv file using the url() function
AMIB_interactionP2 <- read.csv(file=url(filepath), header=TRUE)

Subsetting to Variables of Interest

For the purpose of illustrating different models we use a not-small (to be truly multivariate) and not-large (for pedagogy) subset of variables:

  • Person index:

    • id

  • Time indices:

    • day (1-21)

    • interaction (1-168)

  • Six variables of interest:

    • igaff interpersonal affection

    • igdom interpersonal dominance

    • agval affect valence

    • agarous affect arousal

    • stress stress

    • health health[reverse coded]

Subsetting data to variables of interest:

#subsetting to vaiables of interest
AMIB_interactionP2 <- AMIB_interactionP2 %>%
  select(id, day, interaction, igaff, igdom, 
         agval, agarous, stress, health)

describe(AMIB_interactionP2)
##             vars    n   mean    sd median trimmed    mad min max range  skew
## id             1 4200 310.36 77.12    321  308.40 121.57 203 439   236  0.07
## day            2 4200  10.57  5.99     10   10.50   7.41   1  21    20  0.09
## interaction    3 4200  73.88 45.18     71   72.32  54.86   1 168   167  0.23
## igaff          4 4190   7.39  1.78      8    7.69   1.48   1   9     8 -1.43
## igdom          5 4190   6.92  1.62      7    7.08   1.48   1   9     8 -0.93
## agval          6 4189   6.80  2.09      7    7.08   2.22   1   9     8 -0.93
## agarous        7 4189   5.92  2.11      6    6.04   2.97   1   9     8 -0.51
## stress         8 4191   1.22  1.38      1    1.02   1.48   0   5     5  0.97
## health         9 4192   1.08  1.17      1    0.93   1.48   0   5     5  0.88
##             kurtosis   se
## id             -1.21 1.19
## day            -1.16 0.09
## interaction    -1.00 0.70
## igaff           2.01 0.03
## igdom           0.91 0.02
## agval           0.13 0.03
## agarous        -0.56 0.03
## stress          0.02 0.02
## health          0.00 0.02
#getting length of individual time-series
series_length <- AMIB_interactionP2 %>%
  group_by(id) %>%
  summarize(num_obs = length(unique(interaction)))
            
describe(series_length)
##         vars  n   mean    sd median trimmed    mad min max range  skew kurtosis
## id         1 30 306.97 79.59  318.5  304.04 124.54 203 439   236  0.14    -1.35
## num_obs    2 30 140.00 31.31  152.5  143.96  22.98  64 168   104 -0.75    -0.75
##            se
## id      14.53
## num_obs  5.72

Note that the variables are on different scales (e.g., 1-9, 0-5).

A. N = 1 model: unified structural equation model (uSEM)

The uSEM is a person-specific model, meaning that each person is modeled separately - one model for a person. We illustrate fitting for one person, and then indicate how we loop through many persons to obtain models for each person. Implementation is done using the pompom package by Xiao Yang, Nilam Ram, and Peter Molenaar: https://cran.r-project.org/web/packages/pompom/index.html

Subsetting to One Person’s Data (id 203)

#select one person's data
data_indiv <- AMIB_interactionP2 %>%
  filter(id == 203)

describe(data_indiv)
##             vars   n   mean    sd median trimmed   mad min max range  skew
## id             1 111 203.00  0.00    203  203.00  0.00 203 203     0   NaN
## day            2 111  10.39  5.94     10   10.25  7.41   1  21    20  0.19
## interaction    3 111  56.00 32.19     56   56.00 41.51   1 111   110  0.00
## igaff          4 111   8.13  0.53      8    8.15  0.00   7   9     2  0.11
## igdom          5 111   7.96  0.48      8    7.98  0.00   6   9     3 -0.53
## agval          6 111   8.13  0.53      8    8.14  0.00   7   9     2  0.12
## agarous        7 111   7.91  0.47      8    7.94  0.00   6   9     3 -0.70
## stress         8 111   0.15  0.49      0    0.01  0.00   0   3     3  3.56
## health         9 111   0.23  0.60      0    0.09  0.00   0   3     3  3.07
##             kurtosis   se
## id               NaN 0.00
## day            -1.15 0.56
## interaction    -1.23 3.06
## igaff           0.21 0.05
## igdom           3.10 0.05
## agval           0.30 0.05
## agarous         2.82 0.05
## stress         13.39 0.05
## health         10.05 0.06

Plotting Six-Variate Time-Series for One Person

#plotting intraindividual change 
data_indiv %>%
  ggplot(aes(x = interaction, group= id)) +
  #first variable
  geom_line(aes(y=igaff), color=1) +
  geom_line(aes(y=igdom), color=2) +
  geom_line(aes(y=agval), color=3) +
  geom_line(aes(y=agarous), color=4) +
  geom_line(aes(y=stress), color=5) +
  geom_line(aes(y=health), color=6) +
  #plot layouts
  scale_x_continuous(name="Interaction #") +
  scale_y_continuous(name="Raw Values") +  
  theme_classic() +
  theme(axis.title=element_text(size=14),
        axis.text=element_text(size=14),
        plot.title=element_text(size=14, hjust=.5)) +
  ggtitle("ID # 203")

Here it is easy to note that the variables are on different scales (e.g., 1-9, 0-5).

Standardize the Data

For the uSEM, we standardize them all within-person to keep focus on dynamics. Importantly, the data are considered as stationary time-series that fluctuate around zero.

Standardizing within-person:

# standardize specific data columns (not the id or time variables in first 3 columns)
data_indiv[c(4:9)] <- lapply(data_indiv[c(4:9)], 
                             function(x) c(scale(x, center=TRUE, scale=TRUE)))
describe(data_indiv)
##             vars   n   mean    sd median trimmed   mad    min    max  range
## id             1 111 203.00  0.00 203.00  203.00  0.00 203.00 203.00   0.00
## day            2 111  10.39  5.94  10.00   10.25  7.41   1.00  21.00  20.00
## interaction    3 111  56.00 32.19  56.00   56.00 41.51   1.00 111.00 110.00
## igaff          4 111   0.00  1.00  -0.25    0.03  0.00  -2.13   1.64   3.78
## igdom          5 111   0.00  1.00   0.08    0.03  0.00  -4.09   2.16   6.25
## agval          6 111   0.00  1.00  -0.25    0.02  0.00  -2.15   1.65   3.81
## agarous        7 111   0.00  1.00   0.20    0.07  0.00  -4.02   2.31   6.32
## stress         8 111   0.00  1.00  -0.31   -0.29  0.00  -0.31   5.81   6.12
## health         9 111   0.00  1.00  -0.39   -0.24  0.00  -0.39   4.59   4.98
##              skew kurtosis   se
## id            NaN      NaN 0.00
## day          0.19    -1.15 0.56
## interaction  0.00    -1.23 3.06
## igaff        0.11     0.21 0.09
## igdom       -0.53     3.10 0.09
## agval        0.12     0.30 0.09
## agarous     -0.70     2.82 0.09
## stress       3.56    13.39 0.09
## health       3.07    10.05 0.09

Plotting Six-Variate Time-series for One Person (Standardized Vars)

#plotting intraindividual change 
data_indiv %>%
  ggplot(aes(x = interaction, group= id)) +
  #first variable
  geom_line(aes(y=igaff), color=1) +
  geom_line(aes(y=igdom), color=2) +
  geom_line(aes(y=agval), color=3) +
  geom_line(aes(y=agarous), color=4) +
  geom_line(aes(y=stress), color=5) +
  geom_line(aes(y=health), color=6) +
  #plot layouts
  scale_x_continuous(name="Interaction #") +
  scale_y_continuous(name="Standardized Values") +  
  theme_classic() +
  theme(axis.title=element_text(size=14),
        axis.text=element_text(size=14),
        plot.title=element_text(size=14, hjust=.5)) +
  ggtitle("ID # 203")

Now we see that all the variables are in standardized form.

Check the Data

It is useful to check that there are data in all columns. If any one of the variables is all missing (or has no variance), the model cannot be fit. Missing data on a few observations within a column is ok.

# check column missing
na_col <- 0
for (col in 1:ncol(data_indiv)) {
  if (sum(is.na(data_indiv[,col])) == nrow(data_indiv)){
        na_col <- na_col + 1
  }
}
na_col
## [1] 0

All columns are reported.

Fitting the uSEM Model

In pompom syntax, we need to specify how many variables are included in the analysis (var.number), the dataset without the id and time variables (data), the number of lags to examine (lag.order), whether we want intermediate output (verbose), and whether to remove non-significant coefficients in the last iteration of model construction (trim).

Fitting uSEM model to individual data

usem_indiv <- uSEM(var.number = 6, 
               data = data_indiv[ ,c(4:9)], 
               lag.order = 1, 
               verbose = FALSE,
               trim = FALSE)

For these data, this only took a few seconds.

Save the uSEM Model Summary in an Object

The model summary is put into an object for easy access to the estimated coefficients, the standard errors of the coefficients, and the model fit indices.

Obtaining model summary.

ms_indiv <- model_summary(model.fit = usem_indiv,
                          var.number = 6, 
                          lag.order = 1)

ms_indiv
## $beta
##            [,1]        [,2]      [,3]      [,4]      [,5]      [,6]       [,7]
##  [1,] 0.0000000  0.00000000 0.0000000 0.0000000 0.0000000 0.0000000  0.0000000
##  [2,] 0.0000000  0.00000000 0.0000000 0.0000000 0.0000000 0.0000000  0.0000000
##  [3,] 0.0000000  0.00000000 0.0000000 0.0000000 0.0000000 0.0000000  0.0000000
##  [4,] 0.0000000  0.00000000 0.0000000 0.0000000 0.0000000 0.0000000  0.0000000
##  [5,] 0.0000000  0.00000000 0.0000000 0.0000000 0.0000000 0.0000000  0.0000000
##  [6,] 0.0000000  0.00000000 0.0000000 0.0000000 0.0000000 0.0000000  0.0000000
##  [7,] 0.2179232  0.00000000 0.0000000 0.0000000 0.0000000 0.0000000  0.0000000
##  [8,] 0.0000000 -0.06241162 0.0000000 0.1902340 0.0000000 0.0000000  0.5658826
##  [9,] 0.0000000  0.00000000 0.2048743 0.0000000 0.0000000 0.0000000  0.3757771
## [10,] 0.0000000  0.00000000 0.0000000 0.1843903 0.0000000 0.0000000  0.3460383
## [11,] 0.0000000  0.00000000 0.0000000 0.0000000 0.4339751 0.0000000 -0.1891902
## [12,] 0.0000000  0.00000000 0.0000000 0.0000000 0.0000000 0.6903588 -0.1366965
##       [,8] [,9] [,10] [,11] [,12]
##  [1,]    0    0     0     0     0
##  [2,]    0    0     0     0     0
##  [3,]    0    0     0     0     0
##  [4,]    0    0     0     0     0
##  [5,]    0    0     0     0     0
##  [6,]    0    0     0     0     0
##  [7,]    0    0     0     0     0
##  [8,]    0    0     0     0     0
##  [9,]    0    0     0     0     0
## [10,]    0    0     0     0     0
## [11,]    0    0     0     0     0
## [12,]    0    0     0     0     0
## 
## $beta.se
##             [,1]       [,2]     [,3]       [,4]       [,5]       [,6]
##  [1,] 0.00000000 0.00000000 0.000000 0.00000000 0.00000000 0.00000000
##  [2,] 0.00000000 0.00000000 0.000000 0.00000000 0.00000000 0.00000000
##  [3,] 0.00000000 0.00000000 0.000000 0.00000000 0.00000000 0.00000000
##  [4,] 0.00000000 0.00000000 0.000000 0.00000000 0.00000000 0.00000000
##  [5,] 0.00000000 0.00000000 0.000000 0.00000000 0.00000000 0.00000000
##  [6,] 0.00000000 0.00000000 0.000000 0.00000000 0.00000000 0.00000000
##  [7,] 0.09186633 0.00000000 0.000000 0.00000000 0.00000000 0.00000000
##  [8,] 0.00000000 0.08313478 0.000000 0.08281152 0.00000000 0.00000000
##  [9,] 0.00000000 0.00000000 0.085359 0.00000000 0.00000000 0.00000000
## [10,] 0.00000000 0.00000000 0.000000 0.08791538 0.00000000 0.00000000
## [11,] 0.00000000 0.00000000 0.000000 0.00000000 0.08416846 0.00000000
## [12,] 0.00000000 0.00000000 0.000000 0.00000000 0.00000000 0.06710492
##             [,7] [,8] [,9] [,10] [,11] [,12]
##  [1,] 0.00000000    0    0     0     0     0
##  [2,] 0.00000000    0    0     0     0     0
##  [3,] 0.00000000    0    0     0     0     0
##  [4,] 0.00000000    0    0     0     0     0
##  [5,] 0.00000000    0    0     0     0     0
##  [6,] 0.00000000    0    0     0     0     0
##  [7,] 0.00000000    0    0     0     0     0
##  [8,] 0.07904285    0    0     0     0     0
##  [9,] 0.08640937    0    0     0     0     0
## [10,] 0.08900623    0    0     0     0     0
## [11,] 0.08519014    0    0     0     0     0
## [12,] 0.06790280    0    0     0     0     0
## 
## $psi
##             [,1]        [,2]        [,3]        [,4]       [,5]        [,6]
##  [1,]  0.9994419  0.54618258  0.43583469  0.33952881 -0.1836394 -0.18315409
##  [2,]  0.5461826  0.99994828  0.23527287  0.34397710 -0.1307191 -0.08030321
##  [3,]  0.4358347  0.23527287  0.99943356  0.28759687 -0.1851460 -0.09846587
##  [4,]  0.3395288  0.34397710  0.28759687  0.99963571 -0.2103505 -0.08011198
##  [5,] -0.1836394 -0.13071910 -0.18514600 -0.21035052  0.9991041  0.09186890
##  [6,] -0.1831541 -0.08030321 -0.09846587 -0.08011198  0.0918689  0.99861280
##  [7,]  0.0000000  0.00000000  0.00000000  0.00000000  0.0000000  0.00000000
##  [8,]  0.0000000  0.00000000  0.00000000  0.00000000  0.0000000  0.00000000
##  [9,]  0.0000000  0.00000000  0.00000000  0.00000000  0.0000000  0.00000000
## [10,]  0.0000000  0.00000000  0.00000000  0.00000000  0.0000000  0.00000000
## [11,]  0.0000000  0.00000000  0.00000000  0.00000000  0.0000000  0.00000000
## [12,]  0.0000000  0.00000000  0.00000000  0.00000000  0.0000000  0.00000000
##            [,7]      [,8]     [,9]     [,10]     [,11]     [,12]
##  [1,] 0.0000000 0.0000000 0.000000 0.0000000 0.0000000 0.0000000
##  [2,] 0.0000000 0.0000000 0.000000 0.0000000 0.0000000 0.0000000
##  [3,] 0.0000000 0.0000000 0.000000 0.0000000 0.0000000 0.0000000
##  [4,] 0.0000000 0.0000000 0.000000 0.0000000 0.0000000 0.0000000
##  [5,] 0.0000000 0.0000000 0.000000 0.0000000 0.0000000 0.0000000
##  [6,] 0.0000000 0.0000000 0.000000 0.0000000 0.0000000 0.0000000
##  [7,] 0.9278185 0.0000000 0.000000 0.0000000 0.0000000 0.0000000
##  [8,] 0.0000000 0.6636745 0.000000 0.0000000 0.0000000 0.0000000
##  [9,] 0.0000000 0.0000000 0.761397 0.0000000 0.0000000 0.0000000
## [10,] 0.0000000 0.0000000 0.000000 0.8498415 0.0000000 0.0000000
## [11,] 0.0000000 0.0000000 0.000000 0.0000000 0.7784097 0.0000000
## [12,] 0.0000000 0.0000000 0.000000 0.0000000 0.0000000 0.4926709
## 
## $chisq
##    chisq 
## 37.16166 
## 
## $cfi
## cfi 
##   1 
## 
## $tli
##      tli 
## 1.011639 
## 
## $rmsea
## rmsea 
##     0 
## 
## $srmr
##       srmr 
## 0.04684503

We see all the matrices of the model as well as some summary fit statistics.

Check Fit Indices

To assess if the model reached a good fit we use a “3 out of 4” rule, meaning that if at least 3 of the 4 following fit rules are satisfied, the model is considered as being a good fitting model. Note that this is an emerging approach in the person-specific uSEM literature (not a standard rule for all SEM),

  1. CFI > .95

  2. TLI > .95

  3. RMSEA < .08

  4. SRMR < .08

Checking if the 3 of 4 Rule is Satisfied

#check for good fit
print(sum(ms_indiv$cfi > .95, 
      ms_indiv$tli > .95,
      ms_indiv$rmsea < .08,
      ms_indiv$srmr < .08) >= 3)
## [1] TRUE

Becasue the model fitting is iterative and continues until the data are well-represented, we expect to achieve good fit. Lack of good fit may indicate some data problems.

Visualizing the uSEM model

While the matrices above are extremely useful for understanding the model, network visualizations often facilitate interpretation.

Plot model as a network graph.

#plot model as network
plot_network_graph(ms_indiv$beta, var.number = 6)

## NULL

Nodes 1 to 6 represent the six variables. Red/green edges indicate negative/positive temporal relations (or coefficients). Dashed edges indicate lag-1 relations and solid edges indicate contemporaneous relations. The width of the edge indicates the absolute value of the relation.

Interpretation

When the dimension of the network is larger than 2 (nodes), it is often difficult to interpret the overall dynamics of the network as focus on any single pair of nodes may ignore other connections in the network. Therefore, we have developed new network metrics to examine how the network behaves. In particular we examine recovery time of each node, which describes how fast that variable will return to equilibrium (a stable state) after a perturbation in that node or any other node. We set a threshold for what is considered close-enough to equilibrium (threshold), and can evaluate the recovery times using a bootstrap approach that accommodates uncertainty in the model parameters.

Compute and Examine the Recovery Times Implied by the Estimated Network of Relations

#compute recovery time for all nodes
recovery_time <- iRAM(model.fit = usem_indiv, 
     beta = ms_indiv$beta,
     var.number = 6,
     lag.order = 1, 
     threshold = 0.01, 
     boot = TRUE,
     replication = 200, steps = 100)

#examine mean across the recovery time of all node-pairs
recovery_time$mean
##       [,1]  [,2]  [,3]  [,4] [,5]   [,6]
## [1,] 4.545 4.665 4.815 4.625 5.66  9.240
## [2,] 0.000 3.275 0.000 0.000 0.00  0.000
## [3,] 0.000 0.000 4.440 0.000 0.00  0.000
## [4,] 0.000 4.015 0.000 4.195 0.00  0.000
## [5,] 0.000 0.000 0.000 0.000 7.14  0.000
## [6,] 0.000 0.000 0.000 0.000 0.00 14.455

The recovery times can also be plotted for easier interpretation.

Plot the Time-Profiles for When Each Node Gets a Perturbation (xupper controls length of time axis)

plot_time_profile(recovery_time$time.profile.data, 
                  var.number = 6,
                  threshold = 0.01,
                  xupper = 20)

## NULL

Interpretation of the Recovery Time of Node-Pairs and the Structure of the Network

We can see most nodes have a recovery time of 3 ~ 4 steps, but node5 (“stress”) and node6 (“health”) have longer recovery times (7.1 and 14.5 time steps, respectively). Looking at the network graph, we see stronger self-loops for node5 and node6 (thickness of lines indicates size of auto-regression). These parameters indicate that there is lots of “carry-over” across observations and thus the length of time it takes this variable to return to equilibrium is extended. In this case, the interpretation is that “stress” and “health” operate at a slower time-scale than the “interpersonal behavior” and “affect” variables. This makes sense!

Further details of how recovery time is computed and interpreted can be found in the selected readings.

Looping to Fit Models for Many Persons

Now we can run the model in a loop for multiple participants. Let us try 5 participants…

Preparing and running loop.

#create list of ids
id_list <- unique(AMIB_interactionP2$id)

for (id in id_list[1:5]) # select the first 5 participants as examples
{
  #select data for individual
  data_indiv <- AMIB_interactionP2[AMIB_interactionP2$id == id,]
  #standardize data for individual
  data_indiv[c(4:9)] <- lapply(data_indiv[c(4:9)], 
                             function(x) c(scale(x, center=TRUE, scale=TRUE)))
  #indicate id
  print(id)
  
  # check column missing
  na.col <- 0
  for (col in 1:ncol(data_indiv))
  {
    if (sum(is.na(data_indiv[,col])) == nrow(data_indiv)){
          na.col <- na.col + 1
    }
  }
  #fit uSEM model
  if (na.col == 0) # no NA column
  {
      usem_indiv <- uSEM(var.number = 6, 
                   data = data_indiv[ ,c(4:9)], 
                   lag.order = 1, 
                   verbose = FALSE,
                   trim = FALSE)
  
     ms_indiv <- model_summary(
       model.fit = usem_indiv,
       var.number = 6, 
       lag.order = 1)
     
     print(sum(ms_indiv$cfi > .95, 
           ms_indiv$tli > .95,
           ms_indiv$rmsea < .08,
           ms_indiv$srmr < .08) >= 3)

     plot_network_graph(ms_indiv$beta, var.number = 6)
  }
}
## [1] 203
## [1] TRUE

## [1] 204
## [1] TRUE

## [1] 205
## [1] TRUE

## [1] 208
## [1] TRUE

## [1] 211
## [1] TRUE

We see lots of heterogeneity in the structure of these 5 individuals’ temporal dynamics!

B. N = all model: multilevel vector autoregression model (mlVAR)

Now we turn to fitting of a multilevel VAR model.

The multilevel VAR is a group-level model, in contrast to the uSEM which is a person-specific model. Thus, multiple individuals are fit together in one model.

Describing the Data

describe(AMIB_interactionP2)
##             vars    n   mean    sd median trimmed    mad min max range  skew
## id             1 4200 310.36 77.12    321  308.40 121.57 203 439   236  0.07
## day            2 4200  10.57  5.99     10   10.50   7.41   1  21    20  0.09
## interaction    3 4200  73.88 45.18     71   72.32  54.86   1 168   167  0.23
## igaff          4 4190   7.39  1.78      8    7.69   1.48   1   9     8 -1.43
## igdom          5 4190   6.92  1.62      7    7.08   1.48   1   9     8 -0.93
## agval          6 4189   6.80  2.09      7    7.08   2.22   1   9     8 -0.93
## agarous        7 4189   5.92  2.11      6    6.04   2.97   1   9     8 -0.51
## stress         8 4191   1.22  1.38      1    1.02   1.48   0   5     5  0.97
## health         9 4192   1.08  1.17      1    0.93   1.48   0   5     5  0.88
##             kurtosis   se
## id             -1.21 1.19
## day            -1.16 0.09
## interaction    -1.00 0.70
## igaff           2.01 0.03
## igdom           0.91 0.02
## agval           0.13 0.03
## agarous        -0.56 0.03
## stress          0.02 0.02
## health          0.00 0.02

From a multilevel perspective, the data are considered as 4200 observations nested within 30 persons that are considered as replicates drawn from a population.

Plotting Six-variate Time-series for All Persons

#plotting intraindividual change 
AMIB_interactionP2 %>%
  ggplot(aes(x = interaction, group= id)) +
  #first variable
  geom_line(aes(y=igaff), color=1) +
  geom_line(aes(y=igdom), color=2) +
  geom_line(aes(y=agval), color=3) +
  geom_line(aes(y=agarous), color=4) +
  geom_line(aes(y=stress), color=5) +
  geom_line(aes(y=health), color=6) +
  #plot layouts
  scale_x_continuous(name="Interaction #") +
  scale_y_continuous(name="Raw Values") +  
  theme_classic() +
  theme(axis.title=element_text(size=10),
        axis.text=element_text(size=10),
        plot.title=element_text(size=12, hjust=.5)) +
  facet_wrap(~id, ncol=2) +
  ggtitle("AMIB Phase 2 Data (6 vars)")

We see there is some difficulty in having so many persons with so much data =:-]

Not Standardizing the Data

For the mlVAR, the data are kept in their original form. The estimation package does some standardization, as well as separation of within-person and between-person components in the background. Importantly, the data are considered as stationary time-series that fluctuate around individuals’ mean-level.

Fitting the mlVAR model

In mlVAR syntax, we need to specify the dataset (data), the specific variables included in the analysis (vars), the id variable (idvar), the number of lags to examine (lags), and the time indices (dayvar, beepvar). The two time variables are used to accommodate unequal spacing in EMA studies where no observations are obtained during the nighttime (e.g., the last evening beep is not considered to be 1-lag away from the first morning beep).

Fitting the mlVAR model to all N = 30 time-series data:

#fitting mlVAR model
mlvar_all <- mlVAR(data = AMIB_interactionP2, 
                     vars = c("igaff", "igdom", "agval", "agarous", "stress", "health"),
                     idvar = "id",
                     lags = 1, 
                     dayvar = "day", 
                     beepvar = "interaction")
##   |                                                                              |                                                                      |   0%  |                                                                              |============                                                          |  17%  |                                                                              |=======================                                               |  33%  |                                                                              |===================================                                   |  50%  |                                                                              |===============================================                       |  67%  |                                                                              |==========================================================            |  83%  |                                                                              |======================================================================| 100%
##   |                                                                              |                                                                      |   0%  |                                                                              |============                                                          |  17%  |                                                                              |=======================                                               |  33%  |                                                                              |===================================                                   |  50%  |                                                                              |===============================================                       |  67%  |                                                                              |==========================================================            |  83%  |                                                                              |======================================================================| 100%
##   |                                                                              |                                                                      |   0%  |                                                                              |==                                                                    |   3%  |                                                                              |=====                                                                 |   7%  |                                                                              |=======                                                               |  10%  |                                                                              |=========                                                             |  13%  |                                                                              |============                                                          |  17%  |                                                                              |==============                                                        |  20%  |                                                                              |================                                                      |  23%  |                                                                              |===================                                                   |  27%  |                                                                              |=====================                                                 |  30%  |                                                                              |=======================                                               |  33%  |                                                                              |==========================                                            |  37%  |                                                                              |============================                                          |  40%  |                                                                              |==============================                                        |  43%  |                                                                              |=================================                                     |  47%  |                                                                              |===================================                                   |  50%  |                                                                              |=====================================                                 |  53%  |                                                                              |========================================                              |  57%  |                                                                              |==========================================                            |  60%  |                                                                              |============================================                          |  63%  |                                                                              |===============================================                       |  67%  |                                                                              |=================================================                     |  70%  |                                                                              |===================================================                   |  73%  |                                                                              |======================================================                |  77%  |                                                                              |========================================================              |  80%  |                                                                              |==========================================================            |  83%  |                                                                              |=============================================================         |  87%  |                                                                              |===============================================================       |  90%  |                                                                              |=================================================================     |  93%  |                                                                              |====================================================================  |  97%  |                                                                              |======================================================================| 100%

In the background the multivariate model is fit as a collection of multilevel models (using lmer) and the outputs are combined together to obtain the full set of multivariate within-person and between-person relations. For these data, the estimation only took a couple of minutes.

Examine the model summary.

# Summary of all parameter estimates:
summary(mlvar_all)
## 
## mlVAR estimation completed. Input was:
##      - Variables: igaff igdom agval agarous stress health 
##      - Lags: 1 
##      - Estimator: lmer 
##      - Temporal: correlated
## 
## Information indices:
##      var      aic      bic
##    igaff 9539.195 9792.416
##    igdom 9603.567 9856.787
##    agval 8805.919 9059.140
##  agarous 8843.996 9097.217
##   stress 6603.191 6856.411
##   health 5759.345 6012.565
## 
## 
## Temporal effects:
##     from      to lag  fixed    SE     P ran_SD
##    igaff   igaff   1  0.050 0.027 0.070  0.097
##    igaff   igdom   1 -0.019 0.023 0.417  0.060
##    igaff   agval   1 -0.009 0.020 0.661  0.044
##    igaff agarous   1 -0.039 0.018 0.031  0.026
##    igaff  stress   1  0.017 0.020 0.394  0.075
##    igaff  health   1  0.027 0.014 0.049  0.037
##    igdom   igaff   1 -0.025 0.018 0.181  0.027
##    igdom   igdom   1  0.034 0.026 0.192  0.096
##    igdom   agval   1 -0.036 0.016 0.025  0.017
##    igdom agarous   1  0.036 0.017 0.039  0.034
##    igdom  stress   1  0.026 0.012 0.028  0.015
##    igdom  health   1  0.006 0.012 0.596  0.028
##    agval   igaff   1 -0.006 0.025 0.815  0.060
##    agval   igdom   1  0.015 0.029 0.603  0.089
##    agval   agval   1  0.146 0.034 0.000  0.136
##    agval agarous   1  0.026 0.022 0.246  0.042
##    agval  stress   1  0.029 0.017 0.086  0.041
##    agval  health   1  0.004 0.014 0.769  0.023
##  agarous   igaff   1 -0.008 0.029 0.769  0.115
##  agarous   igdom   1  0.058 0.026 0.028  0.094
##  agarous   agval   1  0.013 0.019 0.489  0.046
##  agarous agarous   1  0.187 0.026 0.000  0.102
##  agarous  stress   1 -0.020 0.015 0.207  0.049
##  agarous  health   1 -0.004 0.014 0.767  0.044
##   stress   igaff   1 -0.067 0.023 0.004  0.023
##   stress   igdom   1  0.003 0.026 0.904  0.065
##   stress   agval   1 -0.174 0.032 0.000  0.124
##   stress agarous   1  0.052 0.041 0.207  0.186
##   stress  stress   1  0.641 0.043 0.000  0.209
##   stress  health   1  0.048 0.018 0.008  0.061
##   health   igaff   1 -0.010 0.031 0.737  0.110
##   health   igdom   1  0.021 0.028 0.459  0.087
##   health   agval   1 -0.054 0.026 0.036  0.083
##   health agarous   1 -0.084 0.026 0.001  0.083
##   health  stress   1  0.069 0.021 0.001  0.072
##   health  health   1  0.657 0.037 0.000  0.182
## 
## 
## Contemporaneous effects (posthoc estimated):
##       v1      v2 P 1->2 P 1<-2   pcor ran_SD_pcor    cor ran_SD_cor
##    igdom   igaff  0.000  0.000  0.184       0.175  0.250      0.211
##    agval   igaff  0.000  0.000  0.385       0.131  0.461      0.153
##    agval   igdom  0.014  0.086  0.057       0.081  0.167      0.140
##  agarous   igaff  0.206  0.073  0.041       0.074  0.102      0.100
##  agarous   igdom  0.000  0.000  0.227       0.080  0.248      0.103
##  agarous   agval  0.456  0.965  0.014       0.129  0.063      0.152
##   stress   igaff  0.002  0.036 -0.060       0.058 -0.239      0.145
##   stress   igdom  0.501  0.606  0.013       0.045 -0.059      0.091
##   stress   agval  0.000  0.000 -0.343       0.151 -0.426      0.174
##   stress agarous  0.001  0.004  0.074       0.058  0.023      0.091
##   health   igaff  0.496  0.614  0.010       0.012 -0.099      0.119
##   health   igdom  0.655  0.622 -0.010       0.032 -0.063      0.057
##   health   agval  0.002  0.000 -0.090       0.068 -0.195      0.149
##   health agarous  0.000  0.000 -0.105       0.062 -0.108      0.059
##   health  stress  0.000  0.000  0.188       0.209  0.244      0.222
## 
## 
## Between-subject effects:
##       v1      v2 P 1->2 P 1<-2   pcor    cor
##    igdom   igaff  0.053  0.012  0.320  0.420
##    agval   igaff  0.000  0.000  0.657  0.797
##    agval   igdom  0.201  0.333 -0.182  0.255
##  agarous   igaff  0.534  0.439 -0.027  0.327
##  agarous   igdom  0.000  0.013  0.472  0.532
##  agarous   agval  0.267  0.267  0.169  0.295
##   stress   igaff  0.511  0.456  0.106 -0.480
##   stress   igdom  0.583  0.776  0.071 -0.105
##   stress   agval  0.252  0.091 -0.223 -0.612
##   stress agarous  0.741  0.778  0.050 -0.100
##   health   igaff  0.013  0.845 -0.172 -0.612
##   health   igdom  0.695  0.494 -0.090 -0.233
##   health   agval  0.361  0.252 -0.161 -0.679
##   health agarous  0.886  0.651  0.022 -0.184
##   health  stress  0.000  0.000  0.640  0.780

We see some summary fit statistics, and all the matrices representing the different parts of the model.

Visualizing the mlVAR model

While the matrices above are extremely useful for understanding the model, network visualizations often facilitate interpretation - and have undoubtedly contributed to the quick uptake of this method.

Plot the Group-Level Networks for Within-Person and Between-Person Relations

# Plot temporal relations:
plot(mlvar_all, "temporal", 
     title = "Within-person temporal (lag-1) relations", 
     layout = "circle", 
     nonsig = "hide")

# Plot contemporaneous partial correlations:
plot(mlvar_all, "contemporaneous", 
     title = "Within-person contemporaneous relations",
     layout = "circle", 
     nonsig = "hide")

# Plot between-persons partial correlations:
plot(mlvar_all, "between", 
     title = "Between-person relations", 
     layout = "circle", 
     nonsig = "hide")

Interpretation

Temporal network (within-person lag-1 relations)

For the prototypical person, health and stress have strong positive auto-regression, meaning they tend to persist over time. stress also has a negative forward influence on agval (affect valence), meaning when stress is higher, affect valence goes lower.

Contemporaneous network (within-person contemporaneous relations)

For the prototypical person, many variables tend to move together, as indicated by the number and thickness of the edges that exist in the network. For example, igaff (interpersonal affiliation) and agval (affect valence) tend to fluctuate together (green positive edge), while stress and agval tend to fluctuate in opposition to one another (red negative edge).

Between-person network (between-person relations)

Across persons, there are some variables covary. For example, individuals with high average igaff (interpersonal affiliation) also tend to have high average agval (affect valence) (green positive edge) and low average health (lower score is better health) (red negative edge).

Conclusion

This tutorial illustrates two different approaches for examining multivariate dynamics - a person-specific (N = 1) network approach based on the unified structural equation model (uSEM, implemented in R using the pompom package); and a group-based (N = all) network approach based on the multilevel vector autoregression model (mlVAR, implemented in R using the mlVAR package).

As we get the requisite data (long time-series), we look forward making use of these approaches, and understanding more about which approach is useful in which setting and for which research questions.

Thanks for being multivariate and dynamic as we all figure it out!

Citations

Bringmann, L. F., Vissers, N., Wichers, M., Geschwind, N., Kuppens, P., Peeters, F., Borsboom, D., & Tuerlinckx, F. (2013). A Network Approach to Psychopathology: New Insights into Clinical Longitudinal Data. PLOS ONE, 8(4), e60188. https://doi.org/10.1371/journal.pone.0060188

Epskamp, S., Deserno, M. K., & Bringmann, L. F. (2024). mlVAR: Multi-Level Vector Autoregression (Version 0.5.2) [Computer software]. https://CRAN.R-project.org/package=mlVAR

R Core Team. (2024). R: A Language and Environment for Statistical Computing. Foundation for Statistical Computing. https://www.R-project.org/

Ram, N., Conroy, D. E., Pincus, A. L., Hyde, A. L., & Molloy, L. E. (2012). Tethering theory to method: Using measures of intraindividual variability to operationalize individuals’ dynamic characteristics. In G. Hancock & J. Harring (Eds.), Advances in longitudinal methods in the social and behavioral sciences (pp. 81-110). New York: Information Age.

Revelle, W. (2024). psych: Procedures for Psychological, Psychometric, and Personality Research. Northwestern University. https://CRAN.R-project.org/package=psych

Wickham, H. (2016). ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag. https://ggplot2.tidyverse.org/

Wickham, H., François, R., Henry, L., Müller, K., & Vaughan, D. (2023). dplyr: A Grammar of Data Manipulation (Version 1.1.4). https://CRAN.R-project.org/package=dplyr

Yang, X., Ram, N., Gest, S. D., Lydon-Staley, D. M., Conroy, D. E., Pincus, A. L., & Molenaar, P. C. M. (2018). Socioemotional Dynamics of Emotion Regulation and Depressive Symptoms: A Person-Specific Network Approach. Complexity, 2018(1), 5094179. https://doi.org/10.1155/2018/5094179

Yang, X., Ram, N., & Molenaar, P. C. M. (2021). pompom: Person-Oriented Method and Perturbation on the Model (Version 0.2.1). https://CRAN.R-project.org/package=pompom