if (FALSE) {
## We specify a time-varying MGM and recover it using tvmgm()
# 1) Specify Model
# a) Define Graph
p <- 6
type = c("c", "c", "g", "g", "p", "p")
level = c(2, 3, 1, 1, 1, 1)
n_timepoints <- 1000
# b) Define Interaction
factors <- list()
factors[[1]] <- matrix(c(1,2,
2,3,
3,4), ncol=2, byrow = T) # no pairwise interactions
factors[[2]] <- matrix(c(1,2,3,
2,3,4), ncol=3, byrow = T) # one 3-way interaction
interactions <- list()
interactions[[1]] <- vector("list", length = 3)
interactions[[2]] <- vector("list", length = 2)
# 3 2-way interactions
interactions[[1]][[1]] <- array(0, dim = c(level[1], level[2], n_timepoints))
interactions[[1]][[2]] <- array(0, dim = c(level[2], level[3], n_timepoints))
interactions[[1]][[3]] <- array(0, dim = c(level[3], level[4], n_timepoints))
# 2 3-way interactions
interactions[[2]][[1]] <- array(0, dim = c(level[1], level[2], level[3], n_timepoints))
interactions[[2]][[2]] <- array(0, dim = c(level[2], level[3], level[4], n_timepoints))
theta <- .3
interactions[[1]][[1]][1, 1, ] <- theta
interactions[[1]][[2]][1, 1, ] <- theta
interactions[[1]][[3]][1, 1, ] <- seq(0, theta, length = n_timepoints)
interactions[[2]][[1]][1, 1, 1, ] <- theta
interactions[[2]][[2]][1, 1, 1, ] <- theta
# c) Define Thresholds
thresholds <- list()
thresholds[[1]] <- matrix(0, nrow = n_timepoints, ncol= level[1])
thresholds[[2]] <- matrix(0, nrow = n_timepoints, ncol= level[2])
thresholds[[3]] <- matrix(0, nrow = n_timepoints, ncol= level[3])
thresholds[[4]] <- matrix(0, nrow = n_timepoints, ncol= level[4])
thresholds[[5]] <- matrix(.1, nrow = n_timepoints, ncol= level[5])
thresholds[[6]] <- matrix(.1, nrow = n_timepoints, ncol= level[6])
# d) define sds
sds <- matrix(.2, ncol=p, nrow=n_timepoints)
# 2) Sample Data
set.seed(1)
d_iter <- tvmgmsampler(factors = factors,
interactions = interactions,
thresholds = thresholds,
sds = sds,
type = type,
level = level,
nIter = 100,
pbar = TRUE)
data <- d_iter$data
head(data)
# delete inf rows:
ind_finite <- apply(data, 1, function(x) if(all(is.finite(x))) TRUE else FALSE)
table(ind_finite) # all fine for this setup & seed
# in case of inf values (no theory on how to keep k-order MGM well-defined)
data <- data[ind_finite, ]
# 3) Recover
mgm_c_cv <- tvmgm(data = data,
type = type,
level = level,
k = 3,
estpoints = seq(0, 1, length=10),
bandwidth = .1,
lambdaSel = "CV",
ruleReg = "AND",
pbar = TRUE,
overparameterize = T,
signInfo = FALSE)
# Look at time-varying pairwise parameter 3-4
mgm_c_cv$pairwise$wadj[3,4,] # recovers increase
# 4) Predict values / compute nodewise Errors
pred_mgm_cv_w <- predict.mgm(mgm_c_cv,
data = data,
tvMethod = "weighted")
pred_mgm_cv_cM <- predict.mgm(mgm_c_cv,
data = data,
tvMethod = "closestModel")
pred_mgm_cv_w$errors
pred_mgm_cv_cM$errors # Pretty similar!
# For more examples see https://github.com/jmbh/mgmDocumentation
}
Run the code above in your browser using DataLab