APA-ATI Intensive Longitudinal Data: Session R
Multivariate Dynamics & Networks
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 multievel 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
D. Selected Readings
for uSEM model: Yang, X., Ram, N., Gest, S., Lydon, D., 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, Article ID 5094179. doi: 10.1155/2018/5094179 [Open Access https://www.hindawi.com/journals/complexity/2018/5094179/]
for mlVAR model: Bringmann LF, Vissers N, Wichers M, Geschwind N, Kuppens P, Peeters F, et al. (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 [Open Access https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0060188]
Preliminaries
Loading Libraries
Loading libraries used in this script.
# Check to see if necessary packages are installed, and install if not
packages <- c("psych", "pompom", "mlVAR")
if (length(setdiff(packages, rownames(installed.packages()))) > 0) {
install.packages(setdiff(packages, rownames(installed.packages())))
}
# Load packages
library(psych) #for data description
library(plyr) #for data manipulation
library(ggplot2) #for data visualization
library(pompom) #for uSEM
library(mlVAR) #for mlVAR modelsLoading Data
We make use of the extended AMIB data. 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://quantdev.ssri.psu.edu/sites/qdev/files/AMIBshare_phase2_interaction_2019_0501.csv"
#read in the .csv file using the url() function
AMIB_interactionP2 <- read.csv(file=url(filepath),header=TRUE)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) and “interaction” (1-168) Six variables of interest: “igaff”, “igdom”, “agval”, “agarous”“,”stress“, and”health" (interpersonal affection, interpersonal dominance, affect valence, affect arousal, stress, health[reverse coded])
Subsetting data to variables of interest.
#subsetting to vaiables of interest
AMIB_interactionP2 <- AMIB_interactionP2[ ,c("id","day","interaction",
"igaff","igdom",
"agval", "agarous",
"stress", "health")]
describe(AMIB_interactionP2)## vars n mean sd median trimmed mad min max range
## id 1 4200 310.36 77.12 321 308.40 121.57 203 439 236
## day 2 4200 10.57 5.99 10 10.50 7.41 1 21 20
## interaction 3 4200 73.88 45.18 71 72.32 54.86 1 168 167
## igaff 4 4190 7.39 1.78 8 7.69 1.48 1 9 8
## igdom 5 4190 6.92 1.62 7 7.08 1.48 1 9 8
## agval 6 4189 6.80 2.09 7 7.08 2.22 1 9 8
## agarous 7 4189 5.92 2.11 6 6.04 2.97 1 9 8
## stress 8 4191 1.22 1.38 1 1.02 1.48 0 5 5
## health 9 4192 1.08 1.17 1 0.93 1.48 0 5 5
## skew kurtosis se
## id 0.07 -1.21 1.19
## day 0.09 -1.16 0.09
## interaction 0.23 -1.00 0.70
## igaff -1.43 2.01 0.03
## igdom -0.93 0.91 0.02
## agval -0.93 0.13 0.03
## agarous -0.51 -0.56 0.03
## stress 0.97 0.02 0.02
## health 0.88 0.00 0.02#getting length of individual time-series
# series_length <- ddply(AMIB_interactionP2, "id", summarize,
# num_obs = length(unique(interaction)))
# describe(series_length)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[AMIB_interactionP2$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.06Plotting six-variate time-series for one person
#plotting intraindividual change
ggplot(data = data_indiv,
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
## id 1 111 203.00 0.00 203.00 203.00 0.00 203.00 203.00
## day 2 111 10.39 5.94 10.00 10.25 7.41 1.00 21.00
## interaction 3 111 56.00 32.19 56.00 56.00 41.51 1.00 111.00
## igaff 4 111 0.00 1.00 -0.25 0.03 0.00 -2.13 1.64
## igdom 5 111 0.00 1.00 0.08 0.03 0.00 -4.09 2.16
## agval 6 111 0.00 1.00 -0.25 0.02 0.00 -2.15 1.65
## agarous 7 111 0.00 1.00 0.20 0.07 0.00 -4.02 2.31
## stress 8 111 0.00 1.00 -0.31 -0.29 0.00 -0.31 5.81
## health 9 111 0.00 1.00 -0.39 -0.24 0.00 -0.39 4.59
## range skew kurtosis se
## id 0.00 NaN NaN 0.00
## day 20.00 0.19 -1.15 0.56
## interaction 110.00 0.00 -1.23 3.06
## igaff 3.78 0.11 0.21 0.09
## igdom 6.25 -0.53 3.10 0.09
## agval 3.81 0.12 0.30 0.09
## agarous 6.32 -0.70 2.82 0.09
## stress 6.12 3.56 13.39 0.09
## health 4.98 3.07 10.05 0.09Plotting six-variate time-series for one person
#plotting intraindividual change
ggplot(data = data_indiv,
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] 0All 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 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]
## [1,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [2,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [3,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [4,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [5,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [6,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [7,] 0.2178887 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [8,] 0.0000000 -0.0623992 0.0000000 0.1902141 0.0000000 0.0000000
## [9,] 0.0000000 0.0000000 0.2048157 0.0000000 0.0000000 0.0000000
## [10,] 0.0000000 0.0000000 0.0000000 0.1843832 0.0000000 0.0000000
## [11,] 0.0000000 0.0000000 0.0000000 0.0000000 0.4339714 0.0000000
## [12,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.6903564
## [,7] [,8] [,9] [,10] [,11] [,12]
## [1,] 0.0000000 0 0 0 0 0
## [2,] 0.0000000 0 0 0 0 0
## [3,] 0.0000000 0 0 0 0 0
## [4,] 0.0000000 0 0 0 0 0
## [5,] 0.0000000 0 0 0 0 0
## [6,] 0.0000000 0 0 0 0 0
## [7,] 0.0000000 0 0 0 0 0
## [8,] 0.5657580 0 0 0 0 0
## [9,] 0.3759418 0 0 0 0 0
## [10,] 0.3459819 0 0 0 0 0
## [11,] -0.1891716 0 0 0 0 0
## [12,] -0.1366822 0 0 0 0 0
##
## $beta.se
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
## [2,] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
## [3,] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
## [4,] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
## [5,] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
## [6,] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
## [7,] 0.09187786 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
## [8,] 0.00000000 0.08313877 0.00000000 0.08281535 0.00000000 0.00000000
## [9,] 0.00000000 0.00000000 0.08536307 0.00000000 0.00000000 0.00000000
## [10,] 0.00000000 0.00000000 0.00000000 0.08791592 0.00000000 0.00000000
## [11,] 0.00000000 0.00000000 0.00000000 0.00000000 0.08416815 0.00000000
## [12,] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.06710476
## [,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.07903760 0 0 0 0 0
## [9,] 0.08640385 0 0 0 0 0
## [10,] 0.08899676 0 0 0 0 0
## [11,] 0.08518044 0 0 0 0 0
## [12,] 0.06789531 0 0 0 0 0
##
## $psi
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.9994470 0.54618107 0.43583976 0.33952475 -0.18363303
## [2,] 0.5461811 0.99994875 0.23527134 0.34397834 -0.13072104
## [3,] 0.4358398 0.23527134 0.99943868 0.28759278 -0.18513958
## [4,] 0.3395248 0.34397834 0.28759278 0.99963899 -0.21035567
## [5,] -0.1836330 -0.13072104 -0.18513958 -0.21035567 0.99911214
## [6,] -0.1831462 -0.08030562 -0.09845788 -0.08011838 0.09187896
## [7,] 0.0000000 0.00000000 0.00000000 0.00000000 0.00000000
## [8,] 0.0000000 0.00000000 0.00000000 0.00000000 0.00000000
## [9,] 0.0000000 0.00000000 0.00000000 0.00000000 0.00000000
## [10,] 0.0000000 0.00000000 0.00000000 0.00000000 0.00000000
## [11,] 0.0000000 0.00000000 0.00000000 0.00000000 0.00000000
## [12,] 0.0000000 0.00000000 0.00000000 0.00000000 0.00000000
## [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] -0.18314617 0.000000 0.0000000 0.000000 0.0000000 0.0000000
## [2,] -0.08030562 0.000000 0.0000000 0.000000 0.0000000 0.0000000
## [3,] -0.09845788 0.000000 0.0000000 0.000000 0.0000000 0.0000000
## [4,] -0.08011838 0.000000 0.0000000 0.000000 0.0000000 0.0000000
## [5,] 0.09187896 0.000000 0.0000000 0.000000 0.0000000 0.0000000
## [6,] 0.99862530 0.000000 0.0000000 0.000000 0.0000000 0.0000000
## [7,] 0.00000000 0.928056 0.0000000 0.000000 0.0000000 0.0000000
## [8,] 0.00000000 0.000000 0.6637395 0.000000 0.0000000 0.0000000
## [9,] 0.00000000 0.000000 0.0000000 0.761495 0.0000000 0.0000000
## [10,] 0.00000000 0.000000 0.0000000 0.000000 0.8498551 0.0000000
## [11,] 0.00000000 0.000000 0.0000000 0.000000 0.0000000 0.7784112
## [12,] 0.00000000 0.000000 0.0000000 0.000000 0.0000000 0.0000000
## [,12]
## [1,] 0.0000000
## [2,] 0.0000000
## [3,] 0.0000000
## [4,] 0.0000000
## [5,] 0.0000000
## [6,] 0.0000000
## [7,] 0.0000000
## [8,] 0.0000000
## [9,] 0.0000000
## [10,] 0.0000000
## [11,] 0.0000000
## [12,] 0.4926718
##
## $chisq
## chisq
## 37.22171
##
## $cfi
## cfi
## 1
##
## $tli
## tli
## 1.066706
##
## $rmsea
## rmsea
## 0
##
## $srmr
## srmr
## 0.04690272We 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),
CFI > .95
TLI > .95
RMSEA < .08
SRMR < .08
Checking if 3 of 4 rule 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] TRUEBecasue 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)
## NULLNodes 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 threshod for what is considered close-enough to equilibrium (threshhold), and can evalaute 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.455The 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)
## NULLInterpretation 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 mulitple 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 (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
## id 1 4200 310.36 77.12 321 308.40 121.57 203 439 236
## day 2 4200 10.57 5.99 10 10.50 7.41 1 21 20
## interaction 3 4200 73.88 45.18 71 72.32 54.86 1 168 167
## igaff 4 4190 7.39 1.78 8 7.69 1.48 1 9 8
## igdom 5 4190 6.92 1.62 7 7.08 1.48 1 9 8
## agval 6 4189 6.80 2.09 7 7.08 2.22 1 9 8
## agarous 7 4189 5.92 2.11 6 6.04 2.97 1 9 8
## stress 8 4191 1.22 1.38 1 1.02 1.48 0 5 5
## health 9 4192 1.08 1.17 1 0.93 1.48 0 5 5
## skew kurtosis se
## id 0.07 -1.21 1.19
## day 0.09 -1.16 0.09
## interaction 0.23 -1.00 0.70
## igaff -1.43 2.01 0.03
## igdom -0.93 0.91 0.02
## agval -0.93 0.13 0.03
## agarous -0.51 -0.56 0.03
## stress 0.97 0.02 0.02
## health 0.88 0.00 0.02From 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
ggplot(data = AMIB_interactionP2,
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)) +
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 nightime (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%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.568 9856.788
## agval 8805.920 9059.140
## agarous 8844.009 9097.229
## stress 6603.004 6856.225
## health 5759.347 6012.568
##
##
## 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.419 0.060
## igaff agval 1 -0.009 0.020 0.658 0.043
## igaff agarous 1 -0.039 0.018 0.032 0.026
## igaff stress 1 0.017 0.020 0.396 0.077
## 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.193 0.096
## igdom agval 1 -0.036 0.016 0.025 0.017
## igdom agarous 1 0.035 0.017 0.040 0.035
## igdom stress 1 0.026 0.012 0.032 0.017
## igdom health 1 0.006 0.012 0.597 0.028
## agval igaff 1 -0.006 0.025 0.815 0.060
## agval igdom 1 0.015 0.029 0.606 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.030 0.017 0.085 0.042
## agval health 1 0.004 0.014 0.766 0.022
## 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.203 0.049
## agarous health 1 -0.004 0.014 0.769 0.044
## stress igaff 1 -0.067 0.023 0.004 0.023
## stress igdom 1 0.003 0.026 0.903 0.066
## stress agval 1 -0.174 0.032 0.000 0.124
## stress agarous 1 0.053 0.041 0.205 0.185
## stress stress 1 0.641 0.042 0.000 0.208
## stress health 1 0.048 0.018 0.008 0.061
## health igaff 1 -0.010 0.031 0.736 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.031 0.251 0.044
## agval igaff 0.000 0.000 0.385 0.017 0.461 0.024
## agval igdom 0.011 0.086 0.058 0.006 0.169 0.019
## agarous igaff 0.207 0.073 0.041 0.006 0.103 0.010
## agarous igdom 0.000 0.000 0.225 0.005 0.246 0.009
## agarous agval 0.457 0.909 0.016 0.017 0.064 0.023
## stress igaff 0.002 0.037 -0.060 0.003 -0.240 0.021
## stress igdom 0.547 0.616 0.012 0.002 -0.060 0.008
## stress agval 0.000 0.000 -0.343 0.023 -0.426 0.031
## stress agarous 0.001 0.004 0.076 0.004 0.024 0.008
## health igaff 0.497 0.617 0.010 0.000 -0.100 0.014
## health igdom 0.587 0.626 -0.010 0.001 -0.064 0.003
## health agval 0.002 0.000 -0.090 0.005 -0.195 0.022
## health agarous 0.000 0.000 -0.105 0.004 -0.108 0.003
## health stress 0.000 0.000 0.188 0.043 0.244 0.049
##
##
## Between-subject effects:
## v1 v2 P 1->2 P 1<-2 pcor cor
## igdom igaff 0.053 0.011 0.320 0.422
## agval igaff 0.000 0.000 0.657 0.799
## agval igdom 0.202 0.332 -0.182 0.258
## agarous igaff 0.533 0.436 -0.027 0.328
## agarous igdom 0.000 0.013 0.472 0.533
## agarous agval 0.267 0.263 0.170 0.297
## stress igaff 0.510 0.493 0.101 -0.487
## stress igdom 0.583 0.778 0.071 -0.110
## stress agval 0.251 0.102 -0.219 -0.614
## stress agarous 0.740 0.780 0.050 -0.103
## health igaff 0.013 0.845 -0.173 -0.617
## health igdom 0.698 0.492 -0.090 -0.237
## health agval 0.361 0.252 -0.161 -0.681
## health agarous 0.891 0.649 0.023 -0.186
## health stress 0.000 0.000 0.641 0.782We 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 toegher (green positive edge), while “stress”" and “agval”" tend to fluctuate in oposition 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 multievel 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!