| tvmgmsampler {mgm} | R Documentation |
Sample from time-varying k-order Mixed Graphical Model
Description
Generates samples from a time-varying k-order Mixed Graphical Model
Usage
tvmgmsampler(factors, interactions, thresholds, sds, type,
level, nIter = 250, pbar = TRUE, ...)
Arguments
factors |
The same object as |
interactions |
The same object as |
thresholds |
A list with p entries for p variables, each of which contains a N x m matrix. The columns contain the m thresholds for m categories (for continuous variables m = 1 and the entry contains the threshold/intercept). The rows indicate how the thresholds change over time. |
sds |
N x p matrix indicating the standard deviations of Gaussians specified in |
type |
p character vector indicating the type of variable for each column in |
level |
p integer vector indicating the number of categories of each variable. For continuous variables set to 1. |
nIter |
Number of iterations in the Gibbs sampler until a sample is drawn. |
pbar |
If |
... |
Additional arguments. |
Details
tvmgmsampler is a wrapper function around mgmsampler. Its input is the same as for mgmsampler, except that each object has an additional dimension for time. The number of time points is specified via entries in the additional time dimension.
Value
A list containing:
call |
Contains all provided input arguments. |
data |
The N x p data matrix of sampled values |
Author(s)
Jonas Haslbeck <jonashaslbeck@gmail.com>
References
Haslbeck, J. M. B., & Waldorp, L. J. (2020). mgm: Estimating time-varying Mixed Graphical Models in high-dimensional Data. Journal of Statistical Software, 93(8), pp. 1-46. DOI: 10.18637/jss.v093.i08
Yang, E., Baker, Y., Ravikumar, P., Allen, G. I., & Liu, Z. (2014, April). Mixed Graphical Models via Exponential Families. In AISTATS (Vol. 2012, pp. 1042-1050).
Examples
## 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(0, 1, length = 15),
bandwidth = .2,
k = 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(0, 1, length = 15),
bandwidth = .2,
k = 3,
lambdaSel = "CV",
ruleReg = "AND")
tvmgm_obj$interactions$indicator
# Seems very difficult to recover this time-varying 3-way binary interaction
# See also the corresponding problems in the examples of ?mgmsampler
# For more examples see https://github.com/jmbh/mgmDocumentation
## End(Not run)