saturated univariate SEM model of a binary variable
diana
Joined: 04/17/2019
| Attachment | Size |
|---|---|
| warning message.png | 36.52 KB |
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
how to calculate tetrachoric correlation using saturated model
Thanks a lot!
Log in or register to post comments
In reply to how to calculate tetrachoric correlation using saturated model by diana
identify variance
Having looked at your script again, I see another reason why your model is unidentified. You also need to fix the variance of the (latent liability underlying the) binary trait. Change both model-expected covariance matrices to be "standardized" matrices instead:
mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=cMZ,labels=laVaMZT,name="expCovMZ"),
#and
mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=cDZ,labels=laVaDZT,name="expCovDZ"),
The off-diagonal elements of those two matrices will contain the MZ and DZ tetrachoric correlations, respectively.
Log in or register to post comments
model identification
I've never seen that warning message before. It looks like it's being raised by R, not by any OpenMx code. If it's just telling you that some rownames aren't being preserved, then it's probably benign.
I can tell at a glance that your model is unidentified. You need to fix either the regression intercepts 'meanMZT' and 'meanDZT', or fix the thresholds 'expThresMZ' and 'expThresDZ'. Whichever pair you fix, it's simplest to fix it to zero.
The unidentification could explain why you're getting status RED, and why the optimizers don't agree.
Log in or register to post comments
In reply to model identification by AdminRobK
different results yielded by different functions
But I find that when fitting models, different functions will return different results. Because the model contains binary variables, I have tried
mxRun(), mxTryHard(), mxTryHard()with argumentfinetuneGradient=FALSEandmxTryHardOrdinal(). And it seems that when I usemxTryHardOrdinal(),the optimizer can not reach an acceptable solution. Because the function don’t stop untill reach the maximum number of attempts I set. The different rusults yielded by different functions were showed in picture uploaded. In this case, I don’t know which result is reliable and which function I should choose.In regard with tetrachoric correlations, I’m sorry that I’m not sure if I understand what you mean correctly.
1. My understanding is that
type="symm"andtype="stand"are both right for model-expected covariance matrices. If I want to get tetrachoric correlations, I can choosetype="stand", otherwise the type of expected covariance matrices can also be"symm",is that right?2. What’s more, if I set
type="symm", the lower bound for variances must be defined aslbVas <- matrix(NA,ntv,ntv),diag(lbVas) <- 0.0001, right? Because I find the estimated parameters will changed if I don’t uselbVas.Many thanks!
Log in or register to post comments
In reply to model identification by AdminRobK
and I subsequently found that
mxTryHard()with argumentfinetuneGradient=FALSE, extraTries=41, and the function stoped before reach the maximum number of attempts I set.)Look forward to your replay. Many thanks!
Log in or register to post comments
mxTryHard*, tetrachorics
By default, `mxTryHardOrdinal()` uses `exhaustive=TRUE`, meaning that it uses all of its extra attempts, and returns the best solution it finds.
The default argument values of `mxTryHardOrdinal()` are tailored toward analysis of threshold data, so it's what I'd recommend. By default, it exhausts all of its extra tries, and returns the solution with the smallest fitfunction value, even if it has status Red.
If you want to get the same result from `mxTryHardOrdinal()` whenever you run your script, set the random-number-generator seed with `set.seed()` before the call to `mxTryHardOrdinal()`. Alternately, after you run `mxTryHardOrdinal()`, use the start values it used to find the solution as the actual start values in your script (see `omxSetParameters()`), and replace the call to `mxTryHardOrdinal()` with a call to `mxRun()`.
Let me explain something of which you might not be aware. When you analyze a threshold variable in OpenMx, the variable is modeled as a discretized indicator of an underlying, normally-distributed continuum (often called the "liability"). The model-expected mean and variance of the variable are the mean and variance of that underlying normal distribution. If the variable only has two categories (one threshold), then in order to identify the model, its variance must be fixed, and one of either its regression intercept or its threshold must also be fixed. In your script, you're fixing the intercept. But the variance must also be fixed to identify the model. So, if you use `type="Symm"` for the model-expected covariance matrices, the diagonal elements must be fixed to the same value (which is arbitrary, but 1.0 is the simplest choice). If you instead use `type="Stand"`, the diagonal elements will be automatically fixed to 1.0. If the diagonal elements are fixed to 1.0 one way or another, then the off-diagonals will be tetrachoric correlations.
If the model is unidentified, then there are infinitely many sets of parameter values that will fit the data equally well.
Log in or register to post comments
In reply to mxTryHard*, tetrachorics by AdminRobK
Which is the most saturated model for binary variable
If I use
type="Symm"for the model-expected covariance matrices and fix the diagonal elements to the same value, then I can get the same result each time I run my script. However, it means that I constrain the variances of twin1 and twin2 to be equal, and the model isn't the most saturated model. In fact, firstly I want to fit the most saturated model, and then fit the constrained saturated models(equate the thresholds and variances asross twin order and zygosity step by step), finally report tetrachoric correlations from the best-fitting saturated model. Therefore, I added these lines to my script and it seemed work.#in the most saturated model
mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=T, values=cDZ,labels=laVaDZT,lbound=lbVas,name="expCovDZ"),
mxAlgebra(expCovDZ[2,1]/sqrt(expCovDZ[1,1]%*%expCovDZ[2,2]),name="CorDZ")
#Sub3Model: equate var across Twin 1 and Twin 2
Sub3Model<- omxSetParameters(Sub2Model, labels="MZ22",newlabels="MZ11", free=T, values=SatFit$MZ.expCovMZ[1,1]$values)
Sub3Model<- omxSetParameters(Sub3Model, labels="DZ22",newlabels="DZ11", free=T, values=SatFit$DZ.expCovDZ[1,1]$values)
#Sub4Model: equate var across zygosity
Sub4Model<- omxSetParameters(Sub3Model, labels="DZ11",newlabels="MZ11", free=T, values=SatFit$MZ.expCovMZ[1,1]$values)
I am not sure if my thought is right ,or for bianry variable, the most saturated model is the model constraining variances across twin order to be equal and then I just need to constrain thresholds.
Thank you very much for your help!
Log in or register to post comments
fix variances for identification
For a binary variable, the variances (i.e., the diagonal elements) need to be fixed to some value. Otherwise, the model will be unidentified. With only two categories, there's not enough information to freely estimate the variances.
Log in or register to post comments
Ok, now I know how to modify
Thank you very much!
Log in or register to post comments