---
title: "Network-Based R-statistics for mixed-effects models"
#author: "Zeus Gracia-Tabuenca"
#date: "July 2020"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{NBR-LME}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
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)`.

```{r setup}
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.

```{r upper triangle}
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`.

```{r input network, fig.align = "center"}
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.


```{r NBR-LME, eval = FALSE}
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.


```{r multicore NBR-LME, eval = FALSE}
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.

```{r NBR-LME results}
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.

```{r component display, fig.align = "center"}
# 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`).

```{r component cum-pval, fig.height = 3, fig.width = 5, fig.align = "center"}
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")
```