Moderation model produces weird results

Posted on
Picture of user. Julia Joined: 03/29/2012
Hi all,

I am running a bivariate moderation model where the phenotype of interest if dichotomous and the moderator is continuous (centered and scaled, but highly skewed). The estimates that I am getting are way off and don't match to the analysis without the moderation, this is in particular true for the correlation estimates. That is, when I compare estimates at the zero level of moderator with the estimates when there is no moderation in the model, they are very different. For example, rE=-0.27 at M=0 and rE=0.06 if moderation is not modelled; rP=-0.21 at M=0 and rP=0.17 if no moderation (also 0.17 observed from the data).
I have three moderators that we are interested in and such discrepancies are observed for two of them.

What could be a reason for that? I'm a bit lost and don't know how to proceed further, whether I can trust the results of the moderation model or not.

A bit of background: I am running ADE model based on the previous results, but CI for D are very wide and D could be dropped out of the model without significant deterioration of the fit. We decided to keep it in the model for now due to the reviews we got for our previous results (since CI's are large).
For one of the moderators, these estimates are off if D is present, but become consistent with the previous analysis if D is dropped. For the other two moderators, dropping D did not change rE and rP values.
Variance estimates for the main phenotype are consistent and in agreement with the previous analysis if D is dropped out of the model.
I did try both Purcell's bivariate moderation model (based on the Cholesky decomposition) and correlated factor solution with moderation of the paths, and they both consistently give weird estimates for the two of three moderators.
Also, there seem to be no significant moderation of any of the paths.

If anyone has any insight and idea what is going on there and why estimates at M=0 don't match with the estimates when there is no moderation, I would be very grateful!!!

Thank you in advance!
Julia

Replied on Tue, 04/16/2019 - 15:34
Picture of user. AdminRobK Joined: 01/24/2014

I wouldn't be able to even guess without at least seeing the script you're working with. `summary()` output might help, too.
Replied on Wed, 04/17/2019 - 07:00
Picture of user. Julia Joined: 03/29/2012

Here are the scripts that I use and the output from full moderation models and models without any moderation. When having lowsupport_s as moderator, rE and rP are not consistent across Purcell moderation model and CF moderation model and are different from the main effects model (and from the observed phenotypic correlation). When removing D, they at least become consistent across Purcell and CF models, but still not equal to the no moderation results (not to the observed rP). And only when removing all the moderation, rP estimated is in agreement with rP observed.
Thank you for looking at it!
Replied on Wed, 04/17/2019 - 11:41
Picture of user. AdminRobK Joined: 01/24/2014

In the Purcell model, have you tried removing the lower bounds on the unmoderated terms of the path coefficients? I notice there's an active bound at the solution:

3 aU MZ.a 2 2 -6.908945e-04 0!

It seems to me that the bounds shouldn't matter, since you mean-centered the moderator, and you identify the liability scale by placing a constraint on the variance at the moderator's zero point. Still, I'm not 100% sure about that. Notice that the point estimate above is actually negative, meaning that the bound is slightly violated. If the bound weren't there, would the optimizer try to push that parameter into the negative region?

Note that you don't need to put `var_constraint` into both the MZ and DZ MxModels. It suffices to put it into only one of them.

Have you tried using a different optimizer?

Consider replacing `mxRun()` with mxTryHardOrdinal().

In the correlated-factors model, why is 'Rd' fixed to -0.7?

Replied on Wed, 04/17/2019 - 16:09
Picture of user. Julia Joined: 03/29/2012

In reply to by AdminRobK

I tried to remove lower bounds for the path estimates, and the value went just below zero, but very little:

aU MZ.a 2 2 -0.0004634254

I also tried another optimizer as you suggested, and SLSQP produced nearly identical results, whereas CSOLNP was just hanging for 10 minutes and I had to terminate the session.

As for rD=-0.7 in the CF model, this was based on the best model in terms of AIC when trying different rD values in the range from -1 to 0 with 0.1 step (previous analysis indicated negative rD).

Also thank you for note about putting the variance constraint just once into the model. I actually think that in some other scripts I had it inside the global model. Don't know why I changed it here.

Replied on Wed, 04/17/2019 - 16:37
Picture of user. AdminRobK Joined: 01/24/2014

In reply to by Julia

I tried to remove lower bounds for the path estimates, and the value went just below zero, but very little

Did the fit value appreciably improve?

I also tried another optimizer as you suggested, and SLSQP produced nearly identical results, whereas CSOLNP was just hanging for 10 minutes and I had to terminate the session.

It's encouraging that SLSQP's results agree with NPSOL's. Did you try CSOLNP with the MxConstraint in both the MZ and DZ models? There is a known issue with CSOLNP and redundant equality constraints. In the next OpenMx release, CSOLNP will at least not freeze uninterruptibly when there are redundant equalities.

As for rD=-0.7 in the CF model, this was based on the best model in terms of AIC when trying different rD values in the range from -1 to 0 with 0.1 step (previous analysis indicated negative rD).

But why isn't it a free parameter?

Have you tried `mxTryHardOrdinal()`?

Replied on Thu, 04/18/2019 - 08:34
Picture of user. Julia Joined: 03/29/2012

In reply to by AdminRobK

Have you tried mxTryHardOrdinal()?

Yes, I did. The the results are still the same with no improvement in the fit.

Did the fit value appreciably improve?

The fit value was left unchanged.

Did you try CSOLNP with the MxConstraint in both the MZ and DZ models?

Yes, I put the constraint just into the MZ model and tried all the optimizers. Here are the fit indices:
Cholesky moderation (Purcell)
NPSOL: -2LL=12907.98 AIC= -2100.015
SLSQP: -2LL=12907.98, AIC=-2100.016
CSONLP: -2LL=12907.98, AIC=-2100.016

CF moderation
NPSOL: -2LL =12910.81, AIC = -2105.186
SLSQP: -2LL=12910.84, AIC=-2105.156 (Mx Status Red)
CSONLP: -2LL=12910.84, AIC=-2105.156

But why isn't it a free parameter?

I thought that rA and rD (just like rA and rC) could not be estimated simultaneously, could they?

Replied on Thu, 04/18/2019 - 11:18
Picture of user. AdminRobK Joined: 01/24/2014

In reply to by Julia

Well, it looks as though the optimizers really are finding the solution. I agree that the results seem odd, but I guess there's something wrong with our intuition!

I thought that rA and rD (just like rA and rC) could not be estimated simultaneously, could they?

I'm not used to thinking of moderation in terms of the correlated-factors parameterization, so I could be mistaken here, but I don't see any reason why they couldn't be estimated simultaneously. After all, you were able to estimate a cross-path for _D_, 'dC', in the Cholesky-parameterized model, right? It should be possible to get a correlated-factors solution equivalent to the Cholesky solution. But, it's possible that the correlated-factors parameterization is harder to optimize.

Replied on Thu, 04/18/2019 - 14:37
Picture of user. Julia Joined: 03/29/2012

In reply to by AdminRobK

Yes, what is puzzling me is that the rP estimated are so different from rP observed! And that estimates at M=0 are totally different from estimates with no moderation. What can be the explanation here? Can we trust the results of moderation here.

The reason why we try moderation in terms of CF is the paper by Rathouz PJ et al:

Rathouz PJ, Van Hulle CA, Rodgers JL, Waldman ID, Lahey BB. Specification, testing, and interpretation of gene-by-measured-environment interaction models in the presence of gene-environment correlation. Behav Genet. 2008;38(3):301–315. doi:10.1007/s10519-008-9193-4

There they say that CF model has more power to detect moderation because it has fewer parameters to estimate. Since our power is quite limited due to a low number of twin pairs and due to a not very prevalent dichotomous outcome, we thought to give CF a try. It doesn't seem to provide any evidence of moderation either, but its fit is better, although the correlation estimates are as weird as in the Purcell's model (for two out of three moderators that we tested).

Replied on Sun, 08/18/2019 - 10:14
No user picture. xiyuesenlinyu Joined: 05/24/2019

Hello
I just run the model posted by Julia(bivChol_Moderation.txt), but here are something wrong when there is no moderation in the model:
Error: The job for model 'MainEffects' exited abnormally with the error message: fit is not finite (Ordinal covariance is not positive definite in data 'DZ.data' row 13703 (loc1))
In addition: Warning message:
In model 'MainEffects' Optimizer returned a non-zero status code 10. Starting values are not feasible. Consider mxTryHard()
However, I don't know the exzact reason, could you give any questions.
Thanks!
Replied on Sun, 08/18/2019 - 10:14
No user picture. xiyuesenlinyu Joined: 05/24/2019

Hello
I just run the model posted by Julia(bivChol_Moderation.txt), but here are something wrong when there is no moderation in the model:
Error: The job for model 'MainEffects' exited abnormally with the error message: fit is not finite (Ordinal covariance is not positive definite in data 'DZ.data' row 13703 (loc1))
In addition: Warning message:
In model 'MainEffects' Optimizer returned a non-zero status code 10. Starting values are not feasible. Consider mxTryHard()
However, I don't know the exzact reason, could you give any questions.
Thanks!
Replied on Sun, 08/18/2019 - 16:25
Picture of user. AdminRobK Joined: 01/24/2014

In reply to by xiyuesenlinyu

Evidently, you need better start values for the free parameters. Did you modify the block of syntax that sets the start values? Values that worked for Julia might not work well with your dataset.

Search this website for "start values". You'll find plenty of advice and discussion about the topic, e.g., this thread. Another thing you could try is to replace `mxRun()` with `mxTryHardOrdinal()` in your script, and see if that helps. I can't really offer any more-specific advice without more details from you.

Replied on Mon, 08/19/2019 - 10:00
No user picture. xiyuesenlinyu Joined: 05/24/2019

In reply to by AdminRobK

I still can’t run my model.
Firstly, I run the following syntax, which induce some warning and error messages .

#=======================================================================#
# PREPARE DATA #
#=======================================================================#
# Load Data

#sink("25+ACE-drink_cur.txt")

mxOption(NULL, "Default optimizer", "NPSOL")
#mxOption( NULL, 'Default optimizer' ,'SLSQP' )
# mxOption(NULL, "Default optimizer", "CSOLNP")
# -----------------------------------------------------------------------
# PREPARE DATA

# Read Data
#twinData <- read.table("D:/FHI/Sosiale forhold og helse/Analyses/IBS/SFH_17_09_2018_subset_pss_discordance_strain_lowsupport_family_format.txt", header=T, sep='\t', dec=',')
twinData <-read.dta13("-25-34581pair-coronary-ssex.dta")
twinData <- twinData [!is.na(twinData $coronary1)&!is.na(twinData $coronary2),]
twinData <- twinData [!is.na(twinData $drink_cur1)&!is.na(twinData $drink_cur2),]

str(twinData)
nobs = dim(twinData)[1]
twinData$age1 = scale(twinData$age1)
twinData$age2 = scale(twinData$age2)
#twinData$Bmi1 = scale(twinData$bmi1)
#twinData$Bmi2 = scale(twinData$bmi2)

# Select Ordinal Variables
nth <- 1 # number of thresholds
varso <- c('coronary') # list of ordinal variables names
nvo <- length(varso) # number of ordinal variables
ntvo <- nvo*2 # number of total ordinal variables
ordVars <- paste(varso,c(rep(1,nvo),rep(2,nvo)),sep="")

# Select Variables for Analysis
Vars <- c("drink_cur","coronary")
nv <- length(Vars) # number of variables
nind <- 2
ntv <- nv*nind # number of total variables
selVars <- paste(Vars,c(rep(1,nv),rep(2,nv)),sep="")
def = c("sex", "age")
ndef= length(def)
defVars = paste(def,c(rep(1,ndef),rep(2,ndef)),sep="")

# Select Data for Analysis
mzData <- subset(twinData, zy_sex1==1|2, c(selVars,defVars))
dzData <- subset(twinData, zy_sex1==3|4, c(selVars,defVars))

mzData = mzData[complete.cases(mzData[,c(defVars, "drink_cur1", "drink_cur2")]),]
dzData = dzData[complete.cases(dzData[,c(defVars, "drink_cur1", "drink_cur2")]),]

mzDataF = mzData
dzDataF = dzData
mzDataF[,ordVars] <- mxFactor( x=mzDataF[,ordVars], levels=c(0,1) )
dzDataF[,ordVars] <- mxFactor( x=dzDataF[,ordVars], levels=c(0,1) )

# Raw data in OpenMx format
dataMZ <- mxData(observed = mzDataF, type = "raw" )
dataDZ <- mxData(observed = dzDataF, type = "raw" )

# ---------------------Cholesky part!------------------------------------

# Set up Cholesky ADE decomposition, with RawData and Matrices Input
# -----------------------------------------------------------------------

## Labeling
aLabs <- c("aM","aC","aU")
cLabs <- c("cM","cC","cU")
eLabs <- c("eM","eC","eU")
meanLabs <- c("meanM","meanP")

aModLabs <- c("aMod11","aMod21","aMod22")
cModLabs <- c("cMod11","cMod21","cMod22")
eModLabs <- c("eMod11","eMod21","eMod22")

threshLabs <- paste(varso,"thresh",sep="_")
betaLabs_age <- paste("beta","Age",Vars, sep="_")
betaLabs_sex <- paste("beta","Sex",Vars, sep="_")

# Set Starting Values
frMV <- c(TRUE, FALSE) # free status for variables
frCvD <- diag(frMV,ntv,ntv) # lower bounds for diagonal of covariance matrix
frCvD[lower.tri(frCvD)] <- TRUE # lower bounds for below diagonal elements
frCvD[upper.tri(frCvD)] <- TRUE # lower bounds for above diagonal elements
frCv <- matrix(as.logical(frCvD),4)
svMe <- c(0,0) # start value for means
svPa <- .4 # start value for path coefficient
svPaD <- vech(diag(svPa,nv,nv)) # start values for diagonal of covariance matrix
svPe <- .8 # start value for path coefficient for e
svPeD <- vech(diag(svPe,nv,nv)) # start values for diagonal of covariance matrix
lbPa <- 0 # start value for lower bounds
lbPaD <- diag(lbPa,nv,nv) # lower bounds for diagonal of covariance matrix
lbPaD[lower.tri(lbPaD)] <- -10 # lower bounds for below diagonal elements
lbPaD[upper.tri(lbPaD)] <- NA # lower bounds for above diagonal elements
svTh <- 1.5 # start value for thresholds

pathModVal = c(0,0.1,0.1)
B_AgeVal = 0.5
B_SexVal = 0.5

## Modeling

# Matrices a, c, and e to store a, c, and e Path Coefficients
pathA <- mxMatrix(name = "a", type = "Lower", nrow = nv, ncol = nv, free=T, labels = aLabs, values=svPaD, lbound=lbPaD)
pathC <- mxMatrix(name = "c", type = "Lower", nrow = nv, ncol = nv, free=T, labels = cLabs, values=svPaD, lbound=lbPaD)
pathE <- mxMatrix(name = "e", type = "Lower", nrow = nv, ncol = nv, free=T, labels = eLabs, values=svPeD, lbound=lbPaD)

modPathA = mxMatrix( "Lower", nrow=nv, ncol=nv, free=c(F,T,T), values=pathModVal, labels=aModLabs, name="aMod" )
modPathC = mxMatrix( "Lower", nrow=nv, ncol=nv, free=c(F,T,T), values=pathModVal, labels=cModLabs, name="cMod" )
modPathE = mxMatrix( "Lower", nrow=nv, ncol=nv, free=c(F,T,T), values=pathModVal, labels=eModLabs, name="eMod" )

# Moderator
mod_tw1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels='data.drink_cur1', name="Mod1")
mod_tw2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels='data.drink_cur2', name="Mod2")

# Matrices generated to hold A, D, and E computed Variance Components
varA1 <- mxAlgebra(name = "A1", expression = (a + Mod1%x%aMod) %*% t(a+ Mod1%x%aMod))
varC1 <- mxAlgebra(name = "C1", expression = (c + Mod1%x%cMod) %*% t(c+ Mod1%x%cMod))
varE1 <- mxAlgebra(name = "E1", expression = (e + Mod1%x%eMod) %*% t(e+ Mod1%x%eMod))

covA12 = mxAlgebra(name = "A12", expression = (a + Mod1%x%aMod) %*% t(a+ Mod2%x%aMod))
covC12 = mxAlgebra(name = "C12", expression = (c + Mod1%x%cMod) %*% t(c+ Mod2%x%cMod))
covE12 = mxAlgebra(name = "E12", expression = (e + Mod1%x%eMod) %*% t(e+ Mod2%x%eMod))

covA21 = mxAlgebra(name = "A21", expression = (a + Mod2%x%aMod) %*% t(a+ Mod1%x%aMod))
covC21 = mxAlgebra(name = "C21", expression = (c + Mod2%x%cMod) %*% t(c+ Mod1%x%cMod))
covE21 = mxAlgebra(name = "E21", expression = (e + Mod2%x%eMod) %*% t(e+ Mod1%x%eMod))

varA2 <- mxAlgebra(name = "A2", expression = (a + Mod2%x%aMod) %*% t(a+ Mod2%x%aMod))
varC2 <- mxAlgebra(name = "C2", expression = (c + Mod2%x%cMod) %*% t(c+ Mod2%x%cMod))
varE2 <- mxAlgebra(name = "E2", expression = (e + Mod2%x%eMod) %*% t(e+ Mod2%x%eMod))

myVarA <- mxAlgebra(name = "A0", expression = a %*% t(a))
myVarC <- mxAlgebra(name = "C0", expression = c %*% t(c))
myVarE <- mxAlgebra(name = "E0", expression = e %*% t(e))

# Algebra to compute total variances and standard deviations (diagonal only) per twin
var1 <- mxAlgebra( A1+C1+E1, name="V1" )
var2 <- mxAlgebra( A2+C2+E2, name="V2" )

myVar <- mxAlgebra( A0+C0+E0, name="V0" )

# Constraint on variance of Binary variables
matUnv <- mxMatrix( type="Unit", nrow=nvo, ncol=1, name="Unv1" )
var_constraint <- mxConstraint(V0[2,2] == Unv1, name="Var1")

### Algebra for expected variance/covariance matrices
expCovMZ <- mxAlgebra(name = "expCovMZ", expression = rbind (cbind(A1+C1+E1, A12+C12),
cbind(A21+C21, A2+C2+E2)))

expCovDZ <- mxAlgebra(name = "expCovDZ", expression = rbind (cbind(A1+C1+E1, 0.5%x%A12+C12),
cbind(0.5%x%A21+C21, A2+C2+E2)))

# Matrices for expected Means for females and males
#setting up the regression
intercept <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=frMV, values=svMe, labels=meanLabs, name="intercept" )
threshold <-mxMatrix( type="Full", nrow=1, ncol=ntvo, free=T, values=svTh, labels=threshLabs, name="Threshold" )

# Regression effects
B_Age <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=.1, labels=betaLabs_age, name="bAge" )
defAge <- mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, labels=c("data.drink_cur1","data.drink_cur2"), name="Age")

B_Sex <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=.1, labels=betaLabs_sex, name="bSex" )
defSex <- mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, labels=c("data.sex1","data.sex2"), name="Sex")

expMean <- mxAlgebra( intercept + (Age %x% bAge) + (Sex %x% bSex) , name="expMean")

inclusions <- list (B_Age, B_Sex, defAge, defSex, intercept, expMean, threshold)

# Objective objects for Multiple Groups
expMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars, thresholds="Threshold", threshnames=ordVars )
expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars ,thresholds="Threshold", threshnames=ordVars )
funML <- mxFitFunctionML()

# MZ and DZ models
modelMZ <- mxModel(pathA, pathC, pathE, modPathA, modPathC, modPathE, mod_tw1, mod_tw2,
varA1, varC1, varE1, covA12, covC12, covE12, covA21, covC21, covE21, varA2, varC2, varE2, var1, var2,
myVarA, myVarC, myVarE, myVar, matUnv, var_constraint,
B_Age, B_Sex, defAge, defSex, intercept, expMean, threshold,
dataMZ, expCovMZ, expMZ, funML, name = "MZ")
modelDZ <- mxModel(pathA, pathC, pathE, modPathA, modPathC, modPathE, mod_tw1, mod_tw2,
varA1, varC1, varE1, covA12, covC12, covE12, covA21, covC21, covE21, varA2, varC2, varE2, var1, var2,
myVarA, myVarC, myVarE, myVar, matUnv, var_constraint,
B_Age, B_Sex, defAge, defSex, intercept, expMean, threshold,
dataDZ, expCovDZ, expDZ, funML, name = "DZ")

#plan <- omxDefaultComputePlan()
#plan$steps$GD <- mxComputeNelderMead()

multi <- mxFitFunctionMultigroup( c("MZ","DZ") )
ACEmodModel <- mxModel( "ACEmod", modelMZ, modelDZ, funML, multi )

#ACEmodFit <- mxAutoStart(ACEmodModel)
ACEmodFit <- mxRun(ACEmodModel)
ACEmodFit <- mxTryHardOrdinal(ACEmodModel)
summary(ACEmodFit)

model <- mxTryHard(ACEmodFit)

mxGetExpected(ACEmodModel, "covariance", "MZ")

MainEffectsModel <- mxModel (ACEmodFit, name='MainEffects')
MainEffectsModel <- omxSetParameters(MainEffectsModel, labels=c('aMod21','aMod22'), values=c(0,0), free=FALSE)
MainEffectsModel <- omxSetParameters(MainEffectsModel, labels=c('cMod21','cMod22'), values=c(0,0), free=FALSE)
MainEffectsModel <- omxSetParameters(MainEffectsModel, labels=c('eMod21','eMod22'), values=c(0,0), free=FALSE)

#ACEmodFit <- mxAutoStart(MainEffectsModel)
#ACEmodFit <- mxRun(MainEffectsModel)
MainEffectsFit <- mxRun(MainEffectsModel)
MainEffectsFit <- mxTryHardOrdinal(MainEffectsModel)
mxCompare(ACEmodFit, MainEffectsFit)
summary(MainEffectsFit)
model <- mxTryHard(MainEffectsFit)

Many of Std.Error are NA in ACEmodModel .

ACEmodFit <- mxRun(ACEmodModel)
Running ACEmod with 22 parameters
Summary of ACEmod
free parameters:
name matrix row col Estimate Std.Error A lbound
1 aM MZ.a 1 1 1.926922e-02 1.902139e-04 ! 0!
2 aC MZ.a 2 1 -4.604157e-01 4.172067e-03 -10
3 aU MZ.a 2 2 1.801983e-02 1.806019e-02 ! 0!
4 cM MZ.c 1 1 1.954307e-02 4.009311e-04 ! 0!
5 cC MZ.c 2 1 -2.463518e-01 1.425125e-02 -10
6 cU MZ.c 2 2 -7.894769e-03 1.938597e-02 ! 0!
7 eM MZ.e 1 1 1.837021e-02 2.447180e-04 ! 0!
8 eC MZ.e 2 1 -8.526084e-01 2.113446e-03 -10
9 eU MZ.e 2 2 -3.436889e-05 1.613466e-02 ! 0!
10 aMod21 MZ.aMod 2 1 -5.222943e+02 NA !
11 aMod22 MZ.aMod 2 2 7.894969e+01 3.810561e+01 !
12 cMod21 MZ.cMod 2 1 -4.150706e+02 NA
13 cMod22 MZ.cMod 2 2 1.734169e+02 NA !
14 eMod21 MZ.eMod 2 1 -1.133761e+03 NA !
15 eMod22 MZ.eMod 2 2 2.499588e+02 NA !
16 beta_Age_drink_cur MZ.bAge 1 1 1.000291e+00 3.464405e-04
17 beta_Age_coronary MZ.bAge 1 2 -1.675897e+03 NA !
18 beta_Sex_drink_cur MZ.bSex 1 1 5.461896e-05 3.377409e-04
19 beta_Sex_coronary MZ.bSex 1 2 -2.344350e-02 1.639405e-02
20 meanM MZ.intercept 1 1 -6.592594e-03 5.280465e-04
21 coronary_thresh MZ.Threshold 1 coronary1 1.178388e+00 2.735566e-02

As for MainEffectsModel(without moderator),

MainEffectsFit <- mxRun(MainEffectsModel)
Running MainEffects with 15 parameters
Error: The job for model 'MainEffects' exited abnormally with the error message: fit is not finite (Ordinal covariance is not positive definite in data 'DZ.data' row 13703 (loc1))
In addition: Warning message:
In model 'MainEffects' Optimizer returned a non-zero status code 10. Starting values are not feasible. Consider mxTryHard()

Secondly, I run mxGetExpected(). However, I just don’t know how to reset the start values.

$MZ
drink_cur1 coronary1 drink_cur2 coronary2
drink_cur1 0.96 0.00 0.32 0.00
coronary1 0.00 0.96 0.00 0.32
drink_cur2 0.32 0.00 0.96 0.00
coronary2 0.00 0.32 0.00 0.96

$DZ
drink_cur1 coronary1 drink_cur2 coronary2
drink_cur1 0.96 0.00 0.24 0.00
coronary1 0.00 0.96 0.00 0.24
drink_cur2 0.24 0.00 0.96 0.00
coronary2 0.00 0.24 0.00 0.96

I also replace mxRun() with MainEffectsFit <- mxTryHardOrdinal(MainEffectsModel). Though the error of NA are removed, the model without moderator didn’t work.

ACEmodFit <- mxTryHardOrdinal(ACEmodModel)
Solution found! Final fit=-448097.41 (started at 171627.54) (11 attempt(s): 8 valid, 3 errors)
Warning message:
In mxTryHard(model = model, greenOK = greenOK, checkHess = checkHess, :
The final iterate satisfies the optimality conditions to the accuracy requested, but the sequence of iterates has not yet converged. Optimizer was terminated because no further improvement could be made in the merit function (Mx status GREEN).

MainEffectsFit <- mxTryHardOrdinal(MainEffectsModel)
All fit attempts resulted in errors - check starting values or model specification.

I read the suggestions you posted, but I failed to use Nelder-Mead implementation. I might be a little stupid, so I need your help.
Please!

Replied on Mon, 08/19/2019 - 11:54
Picture of user. AdminRobK Joined: 01/24/2014

In reply to by xiyuesenlinyu

First off, only put `var_constraint` in one of the MZ or DZ MxModels, not both. That won't matter if you're using NPSOL, but it will be a problem if you use CSOLNP (or Nelder-Mead).

Many of Std.Error are NA in ACEmodModel .

It looks like the optimizer is reaching a solution where the Hessian (as calculated) isn't positive-definite (status code 5). Your phenotype is a threshold trait, and due to the limited accuracy of the algorithm for the multivariate-normal probability integral, sometimes code 5 can occur even when the optimizer has found a minimum. Therefore, you'll want to find the solution with the smallest fitfunction value, even if it has status code 5. Try requesting more attempts from `mxTryHardOrdinal()` via argument `extraTries`, e.g. `extraTries=30`.

As for MainEffectsModel(without moderator),

I suggest running the main-effects model *before* the moderation model. Use `free=FALSE` when creating `modPathA`, `modPathC`, and `modPathE`. Then, the first MxModel you run will be the main effects model. Then, create the moderation model from the fitted main-effects model, and use omxSetParameters() to free the moderation parameters.

Secondly, I run mxGetExpected(). However, I just don’t know how to reset the start values.

What were you trying to do? Use omxSetParameters() to change free parameter values.

I read the suggestions you posted, but I failed to use Nelder-Mead implementation. I might be a little stupid, so I need your help.

If you want to try it, let me suggest some syntax:

plan <- omxDefaultComputePlan()
plan$steps <- list(
NM=mxComputeNelderMead(centerIniSimplex=TRUE),
GD=plan$steps$GD,ND=plan$steps$ND,SE=plan$steps$SE,RD=plan$steps$RD,RE=plan$steps$RD)

Then, put `plan` into the `mxModel()` statement for `MainEffectsModel`, assuming you are creating and running `MainEffectsModel` first. To clear the custom compute plan from an MxModel and go back to the default compute plan, do model@compute <- NULL.

Replied on Tue, 08/20/2019 - 00:43
No user picture. xiyuesenlinyu Joined: 05/24/2019

In reply to by AdminRobK

I change the order of models, and fit the main-effect model first. I also delete the var_constraint in MZ MxModel(of course try this in DZ MxModel). However, there still error massage:

> MainEffectsFit <- mxTryHardOrdinal(MainEffectsModel)

Solution found! Final fit=-444773.78 (started at 169287.04) (11 attempt(s): 5 valid, 6 errors)

It seems the model didn’t run. I want to Try requesting more attempts from mxTryHardOrdinal() via argument extra Tries, e.g. extra Tries=30.

Then, I put plan into the mxModel() statement for MainEffectsModel, errors

plan <- omxDefaultComputePlan()
plan$steps <- list(
NM=mxComputeNelderMead(centerIniSimplex=TRUE),
GD=plan$steps$GD,ND=plan$steps$ND,SE=plan$steps$SE,RD=plan$steps$RD,RE=plan$steps$RD)
multi <- mxFitFunctionMultigroup( c("MZ","DZ") )
#ACEmodModel <- mxModel( "ACEmod", modelMZ, modelDZ, funML, multi
MainEffectsModel<- mxModel( "MainEffects", modelMZ, modelDZ, funML, multi,plan )

MainEffectsFit <- mxTryHardOrdinal(MainEffectsModel)

All fit attempts resulted in errors - check starting values or model specification

As for Try requesting more attempts from mxTryHardOrdinal() via argument extraTries, e.g. extra Tries=30, even 90, errors still.

All fit attempts resulted in errors - check starting values or model specification.

In the main effect model, start values of moderator Path Coefficients are 0, isn’t it? Do I need to change the start values of a, c, and e Path Coefficients.

pathModVal = c(0,0.1,0.1)
B_AgeVal = 0.5
B_SexVal = 0.5

# Matrices a, c, and e to store a, c, and e Path Coefficients
pathA <- mxMatrix(name = "a", type = "Lower", nrow = nv, ncol = nv, free=T, labels = aLabs, values=svPaD, lbound=lbPaD)
pathC <- mxMatrix(name = "c", type = "Lower", nrow = nv, ncol = nv, free=T, labels = cLabs, values=svPaD, lbound=lbPaD)
pathE <- mxMatrix(name = "e", type = "Lower", nrow = nv, ncol = nv, free=T, labels = eLabs, values=svPeD, lbound=lbPaD)

modPathA = mxMatrix( "Lower", nrow=nv, ncol=nv, free=c(F,F,F), values=0, labels=aModLabs, name="aMod" )
modPathC = mxMatrix( "Lower", nrow=nv, ncol=nv, free=c(F,F,F), values=0, labels=cModLabs, name="cMod" )
modPathE = mxMatrix( "Lower", nrow=nv, ncol=nv, free=c(F,F,F), values=0, labels=eModLabs, name="eMod" )

Replied on Tue, 08/20/2019 - 14:27
Picture of user. AdminRobK Joined: 01/24/2014

In reply to by xiyuesenlinyu

First of all, `mxTryHardOrdinal()` uses random-number generation. You can get different results each time you run it, which is especially likely if your model is hard to optimize for whatever reason. If you want reproducible results, precede each call to `mxTryHardOrdinal()` with something to set the random-number generator's seed, e.g. `set.seed(1)`.

I change the order of models, and fit the main-effect model first. I also delete the var_constraint in MZ MxModel(of course try this in DZ MxModel). However, there still error massage:

> MainEffectsFit <- mxTryHardOrdinal(MainEffectsModel)

Solution found! Final fit=-444773.78 (started at 169287.04) (11 attempt(s): 5 valid, 6 errors)

It seems the model didn’t run.

Huh? But it did! `mxTryHardOrdinal()` reported "Solution found!".

Then, I put plan into the mxModel() statement for MainEffectsModel, errors

plan <- omxDefaultComputePlan()
plan$steps <- list(
NM=mxComputeNelderMead(centerIniSimplex=TRUE),
GD=plan$steps$GD,ND=plan$steps$ND,SE=plan$steps$SE,RD=plan$steps$RD,RE=plan$steps$RD)
multi <- mxFitFunctionMultigroup( c("MZ","DZ") )
#ACEmodModel <- mxModel( "ACEmod", modelMZ, modelDZ, funML, multi
MainEffectsModel<- mxModel( "MainEffects", modelMZ, modelDZ, funML, multi,plan )

MainEffectsFit <- mxTryHardOrdinal(MainEffectsModel)

All fit attempts resulted in errors - check starting values or model specification

OK, maybe the custom compute plan was a bad idea.

In the main effect model, start values of moderator Path Coefficients are 0, isn’t it?

Yes, and since they're fixed, they will remain at 0 during optimization.

Do I need to change the start values of a, c, and e Path Coefficients.

You don't have to change them, but maybe there is a better choice of start values than what you're currently using. Sorry I don't have any specific advice about that.

Replied on Tue, 08/20/2019 - 14:59
Picture of user. AdminNeale Joined: 03/01/2013

Hi Julia!

A few comments:

I think that the expectation - that the estimate of A when the moderator value is zero should be the same as when no moderation is specified in the model - is incorrect. I've not read the thread in detail, but basically the estimate of A with no moderation is an average over all levels of Ac + Am * M (A constant plus A moderated * M). Depending on the distribution of M (is it symmetric around zero?), and the size of Am, you should not generally expect Ac = A.

For starting values in moderated models, I usually start moderator effect at zero. Also, if the moderator is something like age in years, it can help a lot to rescale to age in centuries to approximate a 0-1 range for the moderator value.

Finally, as you may already know, with binary data some models are not identified. See Medland, S.E., Neale, M.C., Eaves, L.J., Neale, B.M. (2009). A note on the parameterization of purcell's G x E model for ordinal and binary Data. Behavior Genetics, 39 (2): 220-229.

Replied on Tue, 08/20/2019 - 16:41
Picture of user. tbates Joined: 07/31/2009

I've tried to re-phrase the opening message to understand what's going on.

> I am running a bivariate GxE model with a dichotomous DV, and a continuous but skewed moderator not shared by twins
Is that right? Below, you mention two moderators - that's a valiant script :-)
> Estimates at the zero level of moderator with the estimates when there is no moderation in the model differ.

Non-significant variables happily take (random, drop-able without significant loss of fit) values.

> I am running ADE model even though D can be dropped out of the model but power is low.

I'm going to ask below: what's the N here?

> For one of the moderators...

Hang on... so two moderators?

> For the other two moderators
I'd say "for the other moderator"

What's your n? With a DV, and two moderators, all loading on the DV via three components, and having three moderated loadings on the DV, and moderating three influences on the DV, this will be a hard model to get stable estimates on.

Replied on Wed, 08/21/2019 - 01:53
No user picture. xiyuesenlinyu Joined: 05/24/2019

In reply to by AdminRobK

Thanks!
Now, the model runs well. Another problem, the p-value=1 between the main effect model and model with moderator.

mxCompare(ACEmodFit,MainEffectsFit)
base comparison ep minus2LL df AIC diffLL diffdf p
1 MainEffects 18 -412071.5 159463 -730997.5 NA NA NA
2 MainEffects MainEffects 15 -599113.0 159466 -918045.0 -187041.4 3 1
Replied on Wed, 08/21/2019 - 09:54
Picture of user. AdminRobK Joined: 01/24/2014

In reply to by xiyuesenlinyu

Another problem, the p-value=1 between the main effect model and model with moderator.

If I'm reading the table correctly, then I think the model with 18 free parameters didn't converge. A model with 18 free parameters should have a minus2LL no greater than a model with 15 parameters.

I don't personally find `mxCompare()` very useful. Would you mind posting the `summary()` output you get from both fitted models, as well as your output from mxVersion()?

Replied on Wed, 08/21/2019 - 21:42
No user picture. xiyuesenlinyu Joined: 05/24/2019

In reply to by AdminRobK

If I'm reading the table correctly, then I think the model with 18 free parameters didn't converge. A model with 18 free parameters should have a minus2LL no greater than a model with 15 parameters.

There were no error messages when the extraTries = 10 in the ACEmodModel(the full model with moderator). However, the following warning message appeared when I set extraTries = 30. And the warning messages disappeared when extraTries = 80. As for the output of mxCompare(ACEmodFit,MainEffectsFit) ,there were no change.

> set.seed(1)
> ACEmodFit <- mxTryHardOrdinal(ACEmodModel)
Solution found! Final fit=-476335.1 (started at 171627.54) (11 attempt(s): 7 valid, 4 errors)

> ACEmodFit <- mxTryHardOrdinal(ACEmodModel,extraTries = 30)
Solution found! Final fit=-491704.8 (started at 171627.54) (31 attempt(s): 14 valid, 17 errors)
Warning message:
In mxTryHard(model = model, greenOK = greenOK, checkHess = checkHess, :
The model does not satisfy the first-order optimality conditions to the required accuracy, and no improved point for the merit function could be found during the final linesearch (Mx status RED)

> set.seed(1)
> ACEmodFit <- mxTryHardOrdinal(ACEmodModel,extraTries = 80)
Solution found! Final fit=-476335.1 (started at 171627.54) (81 attempt(s): 35 valid, 46 errors)

mxVersion()

OpenMx version: 2.13.2 [GIT v2.13.2]
R version: R version 3.6.0 (2019-04-26)
Platform: x86_64-w64-mingw32
Default optimizer: CSOLNP
NPSOL-enabled?: Yes
OpenMP-enabled?: No


MainEffectsFit <- mxTryHardOrdinal(MainEffectsModel,extraTries = 30)

Solution found! Final fit=-599112.99 (started at 169287.04) (31 attempt(s): 14 valid, 17 errors)

> summary(MainEffectsFit)
Summary of MainEffects

free parameters:
name matrix row col Estimate Std.Error A lbound
1 aM MZ.a 1 1 1.449491e-02 7.028758e-05 ! 0!
2 aC MZ.a 2 1 8.430153e-01 9.187338e-03 -10
3 aU MZ.a 2 2 8.874368e-06 2.635970e-02 ! 0!
4 cM MZ.c 1 1 1.727313e-02 1.369073e-04 ! 0!
5 cC MZ.c 2 1 3.672863e-01 1.056641e-02 -10
6 cU MZ.c 2 2 4.954755e-06 3.044826e-02 ! 0!
7 eM MZ.e 1 1 6.485712e-04 9.930890e-06 ! 0!
8 eC MZ.e 2 1 3.361104e-01 1.644506e-02 -10
9 eU MZ.e 2 2 2.036068e-01 1.104726e-02 0
10 beta_Age_drink_cur MZ.bAge 1 1 1.000008e+00 1.929002e-05
11 beta_Age_coronary MZ.bAge 1 2 6.779762e-03 1.721201e-02
12 beta_Sex_drink_cur MZ.bSex 1 1 -3.644746e-05 2.361072e-04
13 beta_Sex_coronary MZ.bSex 1 2 -2.464578e-02 1.912105e-02
14 meanM MZ.intercept 1 1 3.261617e-03 3.507051e-04
15 coronary_thresh MZ.Threshold 1 coronary1 1.340007e+00 3.094911e-02
ubound

Model Statistics:
| Parameters | Degrees of Freedom | Fit (-2lnL units)
Model: 15 159466 -599113
Saturated: NA NA NA
Independence: NA NA NA
Number of observations/statistics: 39870/159481

Constraint 'DZ.Var1' contributes 1 observed statistic.

Information Criteria:
| df Penalty | Parameters Penalty | Sample-Size Adjusted
AIC: -918045 -599083.0 -599083.0
BIC: -2288397 -598954.1 -599001.8
CFI: NA
TLI: 1 (also known as NNFI)
RMSEA: 0 [95% CI (NA, NA)]
Prob(RMSEA <= 0.05): NA
To get additional fit indices, see help(mxRefModels)
timestamp: 2019-08-22 08:17:54
Wall clock time: 4.304247 secs
OpenMx version number: 2.13.2
Need help? See help(mxSummary)

> set.seed(1)
> ACEmodFit <- mxTryHardOrdinal(ACEmodModel)

Solution found! Final fit=-476335.1 (started at 171627.54) (11 attempt(s): 7 valid, 4 errors)

> summary(ACEmodFit)
Summary of MainEffects

free parameters:
name matrix row col Estimate Std.Error A lbound
1 aM MZ.a 1 1 7.190948e-04 -7.407524e+00 ! 0!
2 aC MZ.a 2 1 -2.504145e-02 1.653835e-01 -10
3 aU MZ.a 2 2 -5.604461e-08 1.628159e-07 0!
4 cM MZ.c 1 1 1.919975e-04 -2.375641e+01 ! 0!
5 cC MZ.c 2 1 -1.971316e-01 6.797870e-01 -10
6 cU MZ.c 2 2 1.692647e-02 -5.648652e-02 ! 0!
7 eM MZ.e 1 1 1.960369e-02 -1.176272e+01 ! 0!
8 eC MZ.e 2 1 9.799110e-01 -8.964625e-01 -10
9 eU MZ.e 2 2 -1.775966e-08 3.622713e-08 0!
10 aMod21 MZ.aMod 2 1 -2.648286e+01 -2.601806e-06
11 aMod22 MZ.aMod 2 2 -4.271794e-04 -7.992276e-10
12 cMod21 MZ.cMod 2 1 -7.242559e+01 -4.082045e-06
13 cMod22 MZ.cMod 2 2 2.722510e-01 6.434716e-08 !
14 eMod21 MZ.eMod 2 1 4.277475e+02 4.941039e-06
15 eMod22 MZ.eMod 2 2 -1.684815e-04 -4.634743e-10
16 beta_Age_drink_cur MZ.bAge 1 1 1.000149e+00 -1.141624e-01
17 beta_Age_coronary MZ.bAge 1 2 -4.276868e+02 5.793846e-06
18 beta_Sex_drink_cur MZ.bSex 1 1 -1.002649e-04 8.014452e+01
19 beta_Sex_coronary MZ.bSex 1 2 -2.274997e-02 -1.474523e+00
20 meanM MZ.intercept 1 1 4.591315e-03 5.418001e+01
21 coronary_thresh MZ.Threshold 1 coronary1 9.840132e-01 1.179441e+00
ubound

Model Statistics:
| Parameters | Degrees of Freedom | Fit (-2lnL units)
Model: 21 159460 -476335.1
Saturated: NA NA NA
Independence: NA NA NA
Number of observations/statistics: 39870/159481

Constraint 'DZ.Var1' contributes 1 observed statistic.

Information Criteria:
| df Penalty | Parameters Penalty | Sample-Size Adjusted
AIC: -795255.1 -476293.1 -476293.1
BIC: -2165555.4 -476112.6 -476179.4
CFI: NA
TLI: 1 (also known as NNFI)
RMSEA: 0 [95% CI (NA, NA)]
Prob(RMSEA <= 0.05): NA
To get additional fit indices, see help(mxRefModels)
timestamp: 2019-08-22 08:59:39
Wall clock time: 6.842392 secs
OpenMx version number: 2.13.2
Need help? See help(mxSummary)

> mxCompare(ACEmodFit,MainEffectsFit)
base comparison ep minus2LL df AIC diffLL diffdf p
1 MainEffects 21 -476335.1 159460 -795255.1 NA NA NA
2 MainEffects MainEffects 15 -599113.0 159466 -918045.0 -122777.9 6 1

>

Replied on Sat, 08/24/2019 - 11:11
Picture of user. AdminNeale Joined: 03/01/2013

In reply to by xiyuesenlinyu

When a sub model with fewer parameters fits better than its supermodel, the fit of the supermodel has to be wrong, because it should fit at least as well as the sub model. In such situations, poor starting values for the supermodel are often the culprit. I suggest that you take the fitted sub model and free up the parameters (or stop having two parameters equated, estimate them separately) necessary to turn the sub model into the supermodel. Then start fitting the supermodel from the sub model’s solution. At the least, the fit of the supermodel should not get any worse than that of the sub model (which is where the supermodel optimization is beginning).
Replied on Sun, 08/25/2019 - 03:47
No user picture. xiyuesenlinyu Joined: 05/24/2019

In reply to by AdminNeale

Thanks for your suggestions, but I can’t understand well!
“I suggest that you take the fitted sub model and free up the parameters (or stop having two parameters equated, estimate them separately) necessary to turn the sub model into the supermodel.”
Did you mean to reset the three parameters like those
modPathA = mxMatrix( "Lower", nrow=nv, ncol=nv, free=c(F,T,T), values=pathModVal, labels=aModLabs, name="aMod" )
modPathC = mxMatrix( "Lower", nrow=nv, ncol=nv, free=c(F,T,T), values=pathModVal, labels=cModLabs, name="cMod" )
modPathE = mxMatrix( "Lower", nrow=nv, ncol=nv, free=c(F,T,T), values=pathModVal, labels=eModLabs, name="eMod" )
As for “or stop having two parameters equated, estimate them separately”. Though I can estimate them separately(other submodels), I need to estimate all of them in full model.
Another question, the -2LL are negative in some models, are they ok?
Replied on Tue, 08/27/2019 - 11:00
No user picture. xiyuesenlinyu Joined: 05/24/2019

Hello, the model can be run now. however, I got puzzling output: the differences of -2LL are significant between the two models, and an odd p-value(0) for mxcompare().
Could you give some suggestions to me?
Thank you!

> mxCompare(ACEmodFit,MainEffectsFit)
base comparison ep minus2LL df AIC diffLL diffdf p
1 ACEmoderation 21 40711.46 159460 -278208.5 NA NA NA
2 ACEmoderation MainEffects 15 63773.39 159466 -255158.6 23061.93 6 0
Replied on Wed, 08/28/2019 - 10:09
Picture of user. AdminRobK Joined: 01/24/2014

In reply to by xiyuesenlinyu

You're comparing the moderation model to the main-effects model. A nonempty proper subset of the free parameters in the moderation model is fixed in the main-effects model, thus, the main-effects model is said to be a "nested submodel" of the moderation model. The p-value of 0 is telling you that the main-effects model provides much, much worse fit to the data than the moderation model. Or, to put it another way, you can reject the null hypothesis that the six moderation parameters are all equal to zero.
Replied on Wed, 08/28/2019 - 13:38
Picture of user. AdminRobK Joined: 01/24/2014

In reply to by tbates

I don't agree with your use of the word "incorrectly". The p = 0 being discussed in this thread is being reported computationally correctly. Mathematically, the p-value is nonzero, but that value is too small to be represented as a double-precision floating-point value. For that matter--Regarding issue #131, p was also being reported computationally correctly. But there, we decided to sacrifice computational correctness for user experience in an edge case.