Nothing
# jonashalsbeck@gmail.com; August 2019
# Function to plot LHS legend
plotLegend <- function(labels,
TM,
margins,
ylim,
plot_y,
cex.label,
p = p) {
# Generate label vector
if(is.null(labels)) {
label_vec <- paste0(TM[, 1], " - ", TM[, 2])
} else {
tar_mat_label <- TM[ ,1:2]
tar_mat_label <- apply(tar_mat_label, 1:2, as.character)
for(i in 1:p) tar_mat_label[tar_mat_label == i] <- labels[i]
label_vec <- paste0(tar_mat_label[, 1], " - ", tar_mat_label[, 2])
}
# Plot
par(mar=margins)
plot.new()
plot.window(xlim = c(-1, 1), ylim = ylim)
text(0, plot_y, label_vec, cex = cex.label)
} # eoF
# Function to plot data
plotData <- function(TM,
axis.ticks,
ylim,
plot_y,
lwd.qtl,
cex.bg,
cex.mean,
margins,
bgcol = "white") {
# Some settings
xlim <- range(axis.ticks)
# Setup canvas
par(mar=margins)
plot.new()
plot.window(xlim = xlim, ylim = ylim)
# background color
rect(par("usr")[1],par("usr")[3],par("usr")[2],par("usr")[4],col = bgcol, border = FALSE)
# rect(xleft = xlim[1], ybottom = ylim[1], xright = xlim[2], ytop = ylim[2], col = bgcol, border = FALSE)
if(is.null(axis.ticks)) axis.ticks <- round(seq(xlim[1], xlim[2], length = 5), 2)
axis(3, axis.ticks, lwd=0)
abline(h = plot_y, col = "grey")
if(0 %in% axis.ticks) abline(v = 0, lty=2, col = "black") # zero line
# Plot quantiles
segments(x0 = TM[, 4],
y0 = plot_y,
x1 = TM[, 5],
y1 = plot_y, lwd = lwd.qtl)
# Plot prop>0
points(TM[, 3], plot_y, pch=19, col=bgcol, cex = cex.bg)
points(TM[, 3], plot_y, pch=21, col="black", cex = cex.bg)
text(TM[, 3], plot_y, TM[, 6], cex = cex.mean)
} # eoF
plotRes <- function(object,
quantiles = c(.05, .95),
labels = NULL,
decreasing = TRUE,
cut = NULL,
cex.label = .75,
lwd.qtl = 2,
cex.mean = .55,
cex.bg = 2.7,
axis.ticks = c(-.5, -.25, 0, .25, .5, .75, 1),
axis.ticks.mod = NULL,
layout.width.labels = .2,
layout.gap.pw.mod = .15,
table = FALSE)
{
# ----------------------------------------------------------------------
# ---------- Aux variables ---------------------------------------------
# ----------------------------------------------------------------------
# Fill in defaults
# Get vars
mod <- object$call$object$call$moderators
type <- object$call$object$call$type
p <- length(type)
nB <- object$call$nB
binarySign <- object$call$object$call$binarySign
ind_cat <- object$call$object$call$ind_cat
ind_binary <- object$call$object$call$ind_binary
# ----------------------------------------------------------------------
# ---------- Input Checks ----------------------------------------------
# ----------------------------------------------------------------------
# Input checks: basic
if(!("core" %in% class(object))) stop("plotRes() currently only supports resampled mgm() objects.")
if(!("resample" %in% class(object))) stop("PplotRes() only takes resample objects as input (see ?resample).")
# if(missing(quantiles)) stop("No quantiles specified.")
if(object$call$object$call$k > 3) stop("Currently only implemented with mgms with k <= 3.")
# Check whether there is only a single moderator
if(class(mod)[1] == "numeric") if(length(mod) > 1) stop("plotRes() is currently only implemented for MNNs with a single moderator.")
if(class(mod)[1] == "matrix") {
# find whether one variable is in each 3-way interaction; if yes, that's the one we treat as the moderator
n_mod <- nrow(mod)
m_check <- matrix(NA, n_mod, p)
for(i in 1:p) for(m in 1:n_mod) m_check[m, i] <- i %in% mod[m, ]
v_check <- apply(m_check, 2, function(x) all(x == TRUE))
mod_fix <- which(v_check)
if(is.null(mod_fix)) stop("plotRes() is currently only implemented for MNNs with a single moderator.")
# pretend that we have all possible moderation effect of mod; this will result in unmodeled moderation effects to be set to zero (as they should be)
mod <- mod_fix
}
# Define set of continuous and binary variables: for interactions between these we can assign a sign
# Depends on binarySign
if(binarySign) {
set_signdefined <- c(which(type == 'p'), which(type == 'g'), ind_cat[ind_binary])
} else {
set_signdefined <- c(which(type == 'p'), which(type == 'g'))
}
# ----------------------------------------------------------------------
# ---------- A.1) mgm objects (unmoderated) ----------------------------
# ----------------------------------------------------------------------
if(is.null(mod)) {
# ---------- Preprocessing ----------
# Get basic info
dims <- dim(object$bootParameters)
p <- dims[1]
nB <- dims[3]
n_pars <- p*(p-1) / 2
# Collapse into edge x property matrix
tar_mat <- matrix(NA, nrow=n_pars, ncol = 6)
colnames(tar_mat) <- c("Variable A", "Variable B", "Mean", "qtl_low", "qtl_high", "propLtZ")
counter <- 1
for(row in 1:p) {
for(col in row:p) {
if(row!=col){
# Variable ids
tar_mat[counter, 1] <- row
tar_mat[counter, 2] <- col
# Quantiles
qtls <- quantile(object$bootParameters[row, col, ], probs = quantiles)
tar_mat[counter, 3] <- mean(object$bootParameters[row, col, ])
tar_mat[counter, 4] <- qtls[1]
tar_mat[counter, 5] <- qtls[2]
tar_mat[counter, 6] <- mean(abs(object$bootParameters[row, col, ]) > 0) # proportion estimates > 0
# update counter
counter <- counter + 1
}
}
}
tar_mat <- round(tar_mat, 2)
# Order
tar_mat <- tar_mat[order(tar_mat[,3], decreasing = decreasing), ]
# Subset (cut)
if(is.null(cut)) {
TM <- tar_mat
} else {
TM <- tar_mat[cut, ]
}
if(!table) {
# ---------- Plotting ----------
# Compute aux variables for plotting
n_rows <- nrow(TM)
ylim <- c(0, 1)
plot_y <- seq(ylim[2], ylim[1], length = n_rows)
margins <- c(0, .5, 3, .5)
# ----- Setup layout ----
lmat <- matrix(1:2, nrow=1)
lo <- layout(lmat, widths = c(layout.width.labels, 1))
# ----- Part A: Legend ----
plotLegend(labels = labels,
TM = TM,
margins = margins,
ylim = ylim,
plot_y = plot_y,
cex.label = cex.label, p=p)
# ----- Part B: Data ----
plotData(TM = TM,
axis.ticks = axis.ticks,
ylim = ylim,
plot_y = plot_y,
lwd.qtl = lwd.qtl,
cex.bg = cex.bg,
cex.mean = cex.mean,
margins=margins)
# Return table insteaed
} else {
return(TM)
}
} # end if: moderation?
# ----------------------------------------------------------------------
# ---------- A.1) mgm objects (moderated; 1 moderator) -----------------
# ----------------------------------------------------------------------
if(!is.null(mod)) {
# ---------- Preprocessing ----------
# ----- Get estimates out of model object -----
m_ind_allpw <- t(combn(1:p, 2)) # list all possible 2-way interactions
n_pw <- nrow(m_ind_allpw) # how many?
m_pw <- m_mod <- matrix(0, nrow=n_pw, ncol=nB) # Create storage for nB pairwise and moderation effects
# loop through bootstrapped objects
for(b in 1:nB) {
mod_b <- object$models[[b]]
mod_b_pw <- matrix(mod_b$interactions$indicator[[1]], ncol=2)
mod_b_m <- matrix(mod_b$interactions$indicator[[2]], ncol=3)
# Loop pairwise & get pairwise
if(nrow(mod_b_pw) > 0) for(i in 1:nrow(mod_b_pw)) {
ind_pw <- which(apply(m_ind_allpw, 1, function(x) sum(x %in% mod_b_pw[i, ])) == 2)
m_pw[ind_pw, b] <- mod_b$interactions$weightsAgg[[1]][[i]][1]
# Add sign if applicable
if(mod_b$interactions$signs[[1]][[i]] == -1) m_pw[ind_pw, b] <- m_pw[ind_pw, b] * -1
# # add sign if applicable
# if(all(type[mod_b_pw] == "g")) {
# v_est <- unlist(mod_b$interactions$weights[[1]][[i]])
# v_sign <- sign(v_est)
# sign_sel <- as.numeric(names(which.max(table(v_sign)))) # ugly!
# m_pw[ind_pw, b] <- m_pw[ind_pw, b] * sign_sel
# }
} # end: pw int
# Loop moderation & get moderation effects
if(nrow(mod_b_m) > 0) for(i in 1:nrow(mod_b_m)) {
ind_mod <- which(apply(m_ind_allpw, 1, function(x) sum(x %in% mod_b_m[i, ][-which(mod_b_m[i, ]==mod)])) == 2)
m_mod[ind_mod, b] <- mod_b$interactions$weightsAgg[[2]][[i]][1] #* mod_b$interactions$signs[[2]][[i]] # add sign
# Add sign if applicable !!!! THIS DOES NOT WORK !!! because no signs are defined in Reg2Graph() for 3-way interactions; there was a reason I did the complicated thing below
if(mod_b$interactions$signs[[2]][[i]] == -1) m_mod[ind_mod, b] <- m_mod[ind_mod, b] * -1
# add sign if applicable
if(all(mod_b_m[i, ] %in% set_signdefined)) {
v_est <- rep(NA, 3)
for(k in 1:3) {
if(type[mod_b_m[i, k]] == "c") {
v_est[k] <- mod_b$interactions$weights[[2]][[i]][[k]][[2]]
} else {
v_est[k] <- mod_b$interactions$weights[[2]][[i]][[k]]
}
}
v_sign <- sign(v_est)
sign_sel <- as.numeric(names(which.max(table(v_sign)))) # ugly!
m_mod[ind_mod, b] <- m_mod[ind_mod, b] * sign_sel
}
} # end: mod
} # end: bootstrap it
# ----- Compute quantile / median / prop>0 -----
tar_mat_pw <- tar_mat_mod <- matrix(NA, nrow=n_pw, ncol = 6)
colnames(tar_mat_pw) <- colnames(tar_mat_pw) <- c("Variable A", "Variable B", "Mean", "qtl_low", "qtl_high", "propLtZ")
# Pairwise
tar_mat_pw[, 1:2] <- m_ind_allpw
tar_mat_pw[, 3] <- rowMeans(m_pw)
tar_mat_pw[, 4:5] <- t(apply(m_pw, 1, function(x) quantile(x, probs = quantiles)))
tar_mat_pw[, 6] <- apply(m_pw, 1, function(x) mean(x!=0) )
tar_mat_pw <- round(tar_mat_pw, 2)
# Moderation
tar_mat_mod[, 1:2] <- m_ind_allpw
tar_mat_mod[, 3] <- rowMeans(m_mod)
tar_mat_mod[, 4:5] <- t(apply(m_mod, 1, function(x) quantile(x, probs = quantiles)))
tar_mat_mod[, 6] <- apply(m_mod, 1, function(x) mean(x!=0) )
tar_mat_mod <- round(tar_mat_mod, 2)
# Order
ord <- order(tar_mat_pw[,3], decreasing = decreasing)
tar_mat_pw <- tar_mat_pw[ord, ]
tar_mat_mod <- tar_mat_mod[ord, ]
# Subset (cut argument)
if(is.null(cut)) {
TM_pw <- tar_mat_pw
TM_mod <- tar_mat_mod
} else {
TM_pw <- tar_mat_pw[cut, ]
TM_mod <- tar_mat_mod[cut, ]
}
if(!table) {
# ---------- Plotting ----------
n_rows <- nrow(TM_pw)
ylim <- c(0, 1)
plot_y <- seq(ylim[2], ylim[1], length = n_rows)
margins <- c(0, .5, 3, .5)
# ----- Setup layout -----
lmat <- rbind(c(1,2,7,3),
c(4,5,8,6))
lo <- layout(lmat,
widths = c(layout.width.labels, 1, layout.gap.pw.mod, 1),
heights = c(.07, .9))
# ----- 0) Plot Top Legend -----
par(mar=rep(0, 4))
plot.new() # fill 1
plot.new()
plot.window(xlim=c(0, 1), ylim=c(0,1))
text(0.5, 0.5, "Pairwise effects", cex=1.5, adj = 0.5)
plot.new()
plot.window(xlim=c(0, 1), ylim=c(0,1))
text(0.5, 0.5, "Moderation effects", cex=1.5, adj = 0.5)
# ----- A) Plot LHS Legend -----
plotLegend(labels = labels,
TM = TM_pw,
margins = margins,
ylim = ylim,
plot_y = plot_y,
cex.label = cex.label, p=p)
# ----- B) Plot Pairwise effects -----
plotData(TM = TM_pw,
axis.ticks = axis.ticks,
ylim = ylim,
plot_y = plot_y,
lwd.qtl = lwd.qtl,
cex.bg = cex.bg,
cex.mean = cex.mean,
margins = margins)
# ----- C) Plot Moderation effects -----
if(is.null(axis.ticks.mod)) axis.ticks.mod <- axis.ticks
plotData(TM = TM_mod,
axis.ticks = axis.ticks.mod,
ylim = ylim,
plot_y = plot_y,
lwd.qtl = lwd.qtl,
cex.bg = cex.bg,
cex.mean = cex.mean,
bgcol = "white",
margins = margins)
# Return table instead
} else {
# make table
colnames(TM_mod)[3:6] <- c("Mod_Mean", "Mod_qtl_low", "Mod_qtl_high", "Mod_propLtZ")
out_table <- cbind(TM_pw, TM_mod)
return(out_table)
}
} # end if: moderation?
# ----------------------------------------------------------------------
# ---------- B) mvar objects -------------------------------------------
# ----------------------------------------------------------------------
# To do ...
} # eoF
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.