Simulation from:
User-defined helper function for the generate-analyse steps to generate appropriate syntax for lavaan (could be placed in an external .R file and sourced in prior to running runSimulation() instead).
#' @param J number of variables
#' @param analyse logical; is syntax being used to data generation or analysis?
#' @examples
#' cat(genLavaanSyntax(factors=1, indicators=10))
#' cat(genLavaanSyntax(factors=2, indicators=10))
#' cat(genLavaanSyntax(factors=1, indicators=10, analyse = TRUE))
#' cat(genLavaanSyntax(factors=2, indicators=10, analyse = TRUE))
#'
genLavaanSyntax <- function(factors, indicators, analyse = FALSE){
ret <- if(factors == 1){
if(analyse){
paste0(paste0('f1 =~ NA*x1 + ', paste0(paste0('x', 2:indicators),
collapse=' + ')),
'\nf1 ~~ 1*f1')
} else {
paste0(paste0('f1 =~ ', paste0(rep(.7, indicators),
paste0('*x', 1:indicators),
collapse=' + '), ' \n'),
paste0(sprintf('x%s ~~ 0.51*x%s', 1:indicators, 1:indicators),
collapse=' \n'))
}
} else if(factors == 2){
if(analyse){
paste0(paste0('f1 =~ NA*x1 + ',
paste0(paste0('x', 2:indicators), collapse=' + ')),
paste0(sprintf('\nf2 =~ NA*x%s + ', indicators+1),
paste0(paste0('x', 2:indicators + indicators), collapse=' + ')),
'\nf1 ~~ 1*f1 \nf2 ~~ 1*f2 \nf1 ~~ f2')
} else {
paste0(paste0('f1 =~ ', paste0(rep(.7, indicators),
paste0('*x', 1:indicators),
collapse=' + '), ' \n'),
paste0('f2 =~ ', paste0(rep(.7, indicators),
paste0('*x', 1:indicators + indicators),
collapse=' + '), ' \n'),
'f1 ~~ .3*f2 \n',
paste0(sprintf('x%s ~~ 0.51*x%s', 1:indicators, 1:indicators),
collapse=' \n'))
}
} else stop('factors input is incorrect')
ret
}library(SimDesign)
Design <- createDesign(N = c(100, 200, 500, 1000),
categories = c(2, 5),
skewness_kurtosis = list(c(0, 0), c(.75, 1.75), c(.75, 3.75),
c(1.25, 1.75), c(1.25, 3.75)),
factors = c(1, 2),
indicators = c(5, 10),
estimator = c('WLSMV', 'WLS'),
# remove poorly converging combinations
subset = !(estimator == 'WLS' & N %in% c(100, 200) &
factors == 2 & indicators == 10))
######################################################
Generate <- function(condition, fixed_objects = NULL) {
Attach(condition)
syntax <- genLavaanSyntax(factors=factors, indicators=indicators)
cdat <- simulateData(syntax, model.type='cfa', sample.nobs=N,
skewness=skewness_kurtosis[1L],
kurtosis=skewness_kurtosis[2L])
tau <- if(categories == 5)
c(-1.645, -0.643, 0.643, 1.645) else 0
# data generation fix described in Flora's (2002) unpublished dissertation
if(categories == 5 && all(skewness_kurtosis == c(1.25, 1.75)))
tau[1] <- -1.125
dat <- apply(cdat, 2, function(x, tau){
dat <- numeric(length(x))
for(i in 1:length(tau))
dat[x > tau[i]] <- i
dat
}, tau=tau)
# throw error if number of categories not correct
if(!all(apply(dat, 2, function(x) length(unique(x))) == categories))
stop('Number of categories generated is incorrect')
dat
}
Analyse <- function(condition, dat, fixed_objects = NULL) {
Attach(condition)
syntax <- genLavaanSyntax(factors=factors, indicators=indicators, analyse=TRUE)
mod <- cfa(syntax, dat, ordered=colnames(dat), estimator=estimator)
# check that model and coefficients are reasonable
if(!lavInspect(mod, 'converged')) stop('Model did not converge')
pick_lambdas <- matrix(TRUE, indicators*factors, factors)
if(factors == 2)
pick_lambdas[(indicators+1):(indicators*3)] <- FALSE
cfs <- lavInspect(mod, what="std")$lambda[pick_lambdas]
if(any(cfs > 1 | cfs < -1))
stop('Model contains Heywood cases')
if(factors > 2 && abs(lavInspect(mod, what="std")$psi[2,1]) >= 1)
stop('Latent variable psi matrix not positive definite')
# extract desired results
fit <- fitMeasures(mod)
ses <- lavInspect(mod, what="se")$lambda[pick_lambdas]
fitstats <- if(estimator == 'WLS') fit[c('chisq', 'df', 'pvalue')]
else if(estimator == 'WLSMV') fit[c('chisq.scaled', 'df.scaled', 'pvalue.scaled')]
names(fitstats) <- c('chisq', 'df', 'pvalue')
phi21 <- if(factors == 2)
lavInspect(mod, what="std")$psi[1,2] else NULL
ret <- c(fitstats, mean_ses=mean(ses), lambda=cfs, phi21=phi21)
ret
}
Summarise <- function(condition, results, fixed_objects = NULL) {
# model parameters
lambdas <- results[ , grepl('lambda', colnames(results))]
pool_mean_lambdas <- mean(apply(lambdas, 2, mean)) # Equation 10
pool_SD_lambdas <- sqrt(mean(apply(lambdas, 2, var))) # Equation 11
RB_phi21 <- if(condition$factors == 2)
bias(results$phi21, parameter=.3, type='relative', percent=TRUE) else NULL
mean_se <- mean(results$mean_ses)
# goodness-of-fit
edr_05 <- EDR(results$pvalue, alpha = .05)
mean_X2 <- mean(results$chisq)
sd_X2 <- sd(results$chisq)
RB_X2 <- bias(results$chisq, parameter=results$df, type='relative',
percent=TRUE, unname=TRUE)
ret <- c(mean_X2=mean_X2, sd_X2=sd_X2, edr_05=edr_05,
pool_mean_lambdas=pool_mean_lambdas,
pool_SD_lambdas=pool_SD_lambdas, mean_se=mean_se,
RB_X2=RB_X2, RB_phi21=RB_phi21)
ret
}
######################################################
# run simulation
res <- runSimulation(design=Design, replications=500, generate=Generate,
analyse=Analyse, summarise=Summarise, packages='lavaan',
parallel=TRUE, max_errors=100, save_results=TRUE,
filename='FloraCurran2004')
res# A tibble: 300 × 20
N categories skewness_kurtosis factors indicators estimator mean_X2
<dbl> <dbl> <list> <dbl> <dbl> <chr> <dbl>
1 100 2 <dbl [2]> 1 5 WLSMV 4.7961
2 200 2 <dbl [2]> 1 5 WLSMV 4.8579
3 500 2 <dbl [2]> 1 5 WLSMV 5.0052
4 1000 2 <dbl [2]> 1 5 WLSMV 5.0903
5 100 5 <dbl [2]> 1 5 WLSMV 4.8931
6 200 5 <dbl [2]> 1 5 WLSMV 4.9868
7 500 5 <dbl [2]> 1 5 WLSMV 5.1779
8 1000 5 <dbl [2]> 1 5 WLSMV 5.1632
9 100 2 <dbl [2]> 1 5 WLSMV 4.7838
10 200 2 <dbl [2]> 1 5 WLSMV 4.9057
# ℹ 290 more rows
# ℹ 13 more variables: sd_X2 <dbl>, edr_05 <dbl>, pool_mean_lambdas <dbl>,
# pool_SD_lambdas <dbl>, mean_se <dbl>, RB_X2 <dbl>, RB_phi21 <dbl>,
# REPLICATIONS <dbl>, SIM_TIME <chr>, COMPLETED <chr>, SEED <int>,
# ERRORS <int>, WARNINGS <int>RetroSearch is an open source project built by @garambo | Open a GitHub Issue
Search and Browse the WWW like it's 1997 | Search results from DuckDuckGo
HTML:
3.2
| Encoding:
UTF-8
| Version:
0.7.4