---
title: "mlVAR Permutation Test"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{mlVAR Permutation Test}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = ">"
)
```

The `mnet` package can be installed from Github using the `devtools` package:

```{r setup, message=FALSE, warning=FALSE}
# library(devtools)
# install_github("jmbh/mnet")
library(mnet)
```

To illustrate the permutation test for group differences in multilevel VAR models we use simulated data that is included in the `mnet` package. The data consists of two groups with 100 subjects each, with time series of three variables of length 100 for each subject. The data were generated from a mlVAR model in which one of the cross-lagged effects (2->1) is different across groups.

```{r}
data("ExampleData")

dim(ExampleData)
head(ExampleData)

```

The `mlVAR_CG()` function implements the permutation test and has the same input arguments as the `mlVAR()` function from the `mlVAR` package, except that it requires two data arguments for the two groups and requires specifying the number of permutations the test should be based on. The more permutations, the more precisely p-values can be computed.

We provide the two datasets to `mlVAR_CG()`, and as in `mlVAR()`, we specify the column names of the variables that should be modeled with the argument `vars` and indicate the column with the unique identifiers for each subject with the argument `idvar`. Finally, we indicate the column indicating the group membership with the argument `group`,

```{r}
set.seed(1) # reproducibility

timer <- proc.time()[3]

output <- mlVAR_GC(data = ExampleData, 
                   vars = c("V1", "V2", "V3"), 
                   idvar = "id",
                   groups = "group",
                   nCores = 1, # choose cores available on your machine
                   nP = 5, # Should be much more in practice, see preprint!
                   pbar = FALSE)

proc.time()[3] - timer
```

The perhaps most relevant output are the observed group differences and their associated p-values. In the current version, we consider group differences in between-person networks, in VAR networks (both fixed effects and random effects standard deviation) and the contemporaneous networks (both fixed effects and random effects standard deviation). The observed group differences in all these networks can be found in `out$EmpDiffs`. All group differences are Group1 - Group2.

```{r}
output$EmpDiffs
```

The associated p-values can be found in 

```{r}
output$Pval
```

For example, we see that the group difference for the cross lagged effect 2->1 is significant. Note that the p-value is here exactly zero, because we computed p-values based on empirical sampling distribution. If the test-statistic is far away from the mass of the sampling distribution under the null hypothesis (no group difference), then the density of the sampling distribution in this area is extremely low and therefore unlikely to be sampled. This leads to p-values that are exactly zero.

We can also inspect the exact sampling distributions. Let's again consider the cross lagged effect 2->1:

```{r, fig.align="center", fig.width=7, fig.height=5}
hist(output$SampDist$Lagged_fixed[1,2,], 
     xlim=c(-.5, .5), main="", xlab="Parameter Value")
abline(v=output$EmpDiffs$Lagged_fixed[1,2,], col="orange")
legend("topright", legend="Test-statistic", text.col="orange", bty="n")
```

We see that the test-statistic is not even overlapping with the support sampling distribution, which gives a p-value of zero as discussed above.