Attachment | Size |
---|---|
warning message.png [6] | 36.52 KB |
Dear all,
I have some problems of fitting saturated univariate SEM model of a binary variable "whtr_g". I would appreciate it if anyone could help me.
1. I always get warning message when fitting saturated SEM model(see warning message.png), whether using SLSQP, NPSOL or CSOLNP optimizer. What’ more, the result of three optimizers are different. But I don’t know if there are errors in my code. Please help me to check it.
2. The result of mxCheckIdentification() is “model is not locally identifided”, which means the saturated model and estimated parameters are not reliable, right? In this situation, can I compare nested models(such as equating means of twin1 and twin2, or equating means of MZ and DZ in ACE model) against saturated model in order to test whether means/variances/thresholds could be equated between MZ and DZ twins. The code is “mxCompare(SatFit,SubmodelFit)”.
My code are as follows:
selVars <- c('whtr_g1','whtr_g2') covVars <- c("age1","age2","region_type1","region_type2") mzData <- subset(adipodata, zygosity1==0, c(selVars,covVars)) dzData <- subset(adipodata, zygosity1==1, c(selVars,covVars)) mMZ <- mean(unlist(mzData[,1:2]),na.rm=TRUE) #start value for mean mDZ <- mean(unlist(dzData[,1:2]),na.rm=TRUE) #start value for mean cMZ <- cov(mzData[,1:2] , use='complete.obs') #start value for variance covariance cDZ <- cov(dzData[,1:2] , use='complete.obs') #start value for variance covariance nv <- 1 # number of variables ntv <- nv*2 # number of variables per pair nth <- 1 # number of max thresholds laMeMZT <- c("m1MZT","m2MZT") # labels for means for MZT twins laMeDZT <- c("m1DZT","m2DZT") # labels for means for DZT twins laVaMZT <- c("MZ11","MZ12","MZ12","MZ22") # labels for (co)variances for MZT twins laVaDZT <- c("DZ11","DZ12","DZ12","DZ22") # labels for (co)variances for DZT twins # tell mx the data is ordinal mzData$whtr_g1 <- mxFactor(mzData$whtr_g1, levels= c(0:nth)) mzData$whtr_g2 <- mxFactor(mzData$whtr_g2, levels= c(0:nth)) dzData$whtr_g1 <- mxFactor(dzData$whtr_g1, levels= c(0:nth)) dzData$whtr_g2 <- mxFactor(dzData$whtr_g2, levels= c(0:nth)) modelMZ<- mxModel("MZ", # Define definition variables mxMatrix( type="Full", nrow=2, ncol=2, free=F, label=c("data.age1","data.region_type1","data.age2","data.region_type2"), name="MZDefVars"), # define beta mxMatrix( type="Full", nrow=1, ncol=2, free=TRUE, values=0,labels=c('BaTHage','BaTHregion_type'), name="BageTH"), # expected mean Matrices mxMatrix( type="Full", nrow=1, ncol=ntv, free=T, values=mMZ, labels=laMeMZT, name="meanMZT" ), #two labels# mxAlgebra( expression= meanMZT + BageTH%*%MZDefVars, name="expMeanMZT" ), # Define threshold matrices : nth--number of thresholds mxMatrix( type="Full", nrow=nth, ncol=ntv, free=T, values=0.22, lbound=-3, ubound=3,labels=c("Tmz1","Tmz2"), name="expThresMZ"), # expected variance–covariance Matrices mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=T, values=cMZ,labels=laVaMZT,name="expCovMZ"), # Read in the data mxData( observed=mzData, type="raw" ), mxExpectationNormal( covariance="expCovMZ", means="expMeanMZT", dimnames=selVars, thresholds="expThresMZ"), mxFitFunctionML() ) modelDZ<- mxModel("DZ", mxMatrix( type="Full", nrow=2, ncol=2, free=F, label=c("data.age1","data.region_type1","data.age2","data.region_type2"), name="DZDefVars"), mxMatrix( type="Full", nrow=1, ncol=2, free=TRUE, values=0,labels=c('BaTHage','BaTHregion_type'), name="BageTH"), mxMatrix( type="Full", nrow=1, ncol=ntv, free=T, values=mDZ, labels=laMeDZT, name="meanDZT" ), #two labels# mxAlgebra( expression= meanDZT + BageTH%*%DZDefVars, name="expMeanDZT" ), mxMatrix( type="Full", nrow=nth, ncol=ntv, free=T, values=0.20, lbound=-3, ubound=3,labels=c("Tdz1","Tdz2"), name="expThresDZ"), mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=T, values=cDZ,labels=laVaDZT,name="expCovDZ"), mxData( observed=dzData, type="raw"), mxExpectationNormal( covariance="expCovDZ", means="expMeanDZT", dimnames=selVars, thresholds="expThresDZ" ), mxFitFunctionML() ) obj <- mxFitFunctionMultigroup(c('MZ.fitfunction','DZ.fitfunction')) Conf <- mxCI (c ('MZ.expCovMZ[1,1]', 'MZ.expCovMZ[2,1]','DZ.expCovDZ[2,1]', 'MZ.BageTH[1,1]' )) SatModel <- mxModel( "Sat", modelMZ, modelDZ, obj, Conf ) SatFit <- mxRun(SatModel, intervals=TRUE)
Many Thanks!!!!
OpenMx version: 2.12.2
R version: R version 3.3.3