# NOT RUN {
# }
# NOT RUN {
# --------- Example 1: p = 4 dimensional Gaussian ---------
# ----- 1) Specify Model -----
# a) General Graph Info
type = c('g', 'g', 'g', 'g') # Four Gaussians
level = c(1, 1, 1, 1)
n_timepoints = 500 # Number of time points
# b) Define Interaction
factors <- list()
factors[[1]] <- array(NA, dim=c(2, 2)) # two pairwise interactions
factors[[1]][1, 1:2] <- c(3,4)
factors[[1]][2, 1:2] <- c(1,2)
# Two parameters, one linearly increasing from 0 to 0.8, another one lin decreasing from 0.8 to 0
interactions <- list()
interactions[[1]] <- vector('list', length = 2)
interactions[[1]][[1]] <- array(0, dim = c(level[1], level[2], n_timepoints))
interactions[[1]][[1]][1,1, ] <- seq(.8, 0, length = n_timepoints)
interactions[[1]][[2]] <- array(0, dim = c(level[1], level[2], n_timepoints))
interactions[[1]][[2]][1,1, ] <- seq(0, .8, length = n_timepoints)
# c) Define Thresholds
thresholds <- vector('list', length = 4)
thresholds <- lapply(thresholds, function(x) matrix(0, ncol = level[1], nrow = n_timepoints))
# d) Define Standard deviations
sds <- matrix(1, ncol = length(type), nrow = n_timepoints) # constant across variables and time
# ----- 2) Sample cases -----
set.seed(1)
dlist <- tvmgmsampler(factors = factors,
interactions = interactions,
thresholds = thresholds,
sds = sds,
type = type,
level = level,
nIter = 75,
pbar = TRUE)
# ----- 3) Recover model from sampled cases -----
set.seed(1)
tvmgm_obj <- tvmgm(data = dlist$data,
type = type,
level = level,
estpoints = seq(1, n_timepoints, length = 15),
bandwidth = .2,
d = 2,
lambdaSel = 'CV',
ruleReg = 'AND')
# How well did we recover those two time-varying parameters?
plot(tvmgm_obj$pairwise$wadj[3,4,], type='l', ylim=c(0,.8))
lines(tvmgm_obj$pairwise$wadj[1,2,], type='l', col='red')
# Looks quite good
# --------- Example 2: p = 5 binary; one 3-way interaction ---------
# ----- 1) Specify Model -----
# a) General Graph Info
p <- 5 # number of variables
type = rep('c', p) # all categorical
level = rep(2, p) # all binary
n_timepoints <- 1000
# b) Define Interaction
factors <- list()
factors[[1]] <- NULL # no pairwise interactions
factors[[2]] <- array(NA, dim = c(1,3)) # one 3-way interaction
factors[[2]][1, 1:3] <- c(1, 2, 3)
interactions <- list()
interactions[[1]] <- NULL # no pairwise interactions
interactions[[2]] <- vector('list', length = 1) # one 3-way interaction
# 3-way interaction no1
interactions[[2]][[1]] <- array(0, dim = c(level[1], level[2], level[3], n_timepoints))
theta <- 2
interactions[[2]][[1]][1, 1, 1, ] <- seq(0, 2, length = n_timepoints) # fill in nonzero entries
# c) Define Thresholds
thresholds <- list()
for(i in 1:p) thresholds[[i]] <- matrix(0, nrow = n_timepoints, ncol = level[i])
# ----- 2) Sample cases -----
set.seed(1)
dlist <- tvmgmsampler(factors = factors,
interactions = interactions,
thresholds = thresholds,
type = type,
level = level,
nIter = 150,
pbar = TRUE)
# ----- 3) Check Marginals -----
dat <- dlist$data[1:round(n_timepoints/2),]
table(dat[,1], dat[,2], dat[,3])
dat <- dlist$data[round(n_timepoints/2):n_timepoints,]
table(dat[,1], dat[,2], dat[,3])
# Observation: much stronger effect in second hald of the time-series,
# which is what we expect
# ----- 4) Recover model from sampled cases -----
set.seed(1)
tvmgm_obj <- tvmgm(data = dlist$data,
type = type,
level = level,
estpoints = seq(1, n_timepoints, length = 15),
bandwidth = .2,
d = 2,
lambdaSel = 'CV',
ruleReg = 'AND')
tvmgm_obj$rawfactor$indicator
# Seems very difficult to recover this time-varying 3-way binary interaction
# See also the corresponding stationary problems in the examples of ?mgmsampler
# }
# NOT RUN {
# }
Run the code above in your browser using DataCamp Workspace