Moderation model produces weird results

Posted on Tue, 04/16/2019 - 14:26

Julia Joined: 03/29/2012

Forums

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

AdminRobK Joined: 01/24/2014

script?

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

Julia Joined: 03/29/2012

Thank you for the answer!

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!

File attachments

Moderation output.txt

bivChol_Moderation.txt

bivCF_Moderation.txt

Replied on Wed, 04/17/2019 - 11:41

AdminRobK Joined: 01/24/2014

some input

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

Julia Joined: 03/29/2012

I tried

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

AdminRobK Joined: 01/24/2014

mxTryHardOrdinal() ?

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

Julia Joined: 03/29/2012

answers

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

AdminRobK Joined: 01/24/2014

results are probably right

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

Julia Joined: 03/29/2012

Yes, what is puzzling me is

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 Thu, 04/18/2019 - 15:55

AdminRobK Joined: 01/24/2014

It's been a few years since I

It's been a few years since I read that Rathouz et al. paper, so there's a good chance I'm mistaken in what I posted about the correlated-factors parameterization.

Replied on Sun, 08/18/2019 - 10:14

xiyuesenlinyu Joined: 05/24/2019

Hello

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

xiyuesenlinyu Joined: 05/24/2019

Hello

Replied on Sun, 08/18/2019 - 16:25

AdminRobK Joined: 01/24/2014

start values

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

xiyuesenlinyu Joined: 05/24/2019

I still can’t run my model.

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

AdminRobK Joined: 01/24/2014

some suggestions

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

xiyuesenlinyu Joined: 05/24/2019

I change the order of models,

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

AdminRobK Joined: 01/24/2014

First of all,

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

AdminNeale Joined: 03/01/2013

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

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

tbates Joined: 07/31/2009

isolating the question

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 Tue, 08/20/2019 - 18:01

AdminRobK Joined: 01/24/2014

user xiyuesenlinyu?

AdminNeale, tbates: Any advice for xiyuesenlinyu (not the OP)? They're the reason I emailed the list about this thread.

Replied on Wed, 08/21/2019 - 01:53

xiyuesenlinyu Joined: 05/24/2019

Thanks!

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

AdminRobK Joined: 01/24/2014

that doesn't seem right

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

xiyuesenlinyu Joined: 05/24/2019

If I'm reading the table

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 1                  aM 2                  aC 3                  aU 4                  cM 5                  cC 6                  cU 7                  eM 8                  eC 9                  eU 10             aMod21 11             aMod22 12             cMod21 13             cMod22 14             eMod21 15             eMod22 16 beta_Age_drink_cur 17  beta_Age_coronary 18 beta_Sex_drink_cur 19  beta_Sex_coronary 20 21    coronary_thresh ubound

matrix row       col      Estimate     Std.Error A lbound
 MZ.a   1         1  7.190948e-04 -7.407524e+00 !     0!
 MZ.a   2         1 -2.504145e-02  1.653835e-01      -10
 MZ.a   2         2 -5.604461e-08  1.628159e-07       0!
 MZ.c   1         1  1.919975e-04 -2.375641e+01 !     0!
 MZ.c   2         1 -1.971316e-01  6.797870e-01      -10
 MZ.c   2         2  1.692647e-02 -5.648652e-02 !     0!
 MZ.e   1         1  1.960369e-02 -1.176272e+01 !     0!
 MZ.e   2         1  9.799110e-01 -8.964625e-01      -10
 MZ.e   2         2 -1.775966e-08  3.622713e-08       0!
 MZ.aMod   2         1 -2.648286e+01 -2.601806e-06
 MZ.aMod   2         2 -4.271794e-04 -7.992276e-10
 MZ.cMod   2         1 -7.242559e+01 -4.082045e-06
 MZ.cMod   2         2  2.722510e-01  6.434716e-08 !
 MZ.eMod   2         1  4.277475e+02  4.941039e-06
 MZ.eMod   2         2 -1.684815e-04 -4.634743e-10
 MZ.bAge   1         1  1.000149e+00 -1.141624e-01
 MZ.bAge   1         2 -4.276868e+02  5.793846e-06
 MZ.bSex   1         1 -1.002649e-04  8.014452e+01
 MZ.bSex   1         2 -2.274997e-02 -1.474523e+00
 meanM MZ.intercept   1         1  4.591315e-03  5.418001e+01
 MZ.Threshold   1 coronary1  9.840132e-01  1.179441e+00

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 Fri, 08/23/2019 - 05:01

xiyuesenlinyu Joined: 05/24/2019

Mr Robert

Mr Robert
Could you give some suggestions for the output of P-value=1?
Please!

Replied on Fri, 08/23/2019 - 11:04

AdminRobK Joined: 01/24/2014

script?

Please post your current script, preferably as an attachment.

Replied on Sat, 08/24/2019 - 10:37

xiyuesenlinyu Joined: 05/24/2019

ok, thank you!

File attachments

25+ACE-pucell-drink_cur.R

Replied on Mon, 08/26/2019 - 10:13

AdminRobK Joined: 01/24/2014

Create moderation model from fitted main-effects model

Try replacing this (line 249),



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.

Replied on Sat, 08/24/2019 - 11:11

AdminNeale Joined: 03/01/2013

A trick to help sticky optimization

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

xiyuesenlinyu Joined: 05/24/2019

Thanks for your suggestions,

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

xiyuesenlinyu Joined: 05/24/2019

Hello, the model can be run

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 Tue, 08/27/2019 - 12:11

AdminRobK Joined: 01/24/2014

A chi-square test statistic

A chi-square test statistic of 23061.93 on 6df has a p-value computationally equal to zero:



> pchisq(23061.93,df=6,lower.tail=F)

[1] 0

Replied on Tue, 08/27/2019 - 19:10

xiyuesenlinyu Joined: 05/24/2019

I’m sorry, I’m new for OpenMx

I’m sorry, I’m new for OpenMx, even for R. I have no idea for the output. Could you tell me what does that mean with a p-value=0? An Inappropriate model or poor model fitting? Is the model reliable? What can I do for improvement?

Replied on Wed, 08/28/2019 - 10:09

AdminRobK Joined: 01/24/2014

You're comparing the

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 - 10:27

xiyuesenlinyu Joined: 05/24/2019

Thanks for your help. I think

Thanks for your help. I think that P-value <0.05 between two models may be significant, which suggests the moderated effects. My concern is p=0(just right), which seems not quite right.

Replied on Wed, 08/28/2019 - 11:45

AdminRobK Joined: 01/24/2014

underflow

The output of `pchisq()` is underflowing to zero. For the sake of perspective, the 99.99th percentile of a chi-square distribution on 6df is about 27.9:



> qchisq(0.9999,df=6)

[1] 27.85634

Replied on Wed, 08/28/2019 - 12:08

tbates Joined: 07/31/2009

Disconcerting to have p-values, fit value reported incorrectly

This issue is the same: p = 1 being reported as p =0) : https://github.com/OpenMx/OpenMx/issues/131

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

Replied on Wed, 08/28/2019 - 13:38

AdminRobK Joined: 01/24/2014

not "incorrectly"

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.

News

Moderation model produces weird results

News