We created an R package named “pompom” (person-oriented modeling and perturbation on the model), and we will use the functions in “pompom” to compute iRAM (impulse response analysis metric) in this tutorial.
iRAM is built upon a hybrid method that combines intraindividual variability methods and network analysis methods in order to model individuals as high-dimensional dynamic systems. This hybrid method is designed and tested to quantify the extent of interaction in a high-dimensional multivariate system, and applicable on experience sampling data.
Note: please install the package “pompom”, before running the code in this tutorial.
This turorial corresponds to the following paper (under review): Yang et al. Impulse Response Analysis Metric (iRAM): Merging Intraindividual Variability, Network Analysis and Experience Sampling
library(pompom)
library(ggplot2)
library(qgraph)
require(reshape2)
set.seed(1234)
The 3-node network is a hypothetical example used for the purpose of illustration. As explained in the paper, the variables in time-series data are nodes and the temporal relations are edges, in the network terminology.
200 # number of observation
n.obs <- 3 # number of variables
p <- 0.1
noise.level <- cbind(rnorm(n.obs,0, noise.level),
epsilon <-rnorm(n.obs,0, noise.level),
rnorm(n.obs,0, noise.level)) # 3 time-series of noise for 3 variables
matrix(c(0,0,0,0,0,0,
true.beta <-0,0,0,0,0,0,
0,0,0,0,0,0,
0.2,0,0.25,0,0,0.6,
0,0.3,0,-0.2,0,-0.6,
0,-0.2,0.3,0,0,0),
nrow = 2 * p,
ncol = 2 * p,
byrow = TRUE)
matrix(true.beta[(p+1):(2*p),(p+1):(2*p)], nrow = p, ncol = p, byrow = F)
contemporaneous.relations <-1.relations <- matrix(true.beta[(p+1):(2*p),1:p], nrow = p, ncol = p, byrow = F) lag.
True model is:
true.beta
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.0 0.0 0.00 0.0 0 0.0
## [2,] 0.0 0.0 0.00 0.0 0 0.0
## [3,] 0.0 0.0 0.00 0.0 0 0.0
## [4,] 0.2 0.0 0.25 0.0 0 0.6
## [5,] 0.0 0.3 0.00 -0.2 0 -0.6
## [6,] 0.0 -0.2 0.30 0.0 0 0.0
To introduce the terminologies of networks, each of the 3 variables is depicted as a node (circle), and the temporal relations among the nodes are depicted as edges (arrows, and the arrow indicate directionality of the temporal relationship), with color indicating sign of relation (green = positive, red = negative), thickness indicating strength of relation, and line shape indicating temporality of the association (dashed = lag-1, solid = contemporaneous). Lag-1 relations mean the temporal relations between variables from measurement t -1 to the measurement t, and contemporaneous relations means the temporal relations between variables within the same measurement.
plot_network_graph(true.beta, p)
## NULL
Time series was simulated based on the temporal relations and process noise. Our hypohetical 3-node newtork was using simulated data based on a pre-defined temporal relationship matrix and process noise. The temporal relationship is as follows, and process noise is at Mean = 0, SD = .1.
\(\eta(t) = (I-\ A \ )^{-1}\Phi\eta(t-1) + (I-\ A \ )^{-1}\epsilon(t)\)
where \(\ A\) is the block of contemporaneous relations (southeast block of the beta matrix),
\(\Phi\) is the block of lag-1 relations (southwest block of the beta matrix),
\(\epsilon\) is the process noise.
We simulated 200 occasions to represent 200 repeated measurements of the three variables in the real studies.
matrix(rep(0, p * n.obs), nrow = n.obs, ncol = p)
time.series <-1,] <- rnorm(p,0,1)
time.series[
p
row <-for (row in 2:n.obs)
{ solve(diag(1,p, p) - contemporaneous.relations) %*%
time.series[row,] <- lag.1.relations %*% time.series[row-1,] +
solve(diag(1,p, p) - contemporaneous.relations) %*% epsilon[row, ]
} data.frame(time.series)
time.series <-names(time.series) <- c("x", "y", "z")
$time <- seq(1,length(time.series[,1]),1)
time.series
melt(time.series, id="time") ## convert to long format
time.series.long <-
ggplot(data=time.series.long,
aes(x=time, y=value, colour=variable)) +
geom_line() +
facet_wrap( ~ variable , ncol = 1) +
scale_y_continuous(breaks = seq(0, 100, by = 50)) +
theme(
strip.background = element_blank(),
panel.background = element_blank(),
legend.title = element_blank(),
# strip.text.x = element_blank(),
# strip.text.y = element_blank(),
# axis.text.y = element_text(),
# axis.title.y = element_blank()
legend.key = element_blank(),
legend.position = "none",
# legend.title = element_blank(),
# panel.background = element_blank(),
axis.text.y=element_text(color="black",size=12),
axis.text.x=element_text(color="black",size=12),
axis.title.y=element_text(color="black",size=12),
axis.title.x=element_text(color="black",size=12),
axis.line = element_line(color = 'black')
+
) ylim(-1,1)+
xlab("Time") +
ylab("Value")
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Warning: Removed 1 row(s) containing missing values (geom_path).
$time <- NULL time.series
The model fit summary will give and estimates, which are essential information to conduct network analysis later. The model fit summary also gives the fit statistics.
Parameters of the “uSEM” function:
var.number: number of variables, 13 in this case.
time.series: multivariate time-series data in the long format (every row is a measurement, and every column is a variable). Note: scaling is not absolutely necessary, but it can be helpful when some variables have small variances.
lag.order: number of lags in the uSEM model. Default is 1.
verbose: default at FALSE. If TRUE, it will print the top five largest modification indices in the lavaan model of each iteration.
trim: default at TRUE. If TRUE, it will trim the final model, eliminating all the insignificant temporal relations.
p # number of variables
var.number <- 1 # lag order of the model
lag.order <-
uSEM(var.number,
model.fit <-
time.series,
lag.order, # verbose = FALSE, #published code
verbose = TRUE, #test
trim = TRUE)
## [1] "Iteration: 1"
## [1] "Modification indices of beta: "
## lhs op rhs mi epc sepc.lv sepc.all sepc.nox lhs.number
## 115 eta4 ~ eta6 91.546 0.762 0.596 0.596 0.596 4
## 114 eta4 ~ eta5 70.369 -0.539 -0.558 -0.558 -0.558 4
## 113 eta4 ~ eta3 66.586 0.512 0.540 0.540 0.540 4
## 119 eta5 ~ eta6 58.914 -0.565 -0.428 -0.428 -0.428 5
## 123 eta6 ~ eta5 54.353 -0.361 -0.477 -0.477 -0.477 6
## [1] "largest improvement: eta4 ~ eta6"
## [1] "Iteration: 2"
## [1] "Modification indices of beta: "
## lhs op rhs mi epc sepc.lv sepc.all sepc.nox lhs.number
## 119 eta5 ~ eta6 58.914 -0.565 -0.428 -0.428 -0.428 5
## 123 eta6 ~ eta5 54.353 -0.361 -0.477 -0.477 -0.477 6
## 118 eta5 ~ eta4 21.062 -0.273 -0.264 -0.264 -0.264 5
## 114 eta4 ~ eta3 20.100 0.252 0.265 0.265 0.265 4
## 117 eta5 ~ eta3 13.928 -0.257 -0.262 -0.262 -0.262 5
## [1] "largest improvement: eta5 ~ eta6"
## [1] "Iteration: 3"
## [1] "Modification indices of beta: "
## lhs op rhs mi epc sepc.lv sepc.all sepc.nox lhs.number
## 115 eta4 ~ eta3 20.100 0.252 0.265 0.265 0.265 4
## 122 eta6 ~ eta4 10.958 -0.294 -0.377 -0.377 -0.377 6
## 121 eta6 ~ eta2 10.781 -0.189 -0.258 -0.258 -0.258 6
## 114 eta4 ~ eta2 6.883 -0.146 -0.155 -0.155 -0.155 4
## 116 eta4 ~ eta5 6.434 -0.166 -0.167 -0.167 -0.167 4
## [1] "largest improvement: eta4 ~ eta3"
## [1] "Iteration: 4"
## [1] "Modification indices of beta: "
## lhs op rhs mi epc sepc.lv sepc.all sepc.nox lhs.number
## 121 eta6 ~ eta2 10.781 -0.189 -0.258 -0.258 -0.258 6
## 123 eta6 ~ eta5 3.241 -0.188 -0.241 -0.241 -0.241 6
## 116 eta4 ~ eta5 2.478 -0.101 -0.101 -0.101 -0.101 4
## 119 eta5 ~ eta4 1.887 -0.087 -0.087 -0.087 -0.087 5
## 115 eta4 ~ eta2 1.009 -0.063 -0.067 -0.067 -0.067 4
## [1] "largest improvement: eta6 ~ eta2"
## [1] "Iteration: 5"
## [1] "Modification indices of beta: "
## lhs op rhs mi epc sepc.lv sepc.all sepc.nox lhs.number
## 122 eta6 ~ eta4 5.164 -0.357 -0.462 -0.462 -0.462 6
## 121 eta6 ~ eta1 4.110 -0.109 -0.140 -0.140 -0.140 6
## 117 eta4 ~ eta5 2.496 -0.102 -0.104 -0.104 -0.104 4
## 120 eta5 ~ eta4 1.872 -0.087 -0.085 -0.085 -0.085 5
## 116 eta4 ~ eta2 1.053 -0.066 -0.069 -0.069 -0.069 4
## [1] "largest improvement: eta6 ~ eta1"
## [1] "Iteration: 6"
## [1] "Modification indices of beta: "
## lhs op rhs mi epc sepc.lv sepc.all sepc.nox lhs.number
## 118 eta4 ~ eta5 2.510 -0.102 -0.106 -0.106 -0.106 4
## 121 eta5 ~ eta4 1.858 -0.086 -0.083 -0.083 -0.083 5
## 122 eta6 ~ eta4 1.089 -0.276 -0.353 -0.353 -0.353 6
## 117 eta4 ~ eta2 1.089 -0.068 -0.072 -0.072 -0.072 4
## 120 eta5 ~ eta3 0.779 -0.053 -0.054 -0.054 -0.054 5
## [1] "largest improvement: ~ "
Now the uSEM model result is in the object “model.fit”, including beta matrix, psi matrix, and fit statistics. Next, we need to parse the model.fit object into beta matrix and plot the estimated network graph.
parse_beta(var.number = p,
beta.matrix <-model.fit = model.fit,
lag.order = 1,
matrix = TRUE) # parse temporal relations in matrix format
plot_network_graph(beta.matrix$est,
var.number)
## NULL
estimated beta0 =
$est beta.matrix
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.0000000 0.0000000 0.0000000 0 0 0.0000000
## [2,] 0.0000000 0.0000000 0.0000000 0 0 0.0000000
## [3,] 0.0000000 0.0000000 0.0000000 0 0 0.0000000
## [4,] 0.3066036 0.0000000 0.2527896 0 0 0.5761694
## [5,] 0.0000000 0.4078627 0.0000000 0 0 -0.6707407
## [6,] -0.1091465 -0.2513549 0.3145093 0 0 0.0000000
$se beta.matrix
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.00000000 0.00000000 0.00000000 0 0 0.00000000
## [2,] 0.00000000 0.00000000 0.00000000 0 0 0.00000000
## [3,] 0.00000000 0.00000000 0.00000000 0 0 0.00000000
## [4,] 0.04621379 0.00000000 0.05333449 0 0 0.06745900
## [5,] 0.00000000 0.04745960 0.00000000 0 0 0.06460644
## [6,] 0.05327602 0.06321407 0.05600953 0 0 0.00000000
The “model_summary” function will return an object that contains beta, psi and fit statistics, shown in the following code chunks.
model_summary(model.fit,
mdl <-
var.number,
lag.order)
$beta mdl
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.0000000 0.0000000 0.0000000 0 0 0.0000000
## [2,] 0.0000000 0.0000000 0.0000000 0 0 0.0000000
## [3,] 0.0000000 0.0000000 0.0000000 0 0 0.0000000
## [4,] 0.3066036 0.0000000 0.2527896 0 0 0.5761694
## [5,] 0.0000000 0.4078627 0.0000000 0 0 -0.6707407
## [6,] -0.1091465 -0.2513549 0.3145093 0 0 0.0000000
$beta.se mdl
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.00000000 0.00000000 0.00000000 0 0 0.00000000
## [2,] 0.00000000 0.00000000 0.00000000 0 0 0.00000000
## [3,] 0.00000000 0.00000000 0.00000000 0 0 0.00000000
## [4,] 0.04621379 0.00000000 0.05333449 0 0 0.06745900
## [5,] 0.00000000 0.04745960 0.00000000 0 0 0.06460644
## [6,] 0.05327602 0.06321407 0.05600953 0 0 0.00000000
$psi mdl
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.02856458 -0.01829867 0.01245015 0.00000000 0.00000000 0.00000000
## [2,] -0.01829867 0.03239941 -0.02224377 0.00000000 0.00000000 0.00000000
## [3,] 0.01245015 -0.02224377 0.03176525 0.00000000 0.00000000 0.00000000
## [4,] 0.00000000 0.00000000 0.00000000 0.01006503 0.00000000 0.00000000
## [5,] 0.00000000 0.00000000 0.00000000 0.00000000 0.01009431 0.00000000
## [6,] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.01029627
Model fit
Model fit should obey a 3 out of 4 rule (CFI and TLI should be at least .95, and RMSEA and SRMR should be no greater than .08).
$cfi mdl
## cfi
## 0.9828771
$tli mdl
## tli
## 0.9743157
$rmsea mdl
## rmsea
## 0.07733488
$srmr mdl
## srmr
## 0.1207021
Conceptually, impulse response analysis is to perturb the system (the whole psychological network we estimated) by one node (or multiple nodes, but since this model is linear, the multiple nodes scenario is additive of the single node impulse. Thus one can conduct impulse response analysis in either way and get equivalent result). After the system receives such perturbation, or impulse, the impulse will reverbate through the network based on the statistical relationship.
For example, one estimated sadness can be predicted by -0.35 happiness + 0.2 anger - 0.6 self-esteem, if happiness or anger/self-esteem received a perturbation, then one can deduct what is the value of sadness for the next step, and we can also deduct value of other nodes in the same way. So one can have a time profile per node after the perturbation. Time profile is the trajectory of a variable after the system received the perturbation.
This method will give an intuitive view of the dynamic behavior of a network, in complimentary of the static newtork configuration (e.g. density, centrality, clustering, etc).
Since uSEM estimated contemporaneous relations as well, we transformed the contemporaneous relations back to a traditional lagged format, so that one can conduct impulse response analysis. Details of mathematical deduction is shown in Gates et al. (2010).
100 # number of steps to generate time profile
steps <- 200 # number of repilcations in bootstrap
replication <- .01 # setting threshold for approximate asymptote (iRAM calculation) threshold <-
# ponit estimate of iRAM
iRAM(model.fit,
point.estimate.iRAM <-beta = NULL,
var.number = var.number,
lag.order = lag.order,
threshold = threshold,
boot = FALSE,
replication = replication,
steps= steps)
$recovery.time point.estimate.iRAM
## [,1] [,2] [,3]
## [1,] 9 9 7
## [2,] 11 11 10
## [3,] 11 11 9
# point.estimate.iRAM$time.series.data
plot_time_profile(point.estimate.iRAM$time.series.data,
var.number = 3,
threshold = threshold,
xupper = 50)
## NULL
iRAM(model.fit,
bootstrap.iRAM <-beta = NULL,
var.number = var.number,
lag.order = lag.order,
threshold = threshold,
boot = TRUE,
replication = replication,
steps= steps
)
$mean bootstrap.iRAM
## [,1] [,2] [,3]
## [1,] 8.235 8.340 6.965
## [2,] 11.400 11.575 9.990
## [3,] 10.900 11.085 9.460
$upper bootstrap.iRAM
## [,1] [,2] [,3]
## [1,] 12 13 11
## [2,] 17 18 17
## [3,] 18 19 16
$lower bootstrap.iRAM
## [,1] [,2] [,3]
## [1,] 4 5 4
## [2,] 8 7 6
## [3,] 6 7 5
plot_time_profile(bootstrap.iRAM$time.profile.data,
var.number = 3,
threshold = threshold,
xupper = 25)
## NULL
plot_iRAM_dist(bootstrap.iRAM$recovery.time.reps)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## NULL
iRAM(
true.iRAM <-model.fit= NULL,
beta = true.beta,
var.number = var.number,
lag.order = lag.order,
threshold = threshold,
boot = FALSE,
replication = replication,
steps= steps)
plot_time_profile(true.iRAM$time.series.data,
var.number = 3,
threshold = threshold,
xupper = 25)
## NULL
0
sum.diff <-for (row in 1:nrow(bootstrap.iRAM$recovery.time))
{ sum.diff + (bootstrap.iRAM$recovery.time.reps[row,] -
sum.diff <- c(true.iRAM$recovery.time[1,],
$recovery.time[2,],
true.iRAM$recovery.time[3,]))^2
true.iRAM
} sqrt(sum.diff/nrow(bootstrap.iRAM$recovery.time))
RMSE <-
0
sum.diff <-for (row in 1:nrow(bootstrap.iRAM$recovery.time))
{ sum.diff + (bootstrap.iRAM$recovery.time.reps[row,] -
sum.diff <- c(true.iRAM$recovery.time[1,],
$recovery.time[2,],
true.iRAM$recovery.time[3,]))/
true.iRAM (c(true.iRAM$recovery.time[1,],
$recovery.time[2,],
true.iRAM$recovery.time[3,]))
true.iRAM
} 1/nrow(bootstrap.iRAM$recovery.time) * sum.diff
RB <-
NULL
SD <-for (col in 1:(var.number^2))
{ cbind(SD, sd(bootstrap.iRAM$recovery.time.reps[,col]))
SD <- }
true.iRAM: iRAM calculated based on true beta matrix.
boostrap.iRAM$mean: the iRAM calcualted across bootstrapped replications.
RMSE: root mean squared error of the estimated iRAM compared with true iRAM across bootstrapped replications.
RB: relative bias, mean of summed difference between the estimated iRAM and the true iRAM across bootstrapped replications.
SD: standard deviation of the estimated iRAM across bootstrapped replications.
# metrics
$recovery.time # true value true.iRAM
## [,1] [,2] [,3]
## [1,] 8 8 7
## [2,] 10 11 9
## [3,] 11 12 11
$mean # estimated mean bootstrap.iRAM
## [,1] [,2] [,3]
## [1,] 8.235 8.340 6.965
## [2,] 11.400 11.575 9.990
## [3,] 10.900 11.085 9.460
RMSE
## [1] 1.719011 2.111871 1.629417 3.078961 3.007491 2.944486 3.091925 3.378609
## [9] 3.528456
RB
## [1] 0.029375000 0.042500000 -0.005000000 0.140000000 0.052272727
## [6] 0.110000000 -0.009090909 -0.076250000 -0.140000000
SD
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 1.707146 2.089553 1.633129 2.749143 2.95942 2.780026 3.098062 3.260511
## [,9]
## [1,] 3.182616
Sometimes the raw time-series data is non-stationary (e.g., the absolute value of auto-regression is close to or above 1), and one way to model the dynamics is to take the difference score of the time-series and then model the difference scores.
In this case, the temporal relations are describing how the difference score are predicting one another, so the above steps still applies. When conducting impulse response analysis, the equilibrium of the integrated form of time profile can also be computed and plotted.
iRAM_equilibrium(beta = mdl$beta,
iRAM_equilibrium_value <-var.number = var.number,
lag.order = lag.order)
iRAM_equilibrium_value
## e11 e12 e13 e21 e22 e23 e31
## 1 0.08795343 0.3356275 -0.2962957 -0.9552423 1.593893 -0.7990282 1.674173
## e32 e33
## 1 -1.461466 0.8823358
plot_integrated_time_profile(beta.matrix = mdl$beta,
var.number = 3,
lag.order = 1)
## NULL