Hello, I often got NA in confidence interval when fitting univariate or bivariate SEM and changing optimizers didn't help. Here is the code of a bivariate ACE model, and the output when I use 'SLSQP' optimizer . I wonder how can I aviod NA when estimating CI of parameters and algebra. Thank you so much!
vars <-c('tea_5g','smk_num_4g') covars <- c("age","region_type") nv <- 2 # number of variables per twin ncv <- 2 # number of covariates ntv <- nv*2 # number of variables per pair nth <- 3 # number of max thresholds = the number of category - 1 selVars <-paste(vars,c(rep(1,nv),rep(2,nv)),sep="") #c('tea2g1', 'smk_num_4g1', 'tea2g2', 'smk_num_4g2' )# covVars <- paste(covars,c(rep(1,ncv),rep(2,ncv)),sep="") mzData <- subset(tea4g_smk4g_nodrink,zygosity1==0,c(selVars,covVars)) dzData <- subset(tea4g_smk4g_nodrink,zygosity1==1,c(selVars,covVars)) sum(is.na(mzData)) sum(is.na(dzData)) mzData <- na.omit(mzData) dzData <- na.omit(dzData) dim(mzData);dim(dzData) # Generate Descriptive Statistics mzData[,selVars] <-mxFactor(mzData[,selVars], levels=c(0:nth) ) dzData[,selVars] <-mxFactor(dzData[,selVars], levels=c(0:nth) ) # Set Starting Values svb <- 0.1 # sv for beta,get from tryhardordinal svLTh <- -1.5 # start value for first threshold svITh <- 1 # start value for increments svTh <- matrix(rep(c(svLTh,(rep(svITh,nth-1)))),nrow=nth,ncol=nv) # start value for thresholds lbTh <- matrix(rep(c(-3,(rep(0.001,nth-1))),nv),nrow=nth,ncol=nv) # lower bounds for thresholds svCor <- .5 # start value for correlations # Create Labels labTh <- function(lab,vars,nth) { paste(paste("t",1:nth,lab,sep=""),rep(vars,each=nth),sep="") } labdata <- cbind(paste("data.",covars[1],c(1,2),sep = ""),paste("data.",covars[2],c(1,2),sep = "")) #c("data.age1","data.age2","data.region_type1","data.region_type2") labeta <-cbind(paste("B",covars,"_",vars[1],sep = ""),paste("B",covars,"_",vars[2],sep = "")) lathmz <- c(labTh("MZ","tea_5g",nth),labTh("Z","smk_num_4g",nth),labTh("MZ","tea_5g",nth),labTh("Z","smk_num_4g",nth)) lathdz <- c(labTh("DZ","tea_5g",nth),labTh("Z","smk_num_4g",nth),labTh("DZ","tea_5g",nth),labTh("Z","smk_num_4g",nth)) ### Some parameters included in all "submodels": baseACE <- mxModel('Base', mxMatrix( type="Full", nrow=ncv, ncol=nv, free=T, values=svb,labels=labeta, name="Beta"), mxMatrix( type="Zero", nrow=1, ncol=nv, name="meanG" ), # or type="Zero" # Create Matrices for Variance Components mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=c(2,0.1,2), name="A" ), # label=laLower("A",nv), values=valDiag(svPa,nv),# mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=c(2,0.1,2), name="C" ), mxMatrix( type="Symm", nrow=nv, ncol=nv, free=c(F,T,F), values=c(1, 0.1, 1), name="E" ),# fix E to avoid warning,without warning, E can be freely estimated mxAlgebra( expression=A+C+E, name="V" ), mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I"), mxAlgebra( expression=solve(sqrt(I*V)), name="iSD"), mxAlgebra( expression=solve(sqrt(I*V))%&%V, name ="Rph"), mxAlgebra( expression=solve(sqrt(I*A))%&%A, name ="Ra" ), #cov2cor() mxAlgebra( expression=solve(sqrt(I*C))%&%C, name ="Rc" ), mxAlgebra( expression=solve(sqrt(I*E))%&%E, name ="Re" ), # Algebra to compute standardized variance components mxAlgebra( expression=A/V, name="h2"), mxAlgebra( expression=C/V, name="c2"), mxAlgebra( expression=E/V, name="e2"), # the proportion of ace contributing to the overlap mxAlgebra( expression=sqrt(h2[1,1])*sqrt(h2[2,2])*Ra[1,2]/Rph[1,2],name = "pa"), mxAlgebra( expression=sqrt(c2[1,1])*sqrt(c2[2,2])*Rc[1,2]/Rph[1,2],name = "pc"), mxAlgebra( expression=sqrt(e2[1,1])*sqrt(e2[2,2])*Re[1,2]/Rph[1,2],name = "pe") ) modelMZ <- mxModel( "MZ", mxMatrix( type="Full", nrow=2, ncol=ncv, free=F, label=labdata, name="MZDefVars"), mxAlgebra( expression= cbind(Base.meanG + MZDefVars[1,] %*% Base.Beta,Base.meanG + MZDefVars[2,] %*% Base.Beta), name="expMeanMZ" ), mxMatrix( type="Full", nrow=nth, ncol=ntv, free=T, values=svTh, labels=lathmz, name="ThMZ"), mxMatrix( type="Lower", nrow=nth, ncol=nth, free=FALSE, values=1, name="inc" ), mxAlgebra( expression= inc %*% ThMZ, name="threMZ" ), mxAlgebra( expression= rbind( cbind(Base.A+Base.C+Base.E , Base.A+Base.C), cbind(Base.A+Base.C , Base.A+Base.C+Base.E)), name="expCovMZ" ), mxData( observed=mzData, type="raw" ), mxExpectationNormal( covariance="expCovMZ", means="expMeanMZ", dimnames=selVars, thresholds="threMZ" ), mxFitFunctionML() ) modelDZ <- mxModel( "DZ", mxMatrix( type="Full", nrow=2, ncol=ncv, free=F, label=labdata, name="DZDefVars"), mxAlgebra( expression= cbind(Base.meanG + DZDefVars[1,] %*% Base.Beta,Base.meanG + DZDefVars[2,] %*% Base.Beta), name="expMeanDZ" ), mxMatrix( type="Full", nrow=nth, ncol=ntv, free=T, values=svTh,labels=lathdz, name="ThDZ"), mxMatrix( type="Lower", nrow=nth, ncol=nth, free=FALSE, values=1, name="inc" ), mxAlgebra( expression= inc %*% ThDZ, name="threDZ" ), mxAlgebra( expression= rbind( cbind(Base.A+Base.C+Base.E , 0.5%x%Base.A+Base.C), cbind(0.5%x%Base.A+Base.C , Base.A+Base.C+Base.E)), name="expCovDZ" ), mxData( observed=dzData, type="raw" ), mxExpectationNormal( covariance="expCovDZ", means="expMeanDZ", dimnames=selVars, thresholds="threDZ" ), mxFitFunctionML() ) Conf1 <- mxCI (c ('Base.h2[1,1]','Base.c2[1,1]','Base.e2[1,1]','Base.h2[2,2]','Base.c2[2,2]','Base.e2[2,2]')) Conf2 <- mxCI (c ('Base.Ra[2,1]', 'Base.Rc[2,1]', 'Base.Re[2,1]','Base.Rph[2,1]') ) Conf3 <- mxCI (c ("Base.pa","Base.pc","Base.pe") ) AceModel <- mxModel( "ACE", baseACE,modelMZ, modelDZ, Conf1, Conf2,Conf3, mxFitFunctionMultigroup(c('MZ.fitfunction','DZ.fitfunction'))) # ------------------------------------------------------------------------------ # 4) RUN AceModel AceFit <- mxTryHardOrdinal(AceModel, intervals=T,extraTries = 15,OKstatuscodes=c(0,1,5)) (AceSumm <- summary(AceFit,verbose = T))