I am attempting to simulate data using the simstudy package and then calculate coefficient omega with the psych package. I am using for-loops to run through different values for within cluster variation, inter-item correlation, and a number of iterations.
I began testing this with 3-10 iterations and kept getting an error from lavaan saying the model did not converge. When I opened up the output, I could see the for loop was actually reverting back to previous values then stopping. For example, with 10 iterations per condition here is what happened:
var = 0 + all values for rho - ran successfully
var = 0.5 + all values for rho - ran successfully
var = 1.0 + all values for rho - ran successfully
var = 1.5 + rho = 0.1 - 7 successful iterations, then reverted back to var = 1.0 + rho = 0.6 for the remaining 3 iterations in that loop then entire loop stopped prematurely.
I tried updating R, R studio, and R packages. I ran it with and without the tryCatch for errors with Omega. I also tried adjusting the number of iterations and found with more iterations the issue occurs later in the simulations, however, when I put iterations to 1000 (my end goal) and all the output became NAs.
library(simstudy)
library(dirmult)
library(tidyverse)
library(dplyr)
library(readr)
library(psych)
library(GPArotation)
library(multilevel)
library(lavaan)
## creating an empty data frame for results ##
results_df <- data.frame(iteration=numeric(),
within_cluster_variance=numeric(),
inter_item_rho=numeric(),
SL_Omega_raw=numeric())
set.seed(123)
var_value <- c(0.0, 0.5, 1.0, 1.5, 2.0, 2.5)
rho_value <- c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6)
for(v in var_value){
for(r in rho_value){
for(i in 1:5){
# LEVEL 2 #
class.level <- defData(varname = "class_zscore", dist = "normal",
formula = 0, variance = 1, id = "Class_ID")
class.level <- defData(class.level, varname = "Student_Count", dist = "noZeroPoisson", formula = 20)
class.data <- genData(100, class.level)
# LEVEL 1 #
gen.student <- defDataAdd(varname = "student_EXTERNALIZE_score", dist = "normal",
formula = "class_zscore", variance = v)
dtClass <- genCluster(class.data, "Class_ID", numIndsVar = "Student_Count", level1ID = "Student_ID")
dtClass <- addColumns(gen.student, dtClass)
baseprobs <- matrix(c(
0.973, 0.018, .006, 0.003,
0.829, 0.095, .050, 0.026,
0.765, 0.115, .069, 0.051,
0.882, 0.068, .032, 0.018,
0.717, 0.106, .081, 0.096,
0.880, 0.062, .038, 0.020,
0.905, 0.045, .034, 0.016),
nrow = 7, byrow = TRUE)
student.items <- genOrdCat(dtClass , adjVar = "student_EXTERNALIZE_score",
baseprobs, prefix = "Item",
asFactor=FALSE, idname = "Student_ID",
corstr = "cs", rho = r)
student.items->SIM_DATA
SIM_DATA[,6:12]->ITEMS
# Single-level Continuous Omega #
# Initialize SL_Omega_values with "Convergence failed" #
SL_Omega_values <- data.frame(omega.tot = "Convergence failed")
# Attempt to calculate omega total reliability coefficients #
tryCatch({
psych::omega(ITEMS, nfactors = 1, poly = FALSE, plot = FALSE, lavaan = TRUE) -> SL_omega
SL_Omega_values <- as.data.frame(SL_omega$omega.tot)
}, error = function(e) {
# If convergence fails, assign "Convergence failed" to SL_Omega_values[1,1] #
SL_Omega_values[1,1] <- "Convergence failed"
})
results_df[i,] <- (c(i, v, r, SL_Omega_values[1,1]))
write.csv(results_df, paste0("results_df","_",v,"_",r,"_",".csv"))
results_df}}}```
I think the problem is that your tryCatch() isn't saving the values you want in the right way. You need SL_Omega_values to be the output from tryCatch() rather than trying to set that value within the function. Try this:
library(simstudy)
library(dirmult)
library(tidyverse)
library(dplyr)
library(readr)
library(psych)
library(GPArotation)
library(multilevel)
library(lavaan)
## creating an empty data frame for results ##
results_df <- data.frame(iteration=numeric(),
within_cluster_variance=numeric(),
inter_item_rho=numeric(),
SL_Omega_raw=numeric())
set.seed(123)
var_value <- c(0.0, 0.5, 1.0, 1.5, 2.0, 2.5)
rho_value <- c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6)
for(v in var_value){
for(r in rho_value){
for(i in 1:5){
# LEVEL 2 #
class.level <- defData(varname = "class_zscore", dist = "normal",
formula = 0, variance = 1, id = "Class_ID")
class.level <- defData(class.level, varname = "Student_Count", dist = "noZeroPoisson", formula = 20)
class.data <- genData(100, class.level)
# LEVEL 1 #
gen.student <- defDataAdd(varname = "student_EXTERNALIZE_score", dist = "normal",
formula = "class_zscore", variance = v)
dtClass <- genCluster(class.data, "Class_ID", numIndsVar = "Student_Count", level1ID = "Student_ID")
dtClass <- addColumns(gen.student, dtClass)
baseprobs <- matrix(c(
0.973, 0.018, .006, 0.003,
0.829, 0.095, .050, 0.026,
0.765, 0.115, .069, 0.051,
0.882, 0.068, .032, 0.018,
0.717, 0.106, .081, 0.096,
0.880, 0.062, .038, 0.020,
0.905, 0.045, .034, 0.016),
nrow = 7, byrow = TRUE)
student.items <- genOrdCat(dtClass , adjVar = "student_EXTERNALIZE_score",
baseprobs, prefix = "Item",
asFactor=FALSE, idname = "Student_ID",
corstr = "cs", rho = r)
student.items->SIM_DATA
SIM_DATA[,6:12]->ITEMS
# Single-level Continuous Omega #
# Initialize SL_Omega_values with "Convergence failed" #
SL_Omega_values <- data.frame(omega.tot = "Convergence failed")
# Attempt to calculate omega total reliability coefficients #
SL_Omega_values <- tryCatch({
psych::omega(ITEMS, nfactors = 1, poly = FALSE, plot = FALSE, lavaan = TRUE) -> SL_omega
as.data.frame(SL_omega$omega.tot)
}, error = function(e) {
# If convergence fails, assign "Convergence failed" to SL_Omega_values[1,1] #
"Convergence failed"
})
results_df[i,] <- (c(i, v, r, SL_Omega_values[1,1]))
write.csv(results_df, paste0("results_df","_",v,"_",r,"_",".csv"))
results_df}}}