March 15, 2017

Necessary Software

Install & Load Packages

# For Estimation
# install.packages('mgm')
library(mgm)

# For Visualization
# install.packages('qgraph') 
library(qgraph)

The Data

Subset of 16 continuous variables of larger ESM dataset:

  • One individual diagnosed with Major Depression
  • Up to 10 measurements a day
  • Over 238 consecutive days
  • In total 1476 Measurement points

Openly available: https://osf.io/j4fg8/

We select 16 mood related variables: Feeling relaxed, feeling down, feeling irritated, feeling lonely, feeling anxious, feeling satisfied, etc.

The Data

Load Data

data <- read.csv('http://jmbh.github.io/data/mood_data.csv', 
                 header = TRUE)

data[1:3, 1:8]
##       date beepno beeptime Relaxed Down Irritated Satisfied Lonely
## 1 13/08/12      1    08:58       5   -1         1         5     -1
## 2 14/08/12      5    14:32       4    0         3         3      0
## 3 14/08/12      6    16:17       4    0         2         3      0
data_mood <- data[, -c(1:3)] # delete id variables
dim(data_mood)
## [1] 1476   16
p <- ncol(data_mood)

Estimate Mixed Graphical Model

mgm_obj <- mgmfit(data = data_mood, 
               type = rep('g', p), 
               lev = rep(1, p),
               d = 1,
               lambda.sel = 'CV', 
               pbar = FALSE)
## 
##  Note that signs of edge weights (if defined) are stored in fitobject$signs. See ?mgmfit for more info.
# Subset of weighted adjacency matrix
round(mgm_obj$wadj, 2)[1:4, 1:4]
##      [,1] [,2] [,3] [,4]
## [1,] 0.00 0.10 0.25 0.13
## [2,] 0.10 0.00 0.00 0.11
## [3,] 0.25 0.00 0.00 0.06
## [4,] 0.13 0.11 0.06 0.00

Get Signs of Pairwise Interactions

# Add signs to weighted adjacency matrix
wadj <- mgm_obj$wadj * mgm_obj$signs
wadj[is.na(wadj)] <- 0

round(wadj, 2)[1:4, 1:4]
##       [,1]  [,2]  [,3]  [,4]
## [1,]  0.00  0.10 -0.25  0.13
## [2,]  0.10  0.00  0.00 -0.11
## [3,] -0.25  0.00  0.00 -0.06
## [4,]  0.13 -0.11 -0.06  0.00

Visualize as Network

Q <- qgraph(wadj, layout = 'spring', labels = colnames(data_mood))

Compute Prediction Error

pred_obj <- predict(mgm_obj, data_mood, error.continuous = 'VarExpl')

pred_obj$error[1:10,]
##        Variable Error ErrorType
## 1       Relaxed 0.605   VarExpl
## 2          Down 0.761   VarExpl
## 3     Irritated 0.712   VarExpl
## 4     Satisfied 0.731   VarExpl
## 5        Lonely 0.730   VarExpl
## 6       Anxious 0.617   VarExpl
## 7  Enthusiastic 0.674   VarExpl
## 8    Suspicious 0.485   VarExpl
## 9      Cheerful 0.712   VarExpl
## 10       Guilty 0.581   VarExpl

Visualize Network with Prediction Error

qgraph(wadj, layout = 'spring', labels = colnames(data_mood), 
       pie = pred_obj$error$Error, pieColor = 'darkgrey')

Estimate lagged (VAR) model

var_obj <- var.mgm(data = data_mood, 
               type = rep('g', p), 
               lev = rep(1, p),
               lags = 1,
               lambda.sel = 'CV', 
               pbar = FALSE)
## 
##  Note that signs of edge weights (if defined) are stored in fitobject$signs. See ?var.mgm for more info.
round(var_obj$wadj, 2)[1:4, 1:4]
##      [,1] [,2] [,3] [,4]
## [1,] 0.06 0.05    0    0
## [2,] 0.05 0.21    0    0
## [3,] 0.00 0.00    0    0
## [4,] 0.07 0.15    0    0

Get Signs of Lag-1 effects

# Add signs to lagged effects
evar <- var_obj$wadj * var_obj$signs
evar[is.na(evar)] <- 0

round(evar, 2)[1:4, 1:4]
##       [,1]  [,2] [,3] [,4]
## [1,]  0.06 -0.05    0    0
## [2,] -0.05  0.21    0    0
## [3,]  0.00  0.00    0    0
## [4,]  0.07 -0.15    0    0

Visualize VAR Network

qgraph(evar, layout = 'spring', labels = colnames(data_mood))

Assumption of Stationarity

Assumption of Stationarity?

Specification of Time Vector

# Calculate Time Vector
dates <- as.character(data$date)
times <- as.character(data$beeptime)
x <- paste(dates, times)
date_time <- strptime(x, "%d/%m/%y %H:%M")
timeFromZero <- date_time - date_time[1]
timeFromZero <- as.numeric(timeFromZero)
plot(timeFromZero, type = 'l') # how irregular?

Estimating time-varying MGM

tv_mgm <- tv.mgmfit(data = data_mood, 
                    type = rep('g', p), 
                    lev = rep(1, p), 
                    timepoints = timeFromZero,
                    tsteps = 15,
                    bandwidth = .1,
                    d = 1,
                    pbar = FALSE)
## 
##  Note that signs of edge weights (if defined) are stored in fitobject$signs. See ?tv.mgmfit for more info.
# Add signs
tvwadj <- tv_mgm$wadj * tv_mgm$sign
tvwadj[is.na(tvwadj)] <- 0

Visualize tv MGM at t = 1, 8, 15

cln <- colnames(data_mood)

par(mfrow=c(1,3))
qgraph(tvwadj[, , 1], layout = Q$layout, labels = cln, title = 'Time = 1')
qgraph(tvwadj[, , 8], layout = Q$layout, labels = cln, title = 'Time = 8')
qgraph(tvwadj[, , 15], layout = Q$layout, labels = cln, title = 'Time = 15')

Visualize tv MGM: Specific edges

plot(tvwadj[2, 5, ], type = 'l', ylim = c(-.2,.8),
     xlab = 'Time', ylab = 'Parameter value')
lines(tvwadj[5, 15, ], col = 'red')
legend('topleft', c('Down - Lonely', 'Lonely - Worrying'), 
       lwd=c(1, 1), col = c('black', 'red'))

Other features