# NOT RUN {
# }
# NOT RUN {
## We fit a pairwise and 3-order MGM to the mixed Autism dataset (?autism_data)
# 1) Fit Pairwise MGM
# Call mgm()
fit_d2 <- mgm(data = autism_data$data,
type = autism_data$type,
level = autism_data$lev,
k = 2) # ad most pairwise interacitons
# Weighted adjacency matrix
fit_d2$pairwise$wadj
# Visualize using qgraph()
library(qgraph)
qgraph(fit_d2$pairwise$wadj,
edge.color = fit_d2$pairwise$edgecolor,
layout = 'spring',
labels = autism_data$colnames)
# 2) Fit MGM with pairwise & threeway interactions
fit_d3 <- 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_d3$rawfactor$indicator
# Compute Label vector
labels <- c(autism_data$colnames, rep('', sum(fit_d3$factorgraph$nodetype)))
# Visualize factor graph using qgraph()
library(RColorBrewer)
nodeColors <- brewer.pal(3, 'Set1')
qgraph(fit_d3$factorgraph$graph,
color = nodeColors[fit_d3$factorgraph$order+1],
edge.color = fit_d3$factorgraph$edgecolor,
layout = 'spring',
labels = labels,
shape = c('circle', 'square')[fit_d3$factorgraph$nodetype+1],
vsize = c(8, 4)[fit_d3$factorgraph$nodetype+1])
# 3) Predict values
pred_obj <- predict(fit_d3, autism_data$data)
head(pred_obj$predicted) # 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 usind 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 usind 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
# }
# NOT RUN {
# }
Run the code above in your browser using DataCamp Workspace