psychonetrics_multilevel_test.R

# FIML multilevel
# Model described in https://psyarxiv.com/8ha93/
# See also http://statmodel.com/bmuthen/articles/Article_055.pdf

# Load packages:
library("lavaan")
library("psychonetrics")
library("bootnet")
library("mvtnorm")
library("qgraph")

### Input ###

# Number of clusters:
nCluster <- 1000

# Number of people in each cluster (2 = dyad)
nInClusterMin <- 2
nInClusterMax <- 2

# Number of nodes:
nNode <- 8

# (number of people in each cluster) x (number of nodes) should be low! (less than 100 max)

### Script ###

# Simulate sample sizes:
nInCluster <- sample(nInClusterMin:nInClusterMax, nCluster, replace = TRUE)

# Within network:
Wnet <- genGGM(nNode)

# Between network:
Bnet <- genGGM(nNode, p = 1)

# Simulate datasets:
gen <- ggmGenerator()
Means <- gen(nCluster, Bnet)
Datas <- list()
for (i in 1:nCluster){
  Datas[[i]] <- t(t(gen(nInCluster[i], Wnet)) + Means[i,])
  Datas[[i]] <- as.data.frame(Datas[[i]])
  names(Datas[[i]]) <- paste0("V",1:nNode)
  Datas[[i]]$group <- i
  Datas[[i]]$subjectInGroup <- 1:nInCluster[i]
}

# Combine datasets:
Data <- do.call(rbind, Datas)

# To wide format for panel data model:
library("tidyr")
gathered <- Data %>% gather(variable, value, paste0("V",1:nNode))
wideData <- gathered %>% pivot_wider(id_cols = group, names_from = c(variable,subjectInGroup), values_from=value)


# Now make the design matrix:
vars <- names(wideData)[-1]
varsMat <- do.call(rbind,strsplit(vars, split = "_"))
rowVars <- unique(varsMat[,1])
colVars <- unique(varsMat[,2])
design <- matrix(NA, length(rowVars), length(colVars))
for (i in seq_along(rowVars)){
  for (j in seq_along(colVars)){
    varName <- paste0("^",rowVars[i],"_",colVars[j],"$")
    whichVar <- which(grepl(varName, names(wideData)))
    if (length(whichVar) == 1){
      design[i,j] <- paste0(rowVars[i],"_",colVars[j])
    }
  }
}

# The matrix design could also have been formed by hand...

# Form model (make sure not to estimate beta, because this is not longitudinal data):
mod <- panelgvar(wideData, vars = design, estimator = "FIML", beta = "empty")

# Run model:
mod <- mod %>% runmodel

# Prune model and search:
mod <- mod %>% prune(alpha = 0.05, recursive = FALSE) %>% modelsearch

### Output ###

# Fit (DF is not really correct, so interpret with care)
mod %>% fit

# Obtain results:
West <- getmatrix(mod, "omega_zeta_within")
Best <- getmatrix(mod, "omega_zeta_between")

# Plot:
layout(matrix(1:4,2,2,byrow=TRUE))
qgraph(Wnet, title = "True within cluster network", theme = "colorblind")
qgraph(West, title = "Estimated within cluster network", theme = "colorblind")

qgraph(Bnet, title = "True between cluster network", theme = "colorblind")
qgraph(Best, title = "Estimated between cluster network", theme = "colorblind")