1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
|
library(OpenMx)
library(testthat)
context("npdWarningWLS")
set.seed(1)
#
# Step 2: set up simulation parameters
# Note: nVariables>=3, nThresholds>=1, nSubjects>=nVariables*nThresholds (maybe more)
# and model should be identified
#
nVariables<-5
nFactors<-1
nThresholds<-3
nSubjects<-500
isIdentified<-function(nVariables,nFactors) as.logical(1+sign((nVariables*(nVariables-1)/2) - nVariables*nFactors + nFactors*(nFactors-1)/2))
# if this function returns FALSE then model is not identified, otherwise it is.
isIdentified(nVariables,nFactors)
loadings <- matrix(.9,nrow=nVariables,ncol=nFactors)
residuals <- 1 - (loadings * loadings)
sigma <- loadings %*% t(loadings) + vec2diag(residuals)
mu <- matrix(0,nrow=nVariables,ncol=1)
# Step 3: simulate multivariate normal data
set.seed(1234)
continuousData <- mvtnorm::rmvnorm(n=nSubjects,mu,sigma)
# Step 4: chop continuous variables into ordinal data
# with nThresholds+1 approximately equal categories, based on 1st variable
quants<-quantile(continuousData[,1], probs = c((1:nThresholds)/(nThresholds+1)))
ordinalData<-matrix(0,nrow=nSubjects,ncol=nVariables)
for(i in 1:nVariables)
{
ordinalData[,i] <- cut(as.vector(continuousData[,i]),c(-Inf,quants,Inf))
}
# Step 5: make the ordinal variables into R factors
ordinalData <- mxFactor(as.data.frame(ordinalData),levels=c(1:(nThresholds+1)))
dualData <- cbind(continuousData,ordinalData)
colnames(dualData) <- paste0('v', 1:ncol(dualData))
mxd <- omxCheckWarning(omxAugmentDataWithWLSSummary(mxData(dualData, 'raw'), type="ULS", fullWeight=FALSE,
allContinuousMethod="marginals"),
"fake.data: marginal covariance matrix is non-positive definite")
cov <- mxd$observedStats$cov
omxCheckTrue(any(eigen(cov)$val < 0))
|
