Three-node Example to Compute iRAM (impulse response analysis metric)

Overview

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

Outline of Tutorial

Prepare simulation of time-series data
Apply Intraindividual Variability Method (unified Structural SEM, uSEM) on multivariate time-series data
Apply a network analysis method (also called impulse reponse analysis) on person-specific network
Calculate impulse response analysis metric (iRAM)

Step 1. Prepare simulation of time-series data

library(pompom)
library(ggplot2)
library(qgraph)
require(reshape2)
set.seed(1234)

Set up a 3 node model

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.

n.obs <- 200 # number of observation
p <- 3 # number of variables
noise.level <- 0.1
epsilon <- cbind(rnorm(n.obs,0, noise.level), 
             rnorm(n.obs,0, noise.level), 
             rnorm(n.obs,0, noise.level)) # 3 time-series of noise for 3 variables

true.beta <-  matrix(c(0,0,0,0,0,0,
                       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)

contemporaneous.relations <- matrix(true.beta[(p+1):(2*p),(p+1):(2*p)], nrow = p, ncol = p, byrow = F)
lag.1.relations <- matrix(true.beta[(p+1):(2*p),1:p], nrow = p, ncol = p, byrow = F)

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

Network graph

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

Simulate time-series

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.

η(t) = (I− A )⁻¹Φη(t − 1) + (I− A )⁻¹ϵ(t)

where A is the block of contemporaneous relations (southeast block of the beta matrix),

Φ is the block of lag-1 relations (southwest block of the beta matrix),

ϵ is the process noise.

We simulated 200 occasions to represent 200 repeated measurements of the three variables in the real studies.

time.series <- matrix(rep(0, p * n.obs), nrow = n.obs, ncol = p)
time.series[1,] <- rnorm(p,0,1)

row <- p
for (row in 2:n.obs)
{
  time.series[row,] <- solve(diag(1,p, p) - contemporaneous.relations) %*%
    lag.1.relations %*% time.series[row-1,] +
    solve(diag(1,p, p) - contemporaneous.relations) %*% epsilon[row, ]
}
time.series <- data.frame(time.series)
names(time.series) <- c("x", "y", "z")

Plot the simulated time-series

time.series$time <- seq(1,length(time.series[,1]),1)

time.series.long <- melt(time.series, id="time")  ## convert to long format

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 containing missing values or values outside the scale range
## (`geom_line()`).

time.series$time <- NULL

Step 2. Apply Intraindividual Variability Method (unified Structural SEM, uSEM) on multivariate time-series data

Use uSEM to fit the model and parse parameters.

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.

var.number <- p # number of variables
lag.order <- 1 # lag order of the model

model.fit <- uSEM(var.number, 
                  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 90.001  0.761   0.597    0.597    0.597          4
## 114 eta4  ~ eta5 67.805 -0.536  -0.549   -0.549   -0.549          4
## 113 eta4  ~ eta3 64.833  0.505   0.535    0.535    0.535          4
## 119 eta5  ~ eta6 59.407 -0.571  -0.438   -0.438   -0.438          5
## 123 eta6  ~ eta5 53.301 -0.362  -0.472   -0.472   -0.472          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 59.407 -0.571  -0.438   -0.438   -0.438          5
## 123 eta6  ~ eta5 53.301 -0.362  -0.472   -0.472   -0.472          6
## 118 eta5  ~ eta4 20.861 -0.272  -0.267   -0.267   -0.267          5
## 114 eta4  ~ eta3 19.872  0.250   0.263    0.263    0.263          4
## 117 eta5  ~ eta3 13.889 -0.255  -0.264   -0.264   -0.264          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 19.872  0.250   0.263    0.263    0.263          4
## 122 eta6  ~ eta4 11.857 -0.311  -0.398   -0.398   -0.398          6
## 121 eta6  ~ eta2  9.340 -0.178  -0.239   -0.239   -0.239          6
## 114 eta4  ~ eta2  6.805 -0.147  -0.154   -0.154   -0.154          4
## 116 eta4  ~ eta5  6.295 -0.166  -0.165   -0.165   -0.165          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 9.340 -0.178  -0.239   -0.239   -0.239          6
## 123 eta6  ~ eta5 2.553 -0.171  -0.218   -0.218   -0.218          6
## 116 eta4  ~ eta5 2.412 -0.101  -0.100   -0.100   -0.100          4
## 119 eta5  ~ eta4 1.810 -0.085  -0.086   -0.086   -0.086          5
## 115 eta4  ~ eta2 0.963 -0.062  -0.065   -0.065   -0.065          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.382 -0.371  -0.478   -0.478   -0.478          6
## 121 eta6  ~ eta1 4.382 -0.113  -0.144   -0.144   -0.144          6
## 117 eta4  ~ eta5 2.427 -0.101  -0.102   -0.102   -0.102          4
## 120 eta5  ~ eta4 1.797 -0.085  -0.084   -0.084   -0.084          5
## 116 eta4  ~ eta2 1.000 -0.064  -0.067   -0.067   -0.067          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.439 -0.102  -0.104   -0.104   -0.104          4
## 121 eta5  ~ eta4 1.783 -0.084  -0.082   -0.082   -0.082          5
## 117 eta4  ~ eta2 1.033 -0.066  -0.070   -0.070   -0.070          4
## 122 eta6  ~ eta4 1.033 -0.280  -0.357   -0.357   -0.357          6
## 120 eta5  ~ eta3 0.804 -0.053  -0.055   -0.055   -0.055          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.

beta.matrix <- parse_beta(var.number = p, 
                          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 =

beta.matrix$est

##            [,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.3046405  0.0000000 0.2521295    0    0  0.5737546
## [5,]  0.0000000  0.3979234 0.0000000    0    0 -0.6654005
## [6,] -0.1125467 -0.2409437 0.3133049    0    0  0.0000000

beta.matrix$se

##            [,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.04665196 0.00000000 0.05336633    0    0 0.06790949
## [5,] 0.00000000 0.04818624 0.00000000    0    0 0.06459937
## [6,] 0.05317012 0.06361336 0.05582500    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.

mdl <- model_summary(model.fit, 
                     var.number, 
                     lag.order)

mdl$beta

##            [,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.3046405  0.0000000 0.2521295    0    0  0.5737546
## [5,]  0.0000000  0.3979234 0.0000000    0    0 -0.6654005
## [6,] -0.1125467 -0.2409437 0.3133049    0    0  0.0000000

mdl$beta.se

##            [,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.04665196 0.00000000 0.05336633    0    0 0.06790949
## [5,] 0.00000000 0.04818624 0.00000000    0    0 0.06459937
## [6,] 0.05317012 0.06361336 0.05582500    0    0 0.00000000

mdl$psi

##             [,1]        [,2]        [,3]       [,4]       [,5]      [,6]
## [1,]  0.02737529 -0.01662328  0.01122336 0.00000000 0.00000000 0.0000000
## [2,] -0.01662328  0.03003924 -0.02051554 0.00000000 0.00000000 0.0000000
## [3,]  0.01122336 -0.02051554  0.03049976 0.00000000 0.00000000 0.0000000
## [4,]  0.00000000  0.00000000  0.00000000 0.01006035 0.00000000 0.0000000
## [5,]  0.00000000  0.00000000  0.00000000 0.00000000 0.01003422 0.0000000
## [6,]  0.00000000  0.00000000  0.00000000 0.00000000 0.00000000 0.0102251

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).

mdl$cfi

## cfi 
##   1

mdl$tli

##      tli 
## 1.002499

mdl$rmsea

## rmsea 
##     0

mdl$srmr

##       srmr 
## 0.01203164

Step 3. Apply a network analysis method (also called impulse reponse analysis) on person-specific network

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).

Set parameters for impulse response analysis

steps <- 100 # number of steps to generate time profile 
replication <- 200 # number of repilcations in bootstrap 
threshold <- .01 # setting threshold for approximate asymptote (iRAM calculation)

Step 4. Calculate impulse response analysis metric (iRAM)

Point estimation of iRAM

# ponit estimate of iRAM 
point.estimate.iRAM <- iRAM(model.fit, 
                            beta = NULL, 
                            var.number = var.number, 
                            lag.order = lag.order, 
                            threshold = threshold,
                            boot = FALSE, 
                            replication = replication,
                            steps= steps)

point.estimate.iRAM$recovery.time

##      [,1] [,2] [,3]
## [1,]    8    8    7
## [2,]   10   10    9
## [3,]   10   10    8

# point.estimate.iRAM$time.series.data

Plot time profile

plot_time_profile(point.estimate.iRAM$time.series.data, 
                  var.number = 3,
                  threshold = threshold, 
                  xupper = 50)

## NULL

Bootstraped iRAM

bootstrap.iRAM <- iRAM(model.fit, 
                       beta = NULL, 
                       var.number = var.number, 
                       lag.order = lag.order,
                       threshold = threshold,
                       boot = TRUE, 
                       replication = replication,
                       steps= steps
                       )

bootstrap.iRAM$mean

##        [,1]   [,2]  [,3]
## [1,]  7.965  8.110 6.690
## [2,] 10.695 10.805 9.305
## [3,] 10.140 10.265 8.730

bootstrap.iRAM$upper

##      [,1] [,2] [,3]
## [1,]   11   12   10
## [2,]   16   17   15
## [3,]   16   17   15

bootstrap.iRAM$lower

##      [,1] [,2] [,3]
## [1,]    5    6    4
## [2,]    7    7    6
## [3,]    6    6    5

Plot boostrapped time profiles

plot_time_profile(bootstrap.iRAM$time.profile.data, 
                  var.number = 3,
                  threshold = threshold, 
                  xupper = 25)

## NULL

Plot boostrapped iRAM distribution

plot_iRAM_dist(bootstrap.iRAM$recovery.time.reps)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## NULL

Check accuracy of iRAM estimation

true.iRAM <- 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

sum.diff <- 0
for (row in 1:nrow(bootstrap.iRAM$recovery.time))
{
  sum.diff <- sum.diff + (bootstrap.iRAM$recovery.time.reps[row,] -
    c(true.iRAM$recovery.time[1,], 
      true.iRAM$recovery.time[2,], 
      true.iRAM$recovery.time[3,]))^2
}
RMSE <- sqrt(sum.diff/nrow(bootstrap.iRAM$recovery.time))


sum.diff <- 0
for (row in 1:nrow(bootstrap.iRAM$recovery.time))
{
  sum.diff <- sum.diff + (bootstrap.iRAM$recovery.time.reps[row,] -
    c(true.iRAM$recovery.time[1,], 
      true.iRAM$recovery.time[2,], 
      true.iRAM$recovery.time[3,]))/
    (c(true.iRAM$recovery.time[1,], 
      true.iRAM$recovery.time[2,], 
      true.iRAM$recovery.time[3,]))
}
RB <- 1/nrow(bootstrap.iRAM$recovery.time) * sum.diff

SD <- NULL
for (col in 1:(var.number^2))
{
  SD <- cbind(SD, sd(bootstrap.iRAM$recovery.time.reps[,col]))
}

Summary statistics

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
true.iRAM$recovery.time # true value

##      [,1] [,2] [,3]
## [1,]    8    8    7
## [2,]   10   11    9
## [3,]   11   12   11

bootstrap.iRAM$mean # estimated mean

##        [,1]   [,2]  [,3]
## [1,]  7.965  8.110 6.690
## [2,] 10.695 10.805 9.305
## [3,] 10.140 10.265 8.730

RMSE

## [1] 1.528071 1.788854 1.452584 2.486966 2.566125 2.436186 2.849561 3.394849
## [9] 3.622154

RB

## [1] -0.00437500  0.01375000 -0.04428571  0.06950000 -0.01772727  0.03388889
## [7] -0.07818182 -0.14458333 -0.20636364

SD

##          [,1]    [,2]     [,3]     [,4]     [,5]     [,6]     [,7]     [,8]
## [1,] 1.531503 1.78995 1.422681 2.393873 2.565127 2.423083 2.723506 2.925331
##          [,9]
## [1,] 2.829688

Integrated form of impulse respnose analysis

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.

Compute the equilibrium

iRAM_equilibrium_value <- iRAM_equilibrium(beta = mdl$beta,
                         var.number = var.number, 
                         lag.order = lag.order)

iRAM_equilibrium_value

##          e11       e12        e13        e21      e22       e23      e31
## 1 0.09116162 0.3228347 -0.2921117 -0.8578256 1.459135 -0.722253 1.569803
##         e32       e33
## 1 -1.329204 0.8047852

Plot the time profile in the integrated form

plot_integrated_time_profile(beta.matrix = mdl$beta, 
                             var.number = 3, 
                             lag.order = 1)

## NULL