MI_fit_index_sim.R

library("parSim")
library("ggplot2")
library("tidyr")
library("dplyr")

res <- parSim(
  nIndicator = 5,
  nu_diff = c(0,0.5),
  sampleSize = c(50, 100,500,1000,5000),
  reps = 100,
  nCores = 16,
  expression = {
    # Load library:
    library("psychonetrics")
    library("mvtnorm")
    library("dplyr")
    
    # Simulate lambda and theta matrix:
    Lambda <- matrix(runif(nIndicator,0.5,1))
    Theta <- diag(runif(nIndicator,0.5,1))
    
    # Simulate alpha:
    alpha <- array(c(1,runif(1,0.5,2)),c(1,1,2))
    
    # Simulate Psi:
    psi <- array(c(1,runif(1,0.5,2)),c(1,1,2)) 
    
    # Simulate nu:
    nu <- array(0,c(nIndicator,1,2))
    nu[1,1,2] <- nu_diff
    
    # Variable names:
    vars <- paste0("v",seq_len(nIndicator))
    
    # Simulate data:
    data <- do.call(rbind,lapply(1:2,function(g){
      eta <- rnorm(sampleSize, mean = alpha[1,1,g], sd = sqrt(psi[1,1,g]))
      epsilon <- rmvnorm(sampleSize, sigma = Theta)
      Y <- as.data.frame(
        rep(1,sampleSize) %*% t(nu[,,g]) + eta %*% t(Lambda) + epsilon
      )
      names(Y) <- vars
      Y$group <- g
      return(Y)
    }))
    
    # Fit weak invariance model:
    weak <- lvm(data, vars=vars, groups = "group", lambda = matrix(1,nIndicator,1)) %>% groupequal("lambda") %>% runmodel
    fit_weak <- weak@fitmeasures
    
    # fit strong:
    strong <- weak %>% groupequal("nu") %>% runmodel
    fit_strong <- strong@fitmeasures
    
    # Reject strong?
    data.frame(
      chisq_0.05 = pchisq(fit_strong$chisq - fit_weak$chisq, fit_strong$df - fit_weak$df, lower.tail = FALSE) < 0.05,
      chisq_0.01 = pchisq(fit_strong$chisq - fit_weak$chisq, fit_strong$df - fit_weak$df, lower.tail = FALSE) < 0.01,
      AIC = fit_strong$aic.ll > fit_weak$aic.ll,
      BIC = fit_strong$bic > fit_weak$bic,
      RMSEA = fit_strong$rmsea > fit_weak$rmsea,
      CFI = fit_strong$cfi < fit_weak$cfi
    )
  }
)

# Plot results:
res_long_summary <- 
  res %>% 
  pivot_longer(chisq_0.05:CFI,names_to = "test") %>%
  group_by(nIndicator,nu_diff,sampleSize,test) %>%
   summarize(p = mean(value), n = n()) %>%
  mutate(se = sqrt(p * (1-p) / n)) %>% 
  mutate(lower = p - 1.96 * se, upper = p + 1.96 * se)

res_long_summary$nu_diff_factor <- paste0("Bias first item: ",res_long_summary$nu_diff)

ggplot(res_long_summary, aes(x=factor(sampleSize), y = p, min=lower, max=upper)) + 
  facet_grid(nu_diff_factor ~ test) + 
  scale_y_continuous(limits=c(0,1),breaks = seq(0,1,by=0.2), minor_breaks = seq(0,1,by=0.1)) +
  geom_hline(yintercept = 0.05) + 
  geom_hline(yintercept = 0.01, lty = 2) + 
  geom_errorbar() + 
  xlab("Sample Size") + 
  ylab("Probability of rejecting strong invariance") + 
  theme_bw() + 
  theme(axis.text.x = element_text(angle = 90))