if (FALSE) {
# --------- Example 1: p = 10 dimensional Gaussian ---------
# ----- 1) Specify Model -----
# a) General Graph Info
p <- 10 # number of variables
type = rep("g", p) # type of variables
level = rep(1, 10) # number of categories for each variable (1 = convention for continuous)
# b) Define interactions
factors <- list()
factors[[1]] <- matrix(c(1,2,
1,3,
4,5,
7,8), ncol=2, byrow = T) # 4 pairwise interactions
interactions <- list()
interactions[[1]] <- vector("list", length = 4)
# all pairwise interactions have value .5
for(i in 1:4) interactions[[1]][[i]] <- array(.5, dim=c(1, 1))
# c) Define Thresholds
thresholds <- vector("list", length = p)
thresholds <- lapply(thresholds, function(x) 0 ) # all means are zero
# d) Define Variances
sds <- rep(1, p) # All variances equal to 1
# ----- 2) Sample cases -----
data <- mgmsampler(factors = factors,
interactions = interactions,
thresholds = thresholds,
sds = sds,
type = type,
level = level,
N = 500,
nIter = 100,
pbar = TRUE)
# ----- 3) Recover model from sampled cases -----
set.seed(1)
mgm_obj <- mgm(data = data$data,
type = type,
level = level,
k = 2,
lambdaSel = "EBIC",
lambdaGam = 0.25)
mgm_obj$interactions$indicator # worked!
# --------- Example 2: p = 3 Binary model with one 3-way interaction ---------
# ----- 1) Specify Model -----
# a) General Graph Info
type = c("c", "c", "c")
level = c(2, 2, 2)
# b) Define Interaction
factors <- list()
factors[[1]] <- NULL # no pairwise interactions
factors[[2]] <- matrix(c(1,2,3), ncol=3, byrow = T) # one 3-way interaction
interactions <- list()
interactions[[1]] <- NULL
interactions[[2]] <- vector("list", length = 1)
# threeway interaction no1
interactions[[2]][[1]] <- array(0, dim = c(level[1], level[2], level[3]))
theta <- 2
interactions[[2]][[1]][1, 1, 1] <- theta # fill in nonzero entries
# thus: high probability for the case that x1 = x2 = x3 = 1
# c) Define Thresholds
thresholds <- list()
thresholds[[1]] <- rep(0, level[1])
thresholds[[2]] <- rep(0, level[2])
thresholds[[3]] <- rep(0, level[3])
# ----- 2) Sample cases -----
set.seed(1)
dlist <- mgmsampler(factors = factors,
interactions = interactions,
thresholds = thresholds,
type = type,
level = level,
N = 500,
nIter = 100,
pbar = TRUE)
# ----- 3) Check: Contingency Table -----
dat <- dlist$data
table(dat[,1], dat[,2], dat[,3]) # this is what we expected
# ----- 4) Recover model from sampled cases -----
mgm_obj <- mgm(data = dlist$data,
type = type,
level = level,
k = 3,
lambdaSel = "EBIC",
lambdaGam = 0.25,
overparameterize = TRUE)
mgm_obj$interactions$indicator # recovered, plus small spurious pairwise 1-2
# --------- Example 3: p = 5 Mixed Graphical Model with two 3-way interaction ---------
# ----- 1) Specify Model -----
# a) General Graph Info
type = c("g", "c", "c", "g")
level = c(1, 3, 5, 1)
# b) Define Interaction
factors <- list()
factors[[1]] <- NULL # no pairwise interactions
factors[[2]] <- matrix(c(1,2,3,
2,3,4), ncol=3, byrow = T) # no pairwise interactions
interactions <- list()
interactions[[1]] <- NULL
interactions[[2]] <- vector("list", length = 2)
# 3-way interaction no1
interactions[[2]][[1]] <- array(0, dim = c(level[1], level[2], level[3]))
interactions[[2]][[1]][,,1:3] <- rep(.8, 3) # fill in nonzero entries
# 3-way interaction no2
interactions[[2]][[2]] <- array(0, dim = c(level[2], level[3], level[4]))
interactions[[2]][[2]][1,1,] <- .3
interactions[[2]][[2]][2,2,] <- .3
interactions[[2]][[2]][3,3,] <- .3
# c) Define Thresholds
thresholds <- list()
thresholds[[1]] <- 0
thresholds[[2]] <- rep(0, level[2])
thresholds[[3]] <- rep(0, level[3])
thresholds[[4]] <- 0
# d) Define Variances
sds <- rep(.1, length(type))
# ----- 2) Sample cases -----
set.seed(1)
data <- mgmsampler(factors = factors,
interactions = interactions,
thresholds = thresholds,
sds = sds,
type = type,
level = level,
N = 500,
nIter = 100,
pbar = TRUE)
# ----- 3) Check: Conditional Means -----
# We condition on the categorical variables and check whether
# the conditional means match what we expect from the model:
dat <- data$data
# Check interaction 1
mean(dat[dat[,2] == 1 & dat[,3] == 1, 1]) # (compare with interactions[[2]][[1]])
mean(dat[dat[,2] == 1 & dat[,3] == 5, 1])
# first mean higher, ok!
# Check interaction 2
mean(dat[dat[,2] == 1 & dat[,3] == 1, 4]) # (compare with interactions[[2]][[2]])
mean(dat[dat[,2] == 1 & dat[,3] == 2, 4])
# first mean higher, ok!
}
Run the code above in your browser using DataLab