Psychonetrics RNM example

# Install psychonetrics:
# devtools::install_github("sachaepskamp/psychonetrics")
library("psychonetrics")

# Set the seed:
set.seed(1)

# Let's simulate some data, 2-factor model with 10 indicators each and residual chain graph.
# Generate factor loadings:
lambda <- matrix(0,20,2)
lambda[1:10,1] <- lambda[11:20,2] <- 1
lambda <- lambda * runif(40)

# Generate latent varcov:
psi <- matrix(0.5,2,2)
diag(psi) <- 1

# Generate residual network:
library("bootnet")
ResidNet <- genGGM(20)
Theta <- cov2cor(solve(diag(20) - ResidNet))

library("mvtnorm")
Data <- rmvnorm(1000, sigma = lambda %*% psi %*% t(lambda) + Theta)

# Let's put it in psychonetrics. We need only a lambda matrix with 1s indicating which 
# parameters are free and zeroes indicating which are fixed (higher integers can imply 
# equality constraints). 
modelLambda <- 1*(lambda!=0)

# We can use the raw data, but let's use the cov matrix instead.
# Note though, cov matrix has to be max likelihood estimated. Divide 
# by n not n-1, this can be corrected using n/(n-1) as multiplyer to 
# R's default function:
covMat <- (1000/999) * cov(Data)

# Form the model:
model <- rnm(lambda = modelLambda, covs = covMat, nobs = 1000)

# Now we can run the model:
model <- model %>% runmodel

# And inspect fit:
print(model)
model %>% fit

# pretty lousy. We can look at MIs:
model %>% MIs # Indicates quite some edges with V17 -- V16 the most important. Let's free it:
model2 <- model %>% freepar("omega_epsilon","V17","V16")

# And run again:
model2 <- model2 %>% runmodel

# Compare to previous model:
compare(model,model2)

# Much better fit. This can be tedious manually, so let's use automated search by freeing the top MI (network only).
# By default, the search stops when no MI is significant at alpha = 0.01 or if EBIC (gamma = 0.25) is not improved
model3 <- model %>% stepup(criterion = "ebic") # Default criterion is bic, still have to test which is the best.

# Let's plot the result (this is ugly still, no plot method yet):
library("qgraph")
layout(t(1:2))
qgraph(ResidNet, cut = 0, fade = FALSE,
       title = "True residual network")
qgraph(model3@modelmatrices$`1`$omega_epsilon, cut = 0, fade = FALSE,
       title = "Estimated residual network")

# Looks great! Only missing one edge (BIC selection adds one false positive, also good). the model fits very well:
print(model3)
model3 %>% fit
compare(original = model, searched = model3)

# We can also look at parameter estimates:
model3 %>% parameters