if (FALSE) {
## We fit a pairwise and 3-order MGM to the mixed Autism dataset (?autism_data)
# 1) Fit Pairwise MGM
# Call mgm()
fit_k2 <- mgm(data = autism_data$data,
type = autism_data$type,
level = autism_data$lev,
k = 2) # ad most pairwise interacitons
# Weighted adjacency matrix
fit_k2$pairwise$wadj
# Visualize using qgraph()
library(qgraph)
qgraph(fit_k2$pairwise$wadj,
edge.color = fit_k2$pairwise$edgecolor,
layout = "spring",
labels = autism_data$colnames)
# 2) Fit MGM with pairwise & three-way interactions
fit_k3 <- mgm(data = autism_data$data,
type = autism_data$type,
level = autism_data$lev,
k = 3) # include all interactions up to including order 3
# List of estimated interactions
fit_k3$interactions$indicator
# 3) Plot Factor Graph
FactorGraph(object = fit_k3,
PairwiseAsEdge = FALSE,
labels = autism_data$colnames)
# 4) Predict values
pred_obj <- predict(fit_k3, autism_data$data)
head(pred_obj$predicted) # first six rows of predicted values
pred_obj$errors # Nodewise errors
## Here we illustrate why we need to overparameterize the design matrix to
## recover higher order interactions including categorical variables
# 1) Define Graph (one 3-way interaction between 3 binary variables)
# 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 <- .7
interactions[[2]][[1]][1, 1, 1] <- theta #weight theta for conf (1,1,1), weight 0 for all others
# c) Define Thresholds
thresholds <- list()
thresholds[[1]] <- c(0, 0)
thresholds[[2]] <- c(0, 0)
thresholds[[3]] <- c(0, 0)
# 2) Sample from Graph
iter <- 1
set.seed(iter)
N <- 2000
d_iter <- mgmsampler(factors = factors,
interactions = interactions,
thresholds = thresholds,
type = type,
level = level,
N = N,
nIter = 50,
pbar = TRUE)
# 3.1) Estimate order 3 MGM using standard parameterization
d_est_stand <- mgm(data = d_iter$data,
type = type,
level = level,
k = 3,
lambdaSel = "CV",
ruleReg = "AND",
pbar = TRUE,
overparameterize = FALSE,
signInfo = FALSE)
# We look at the nodewise regression for node 1 (same for all)
coefs_stand <- d_est_stand$nodemodels[[1]]$model
coefs_stand
# We see that nonzero-zero pattern of parameter vector does not allow us to infer whether
# interactions are present or not
# 3.2) Estimate order 3 MGM using overparameterization
d_est_over <- mgm(data = d_iter$data,
type = type,
level = level,
k = 3,
lambdaSel = "CV",
ruleReg = "AND",
pbar = TRUE,
overparameterize = TRUE,
signInfo = FALSE)
# We look at the nodewise regression for node 1 (same for all)
coefs_over <- d_est_over$nodemodels[[1]]$model
coefs_over # recovers exactly the 3-way interaction
# For more examples see https://github.com/jmbh/mgmDocumentation
}
Run the code above in your browser using DataLab