Gavagai.R

# Load libraries:
library("dplyr")
library("psychonetrics")

# Read table:
data <- read.table("fwbmlqsgavagai.dat")
names(data) <- c('fwb1', 'fwb2', 'fwb3', 'fwb4', 'fwb5', 
                 'mlqs1', 'mlqs2', 'mlqs3', 'mlqs4', 'mlqs5',
                 'gavagai')

# Factor structure:
lambda <- simplestructure(c(rep("fwb",5),rep("mlqs",5)))

# Variables:
vars <- c('fwb1', 'fwb2', 'fwb3', 'fwb4', 'fwb5', 
  'mlqs1', 'mlqs2', 'mlqs3', 'mlqs4', 'mlqs5')

# Latents:
lats <- c("fwb","mlqs")

# Configural model:
conf <- lvm(
  data,
  lambda = lambda,
  vars = vars,
  latents = lats,
  estimator = "FIML",
  group = "gavagai"
) %>% runmodel

# Inspect fit:
fit(conf) # OK

# Weak invariance:
weak <- conf %>% groupequal("lambda") %>% runmodel

# Compare:
compare(
  configural = conf, 
  weak = weak)

# Weak invariance should be rejected according to AIC and chi-square test

# We can try to go for partial week. First look at MIs:
weak %>% MIs(type = "free",matrices = "lambda",nonZero = TRUE)

# Free one factor loading:
partial_weak_1 <- weak %>% groupfree("lambda","mlqs4","mlqs") %>% runmodel

# Compare:
compare(
  configural = conf, 
  weak = weak,
  partial_weak = partial_weak_1)

# Still rejected according to chi-square and AIC..
partial_weak_1 %>% MIs(type = "free",matrices = "lambda",nonZero = TRUE)

# Free one factor more loading:
partial_weak_2 <- partial_weak_1 %>% groupfree("lambda","mlqs5","mlqs") %>% runmodel

# Compare:
compare(
  configural = conf, 
  weak = weak,
  partial_weak = partial_weak_2)

# still questionable... Let's try to use this model for strong invariance though:
strong <- partial_weak_2 %>% groupequal("nu") %>% runmodel

# Compare groups:
compare(
  configural = conf, 
  partial_weak = partial_weak_2,
  strong=strong)
# Also rejected