Attachment | Size |
---|---|
Model.pdf [6] | 45.14 KB |
Hello, I want to estimate a modified CP model, where I have a Cholesky model with three latent traits (see picture attached).
The d- and e-variables are ordinal (d = 2 categories; e = 3 categories). The s-variable is continuous and in order to integrate it into the LISREL-matrices I define a single factor measurement model with the residual variance set to 0 and the factor loading to 1.
For the ordinal variables I fix the means to zero and the residual variances to 1. I followed the approach proposed here [7] to fix the residual variances of the categorical variables to one instead of the variance of the underlying continuous variable. So I can set the thresholds free.
Do I need an additional constraint for the second order factors?
I get the following error:
All fit attempts resulted in errors - check starting values or model specification
since I have tried to estimate the model with different starting values and I am using mxTryHardOrdinal() I suppose that the specification and not the starting values is the cause of the error.
I used WLS because of the high number of categorical variables: I tried the estimation with FIML but after waiting a long time without results I stopped the process. Since the s-variable has a relative high number of NAs, this may be an issue?
Here is my code:
rm(list=ls()) library(OpenMx) library(dplyr) # Import Data data <- read.csv(file = "data_wide.csv", header = TRUE) # Number of variables m <- 6 # lat. end. V. n <- 18 # lat. ex. V. p <- 32 # Ind. of lat. end. V. q <- 0 # Ind. of lat. ex. V. # Vectors of variables enLV <- c("E","D","S") mlab <- c(paste0(enLV,1),paste0(enLV,2)) exLV <- c("A","C","E") nlab <- c(paste0(exLV[1],enLV,1),paste0(exLV[2],enLV,1),paste0(exLV[3],enLV,1), paste0(exLV[1],enLV,2),paste0(exLV[2],enLV,2),paste0(exLV[3],enLV,2)) plab <- c(paste0("e",1:5,"t",1),paste0("d",1:10,"t",1),"st1",paste0("e",1:5,"t",2),paste0("d",1:10,"t",2),"st2") del <- c(plab[c(6:15,22:31)]) ext <- c(plab[c(1:5,17:21)]) mzData <- subset(data, zyg==1, plab) dzData <- subset(data, zyg==2, plab) # define categorical items as mxfactors # 1. externalizing items (3 levels) mzData[ext] <- mxFactor(mzData[ext], levels = 1:3) dzData[ext] <- mxFactor(dzData[ext], levels = 1:3) # 2. delinquency items (2 levels) mzData[del] <- mxFactor(mzData[del], levels = 0:1) dzData[del] <- mxFactor(dzData[del], levels = 0:1) # Define the expected covariance matrix using hte LISREL approach pathA <- mxMatrix(type = "Lower", nrow = 3, ncol = 3, byrow = TRUE, free = c(TRUE,FALSE,FALSE, TRUE,TRUE,FALSE, TRUE,TRUE,TRUE), values = c(.5,0,0, 0,.5,0, 0,0,.5), labels = c("a11",NA,NA, "a21","a22",NA, "a31","a32","a33"), lbound = c(.0001,NA,NA, -10,.0001,NA, -10,-10,.0001), name = "A") pathC <- mxMatrix(type = "Lower", nrow = 3, ncol = 3, byrow = TRUE, free = c(TRUE,FALSE,FALSE, TRUE,TRUE,FALSE, TRUE,TRUE,TRUE), values = c(.5,0,0, 0,.5,0, 0,0,.5), labels = c("c11",NA,NA, "c21","a22",NA, "c31","c32","c33"), lbound = c(.0001,NA,NA, -10,.0001,NA, -10,-10,.0001), name = "C") pathE <- mxMatrix(type = "Lower", nrow = 3, ncol = 3, byrow = TRUE, free = c(TRUE,FALSE,FALSE, TRUE,TRUE,FALSE, TRUE,TRUE,TRUE), values = c(1,0,0, 0,1,0, 0,0,1), labels = c("e11",NA,NA, "e21","e22",NA, "e31","e32","e33"), lbound = c(.0001,NA,NA, -10,.0001,NA, -10,-10,.0001), name = "E") Null <- mxMatrix(type = "Zero", nrow = 3, ncol = 9, name = "Zeros") Gamma <- mxAlgebra(rbind(cbind(A,C,E,Zeros), cbind(Zeros,A,C,E)), name = "GA") Beta <- mxMatrix(type = "Zero", nrow = m, ncol = m, name = "B") Id <- mxMatrix(type = "Iden", nrow = m, ncol = m, name = "I") LaExt <- mxMatrix(type = "Full", nrow = 5, ncol = 1, byrow = TRUE, free = c(FALSE,TRUE,TRUE,TRUE,TRUE), values = c(1,.8,.8,.8,.8), labels = c("la_ex1","la_ex2","la_ex3","la_ex4","la_ex5"), name = "LExt") NullExt <- mxMatrix(type = "Zero", nrow = 16-5, ncol = 1, name = "NExt") LaDel <- mxMatrix(type = "Full", nrow = 11, ncol = 1, byrow = TRUE, free = c(FALSE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,FALSE), values = c(1,.8,.8,.8,.8,.8,.8,.8,.8,.8,0), labels = c("la_de1","la_de2","la_de3","la_de4","la_de5","la_de6","la_de7","la_de8","la_de9","la_de10", NA), name = "LDel") NullDel <- mxMatrix(type = "Zero", nrow = 5, ncol = 1, name = "NDel") LaSt <- mxMatrix(type = "Full", nrow = 1, ncol = 1, free = FALSE, values = 1, labels = "la_st", name = "LSt") NullSt <- mxMatrix(type = "Zero", nrow = 16-1, ncol = 1, name = "NSt") NullFull <- mxMatrix(type = "Zero", nrow = 16, ncol = 1, name = "NFull") LaY <- mxAlgebra(expression = cbind(rbind(LExt, NExt, NFull), rbind(NDel, LDel, NFull), rbind(NSt, LSt, NFull), rbind(NFull, LExt, NExt), rbind(NFull, NDel, LDel), rbind(NFull, NSt, LSt)), name = "LY") errorval <- rep(c(rep(0,15),0),2) errorlab <- rep(c("eex1","eex2","eex3","eex4","eex5", "ed1","ed2","ed3","ed4","ed5","ed6","ed7","ed8","ed9","ed10", "es"),2) Thee <- mxMatrix(type = "Diag", nrow = 32, ncol = 32, byrow = FALSE, free = FALSE, values = errorval, lbound = 0, labels = errorlab, name = "THE") Psi <- mxMatrix(type = "Diag", nrow = 6, ncol = 6, free = FALSE, values = 0, lbound = 0, labels = rep(c("eEx","eD","eS"),2), name = "PS") PhiVar <- mxMatrix(type = "Iden", ncol = 9, nrow = 9, name = "PHV") PhiCMZ <- mxMatrix(type = "Diag", ncol = 9, nrow = 9, values = c(1,1,1,1,1,1,0,0,0), name = "PHCMZ") PhiCDZ <- mxMatrix(type = "Diag", ncol = 9, nrow = 9, values = c(.5,.5,.5,1,1,1,0,0,0), name = "PHCDZ") PhiMZ <- mxAlgebra(rbind(cbind(PHV,PHCMZ), cbind(PHCMZ,PHV)), name = "PHMZ") PhiDZ <- mxAlgebra(rbind(cbind(PHV,PHCDZ), cbind(PHCDZ,PHV)), name = "PHDZ") eCovMZ <- mxAlgebra(expression= LY%*%solve(I-B)%*%(GA%*%PHMZ%*%t(GA)+PS)%*%t(solve(I-B))%*%t(LY)+THE, name="expCovMZ") # Formula: Bollen (1989) eCovDZ <- mxAlgebra(expression= LY%*%solve(I-B)%*%(GA%*%PHDZ%*%t(GA)+PS)%*%t(solve(I-B))%*%t(LY)+THE, name="expCovDZ") #eCovMZ <- mxAlgebra(expression= LY%*%I%*%(GA%*%PHMZ%*%t(GA)+PS)%*%I%*%t(LY)+THE, name="expCovMZ") #eCovDZ <- mxAlgebra(expression= LY%*%I%*%(GA%*%PHDZ%*%t(GA)+PS)%*%I%*%t(LY)+THE, name="expCovDZ") # Defining the mean matrix # mean of isei mean <- mxMatrix(type = "Full", nrow = 1, ncol = 32, free = rep(c(rep(FALSE,15),TRUE),2), values = rep(c(rep(0,15),40),2), labels = rep(c(rep(NA,15),"mst"),2), name = "expMean") # Defining the threshold matrix # Example of logic of premultiplication of thinG with inc ! -> Ensures ordering of thresholds for items with >1 thresholds #start <- matrix(c(rep(-1.5,30),rep(1,30)),nrow = 2, ncol = 30, byrow = TRUE) #one <- matrix(c(1,0,1,1), nrow = 2, ncol = 2, byrow = TRUE) #one%*%start thDim <- plab[c(1:5, 17:21, 6:15, 22:31)] thLab <- c(rep(paste0("th1","e",1:5),2),rep(paste0("th1","d",1:10),2),rep(paste0("th2","e",1:5),2),rep(NA,20)) thFree <- c(rep(TRUE,40),rep(FALSE,20)) thVal <- c(rep(-1.5,30),rep(1,10),rep(NA,20)) thBound <- c(rep(-3,30),rep(0.001,10),rep(NA,20)) thinG <- mxMatrix( type="Full", nrow=2, ncol=30, byrow = TRUE, free=thFree, values=thVal, lbound=thBound, labels=thLab, name="thinG" ) # Incremental matrix: 1 row = start values of 1. threshold; next rows = increments departing from upper row (lbound necessary so that there is at least a 0.0001 increase! nad no decrease to insure the order!) inc <- mxMatrix( type="Lower", nrow=2, ncol=2, free=FALSE, values=1, name="inc" ) threG <- mxAlgebra( expression= inc %*% thinG, name="expTh" ) # Define data object dataMZ <- mxData( observed=mzData, type="raw" ) dataDZ <- mxData( observed=dzData, type="raw" ) # Define expectation objects expMZ <- mxExpectationNormal(covariance="expCovMZ", means="expMean", thresholds="expTh", dimnames=plab, threshnames = thDim) expDZ <- mxExpectationNormal(covariance="expCovDZ", means="expMean", thresholds="expTh", dimnames=plab, threshnames = thDim) # Fit function (WLS) fitfun <- mxFitFunctionWLS() # parameters pars <- list(pathA, pathC, pathE, Null, Gamma, Beta, Id, LaExt, NullExt, LaDel, NullDel, LaSt, NullSt, NullFull, LaY, Thee, Psi, mean, thinG, inc, threG, PhiVar) # group specific model objects modelMZ <- mxModel(pars, PhiMZ, PhiCMZ, eCovMZ, dataMZ, expMZ, fitfun, name="MZ") modelDZ <- mxModel(pars, PhiDZ, PhiCDZ, eCovDZ, dataDZ, expDZ, fitfun, name="DZ") multi <- mxFitFunctionMultigroup( c("MZ","DZ") ) # overall model object modelACE <- mxModel( "ACE", pars, modelMZ, modelDZ, multi) # run model mxOption( NULL , 'Default optimizer' , 'SLSQP' ) set.seed(1) fitACE <- mxTryHardOrdinal(modelACE, extraTries = 10, exhaustive = TRUE) sumACE <- summary( fitACE ) sumACE
Here is the output:
Summary of ACE Starting values are not feasible. Consider mxTryHard() free parameters: name matrix row col Estimate lbound ubound 1 a11 A 1 1 0.66719606 1e-04 2 a21 A 2 1 0.24639404 -10 3 a31 A 3 1 0.01989468 -10 4 a22 A 2 2 0.58462127 1e-04 5 a32 A 3 2 0.23233684 -10 6 a33 A 3 3 0.20895600 1e-04 7 c11 C 1 1 0.38163935 1e-04 8 c21 C 2 1 -0.23787426 -10 9 c31 C 3 1 0.05621162 -10 10 c32 C 3 2 0.09830437 -10 11 c33 C 3 3 0.65105809 1e-04 12 e11 E 1 1 0.78303864 1e-04 13 e21 E 2 1 -0.17010759 -10 14 e31 E 3 1 0.24441680 -10 15 e22 E 2 2 0.81938448 1e-04 16 e32 E 3 2 -0.15604409 -10 17 e33 E 3 3 1.10636594 1e-04 18 la_ex2 LExt 2 1 0.92273536 19 la_ex3 LExt 3 1 0.88631576 20 la_ex4 LExt 4 1 1.03505741 21 la_ex5 LExt 5 1 0.95816739 22 la_de2 LDel 2 1 0.54168182 23 la_de3 LDel 3 1 0.65248564 24 la_de4 LDel 4 1 1.02872796 25 la_de5 LDel 5 1 0.47255096 26 la_de6 LDel 6 1 0.58154706 27 la_de7 LDel 7 1 0.92095349 28 la_de8 LDel 8 1 0.83431902 29 la_de9 LDel 9 1 0.75920953 30 la_de10 LDel 10 1 0.74839136 31 mst expMean 1 16 48.11579465 32 th1e1 thinG 1 1 -1.72455967 -3 33 th2e1 thinG 2 1 0.88925538 0.001 34 th1e2 thinG 1 2 -1.31678786 -3 35 th2e2 thinG 2 2 1.29860129 0.001 36 th1e3 thinG 1 3 -1.97153396 -3 37 th2e3 thinG 2 3 0.81570670 0.001 38 th1e4 thinG 1 4 -1.31957163 -3 39 th2e4 thinG 2 4 1.23111741 0.001 40 th1e5 thinG 1 5 -1.09611506 -3 41 th2e5 thinG 2 5 1.07165598 0.001 42 th1d1 thinG 1 11 -0.91010120 -3 43 th1d2 thinG 1 12 -1.38949672 -3 44 th1d3 thinG 1 13 -1.66416074 -3 45 th1d4 thinG 1 14 -0.97191166 -3 46 th1d5 thinG 1 15 -1.64255135 -3 47 th1d6 thinG 1 16 -1.15400035 -3 48 th1d7 thinG 1 17 -1.11629872 -3 49 th1d8 thinG 1 18 -1.49572992 -3 50 th1d9 thinG 1 19 -1.97838685 -3 51 th1d10 thinG 1 20 -1.98501499 -3 Model Statistics: | Parameters | Degrees of Freedom | Fit ( units) Model: 51 62399 NA Saturated: NA NA NA Independence: NA NA NA Number of observations/statistics: 2042/62450 timestamp: 2020-11-13 08:57:16 Wall clock time: 2.551282 secs OpenMx version number: 2.18.1 Need help? See help(mxSummary)