knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
This vignette documents the implementation of NBR 0.1.3 for linear mixed effect (LME) models.
We will analyze the voles
dataset, which contains a matrix of 96 rows (sessions) and 123 columns (variables). The first three variables include phenotypic information of the subjects/sessions (1: subject ID; 2: Sex; 3: Session 1-3), the remaining 120 variables include the upper triangle edges of a network of 16 brain regions (fMRI functional connectivity).
NOTE: for more detail of the dataset execute help(voles)
.
library(NBR) data("voles") brain_labs <- NBR:::voles_roi dim(voles) head(voles)[1:8]
Here we can obtain the corresponding pairwise interaction of the brain network for each edge.
nnodes <- length(brain_labs) tri_pos <- which(upper.tri(matrix(nrow = nnodes, ncol = nnodes)), arr.ind = T) head(tri_pos)
IT'S VERY IMPORTANT that the order of the columns containing the network data matches with the order of the upper triangle of the network matrix.
Let's plot the average network with lattice::levelplot
.
library(lattice) avg_mx <- matrix(0, nrow = nnodes, ncol = nnodes) avg_mx[upper.tri(avg_mx)] <- apply(voles[-(1:3)], 2, function(x) mean(x, na.rm=TRUE)) avg_mx <- avg_mx + t(avg_mx) # Set max-absolute value in order to set a color range centered in zero. flim <- max(abs(avg_mx)) levelplot(avg_mx, main = "Average", ylab = "ROI", xlab = "ROI", at = seq(-flim, flim, length.out = 100))
The next step is to check the dataset to be tested edgewise. In this case we are going to test if the variables Sex
, Session
, and their interaction (Sex:Session
) have any effect related to the brain networks. Since every subject was assessed in three different sessions, we should add the intercept and the Session
term as random effects adding the random formula ~ 1+Session|id
, where id
accounts for the subject label.
set.seed(18900217) before <- Sys.time() library(nlme) nbr_result <- nbr_lme_aov(net = voles[,-(1:3)], nnodes = 16, idata = voles[,1:3], nperm = 5, mod = "~ Session*Sex", rdm = "~ 1+Session|id", na.action = na.exclude) after <- Sys.time() show(after-before)
Although five permutations is quite low to obtain a proper null distribution, we can see that they take several seconds to be performed. So we suggest paralleling to multiple CPU cores with the cores
argument.
set.seed(18900217) before <- Sys.time() library(nlme) library(parallel) nbr_result <- nbr_lme_aov( net = voles[,-(1:3)], nnodes = 16, idata = voles[,1:3], nperm = 1000, nudist = T, mod = "~ Session*Sex", rdm = "~ 1+Session|id", cores = detectCores(), na.action = na.exclude ) after <- Sys.time() show(after-before)
This may elapse approximately 15 minutes in an Intel(R) Core(TM) i7-8700 CPU @ 3.20GHz with 12 cores. But we can load those results instead of running them again.
nbr_result <- NBR:::voles_nbr show(nbr_result$fwe)
If we observed the Family-Wise Error (FWE) probabilities of the observed components, only the component 1 in the Session
term is lower than the nominal alpha of p < 0.05. The table shows the probabilities associated with: 1) the number of connected edges, and 2) the sum of the strength of the edges. In this case, we will use the sum of strengths, but you can choose depending on your research question.
Let's display the FWE-corrected component.
# Plot significant edges edge_mat <- array(0, dim(avg_mx)) edge_mat[nbr_result$components$Session[,2:3]] <- 1 levelplot(edge_mat, col.regions = rev(heat.colors(100)), main = "Component", ylab = "ROI", xlab = "ROI")
Lastly, if we are not sure if 1000 permutations are enough we can plot the cumulative p-value (black line) with its corresponding binomial marginal error (green lines). To do so, you just need to set TRUE for the return null distribution argument (nudist
).
null_ses_str <- nbr_result$nudist[,2] # Null distribution for Session strength obs_ses_str <- nbr_result$fwe$Session[,4] # Observed Session strength nperm <- length(null_ses_str) cumpval <- cumsum(null_ses_str >= obs_ses_str)/(1:nperm) # Plot p-value stability plot(cumpval, type="l", ylim = c(0,0.06), las = 1, xlab = "Permutation index", ylab = "p-value", main = "Cumulative p-value for Session strength") abline(h=0.05, col="red", lty=2) # Add binomial marginal error mepval <- 2*sqrt(cumpval*(1-cumpval)/1:nperm) lines(cumpval+mepval, col = "chartreuse4") lines(cumpval-mepval, col = "chartreuse4")
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.