# Moderation model produces weird results

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

## script?

Log in or register to post comments

## Thank you for the answer!

Thank you for looking at it!

Log in or register to post comments

## some input

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?

Log in or register to post comments

In reply to some input by AdminRobK

## I tried

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.

Log in or register to post comments

In reply to I tried by Julia

## mxTryHardOrdinal() ?

Did the fit value appreciably improve?

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.

But why isn't it a free parameter?

Have you tried `mxTryHardOrdinal()`?

Log in or register to post comments

In reply to mxTryHardOrdinal() ? by AdminRobK

## answers

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?

Log in or register to post comments

In reply to answers by Julia

## results are probably right

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.

Log in or register to post comments

In reply to results are probably right by AdminRobK

## Yes, what is puzzling me is

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).

Log in or register to post comments

In reply to Yes, what is puzzling me is by Julia

## It's been a few years since I

Log in or register to post comments

## 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!

Log in or register to post comments

## 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!

Log in or register to post comments

In reply to Hello by xiyuesenlinyu

## start values

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.

Log in or register to post comments

In reply to start values 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!

Log in or register to post comments

In reply to I still can’t run my model. by xiyuesenlinyu

## some suggestions

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`.

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.What were you trying to do? Use

`omxSetParameters()`

to change free parameter values.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`

.Log in or register to post comments

In reply to some suggestions by AdminRobK

## I change the order of models,

> 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" )

Log in or register to post comments

In reply to I change the order of models, by xiyuesenlinyu

## First of all,

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

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

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

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.

Log in or register to post comments

## Premise: A without moderation = A when moderated by zero?

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.

Log in or register to post comments

## isolating the question

> 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.

Log in or register to post comments

## user xiyuesenlinyu?

xiyuesenlinyu(not the OP)? They're the reason I emailed the list about this thread.Log in or register to post comments

In reply to user xiyuesenlinyu? 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

2 MainEffects MainEffects 15 -599113.0 159466 -918045.0 -187041.4 3 1

Log in or register to post comments

In reply to Thanks! by xiyuesenlinyu

## that doesn't seem right

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()`

?Log in or register to post comments

In reply to that doesn't seem right by AdminRobK

## If I'm reading the table

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) 21 -476335.1 159460 -795255.1 NA NA NA

base comparison ep minus2LL df AIC diffLL diffdf p

1 MainEffects

2 MainEffects MainEffects 15 -599113.0 159466 -918045.0 -122777.9 6 1

`>`

Log in or register to post comments

In reply to If I'm reading the table by xiyuesenlinyu

## Mr Robert

Could you give some suggestions for the output of P-value=1?

Please!

Log in or register to post comments

In reply to Mr Robert by xiyuesenlinyu

## script?

Log in or register to post comments

In reply to script? by AdminRobK

## ok, thank you!

Log in or register to post comments

In reply to ok, thank you! by xiyuesenlinyu

## Create moderation model from fitted main-effects model

ACEmodModel <- mxModel (MainEffectsModel, name='MainEffects')

, with this,

ACEmodModel <- mxModel (MainEffectsFit, name='ACEmoderation')

. That way, you'll be starting the moderation model from the solution for the main-effects model.

Log in or register to post comments

In reply to Mr Robert by xiyuesenlinyu

## A trick to help sticky optimization

Log in or register to post comments

In reply to A trick to help sticky optimization by AdminNeale

## Thanks for your suggestions,

“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?

Log in or register to post comments

## Hello, the model can be run

Could you give some suggestions to me?

Thank you!

> mxCompare(ACEmodFit,MainEffectsFit)

base comparison ep minus2LL df AIC diffLL diffdf p

1 ACEmoderation

2 ACEmoderation MainEffects 15 63773.39 159466 -255158.6 23061.93 6 0

Log in or register to post comments

In reply to Hello, the model can be run by xiyuesenlinyu

## A chi-square test statistic

dfhas ap-value computationally equal to zero:> pchisq(23061.93,df=6,lower.tail=F)

[1] 0

Log in or register to post comments

In reply to A chi-square test statistic by AdminRobK

## I’m sorry, I’m new for OpenMx

Log in or register to post comments

In reply to I’m sorry, I’m new for OpenMx by xiyuesenlinyu

## You're comparing 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.Log in or register to post comments

In reply to You're comparing the by AdminRobK

## Thanks for your help. I think

Log in or register to post comments

In reply to Thanks for your help. I think by xiyuesenlinyu

## underflow

dfis about 27.9:> qchisq(0.9999,df=6)

[1] 27.85634

Log in or register to post comments

## Disconcerting to have p-values, fit value reported incorrectly

Also this issue where RMSEA and TLI are misreported) https://github.com/OpenMx/OpenMx/issues/221

Log in or register to post comments

In reply to Disconcerting to have p-values, fit value reported incorrectly by tbates

## not "incorrectly"

p= 0 being discussed in this thread is being reported computationally correctly. Mathematically, thep-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,pwas also being reported computationally correctly. But there, we decided to sacrifice computational correctness for user experience in an edge case.Log in or register to post comments