how to avoid NA in confidence interval

Posted on Thu, 10/10/2019 - 10:06

diana Joined: 04/17/2019

Attachment

Size

3.png

33.17 KB

Forums

General SEM Discussions

Hello, I often got NA in confidence interval when fitting univariate or bivariate SEM and changing optimizers didn't help. Here is the code of a bivariate ACE model, and the output when I use 'SLSQP' optimizer . I wonder how can I aviod NA when estimating CI of parameters and algebra. Thank you so much!



vars	  <-c('tea_5g','smk_num_4g')

covars  <- c("age","region_type")

nv 	    <- 2		                                         # number of variables per twin

ncv     <- 2                         # number of covariates

ntv     <- nv*2	                                         # number of variables per pair

nth	    <- 3		                                         # number of max thresholds = the number of category - 1

selVars	<-paste(vars,c(rep(1,nv),rep(2,nv)),sep="")  #c('tea2g1', 'smk_num_4g1', 'tea2g2', 'smk_num_4g2' )#

covVars <- paste(covars,c(rep(1,ncv),rep(2,ncv)),sep="")
mzData <- subset(tea4g_smk4g_nodrink,zygosity1==0,c(selVars,covVars))

dzData <- subset(tea4g_smk4g_nodrink,zygosity1==1,c(selVars,covVars))

sum(is.na(mzData))

sum(is.na(dzData))

mzData <- na.omit(mzData)

dzData <- na.omit(dzData)

dim(mzData);dim(dzData)

# Generate Descriptive Statistics

mzData[,selVars] 	<-mxFactor(mzData[,selVars], levels=c(0:nth) )

dzData[,selVars] 	<-mxFactor(dzData[,selVars], levels=c(0:nth) )
# Set Starting Values

svb  <- 0.1                                    # sv for beta,get from tryhardordinal

svLTh <- -1.5                                  # start value for first threshold

svITh <- 1                                     # start value for increments

svTh <- matrix(rep(c(svLTh,(rep(svITh,nth-1)))),nrow=nth,ncol=nv) # start value for thresholds

lbTh <- matrix(rep(c(-3,(rep(0.001,nth-1))),nv),nrow=nth,ncol=nv) # lower bounds for thresholds

svCor <- .5                                        # start value for correlations

# Create Labels

labTh <- function(lab,vars,nth) { paste(paste("t",1:nth,lab,sep=""),rep(vars,each=nth),sep="") }

labdata <- cbind(paste("data.",covars[1],c(1,2),sep = ""),paste("data.",covars[2],c(1,2),sep = "")) #c("data.age1","data.age2","data.region_type1","data.region_type2") labeta <-cbind(paste("B",covars,"_",vars[1],sep = ""),paste("B",covars,"_",vars[2],sep = "")) lathmz <- c(labTh("MZ","tea_5g",nth),labTh("Z","smk_num_4g",nth),labTh("MZ","tea_5g",nth),labTh("Z","smk_num_4g",nth)) lathdz <- c(labTh("DZ","tea_5g",nth),labTh("Z","smk_num_4g",nth),labTh("DZ","tea_5g",nth),labTh("Z","smk_num_4g",nth)) ### Some parameters included in all "submodels": baseACE <- mxModel('Base', mxMatrix( type="Full", nrow=ncv, ncol=nv, free=T, values=svb,labels=labeta, name="Beta"), mxMatrix( type="Zero", nrow=1, ncol=nv, name="meanG" ), # or type="Zero" # Create Matrices for Variance Components mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=c(2,0.1,2), name="A" ), # label=laLower("A",nv), values=valDiag(svPa,nv),# mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=c(2,0.1,2), name="C" ), mxMatrix( type="Symm", nrow=nv, ncol=nv, free=c(F,T,F), values=c(1, 0.1, 1), name="E" ),# fix E to avoid warning,without warning, E can be freely estimated mxAlgebra( expression=A+C+E, name="V" ), mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I"), mxAlgebra( expression=solve(sqrt(I*V)), name="iSD"), mxAlgebra( expression=solve(sqrt(I*V))%&%V, name ="Rph"), mxAlgebra( expression=solve(sqrt(I*A))%&%A, name ="Ra" ), #cov2cor() mxAlgebra( expression=solve(sqrt(I*C))%&%C, name ="Rc" ), mxAlgebra( expression=solve(sqrt(I*E))%&%E, name ="Re" ), # Algebra to compute standardized variance components mxAlgebra( expression=A/V, name="h2"), mxAlgebra( expression=C/V, name="c2"), mxAlgebra( expression=E/V, name="e2"), # the proportion of ace contributing to the overlap mxAlgebra( expression=sqrt(h2[1,1])*sqrt(h2[2,2])*Ra[1,2]/Rph[1,2],name = "pa"), mxAlgebra( expression=sqrt(c2[1,1])*sqrt(c2[2,2])*Rc[1,2]/Rph[1,2],name = "pc"), mxAlgebra( expression=sqrt(e2[1,1])*sqrt(e2[2,2])*Re[1,2]/Rph[1,2],name = "pe") ) modelMZ <- mxModel( "MZ", mxMatrix( type="Full", nrow=2, ncol=ncv, free=F, label=labdata, name="MZDefVars"), mxAlgebra( expression= cbind(Base.meanG + MZDefVars[1,] %*% Base.Beta,Base.meanG + MZDefVars[2,] %*% Base.Beta), name="expMeanMZ" ), mxMatrix( type="Full", nrow=nth, ncol=ntv, free=T, values=svTh, labels=lathmz, name="ThMZ"), mxMatrix( type="Lower", nrow=nth, ncol=nth, free=FALSE, values=1, name="inc" ), mxAlgebra( expression= inc %*% ThMZ, name="threMZ" ), mxAlgebra( expression= rbind( cbind(Base.A+Base.C+Base.E , Base.A+Base.C), cbind(Base.A+Base.C , Base.A+Base.C+Base.E)), name="expCovMZ" ), mxData( observed=mzData, type="raw" ), mxExpectationNormal( covariance="expCovMZ", means="expMeanMZ", dimnames=selVars, thresholds="threMZ" ), mxFitFunctionML() ) modelDZ <- mxModel( "DZ", mxMatrix( type="Full", nrow=2, ncol=ncv, free=F, label=labdata, name="DZDefVars"), mxAlgebra( expression= cbind(Base.meanG + DZDefVars[1,] %*% Base.Beta,Base.meanG + DZDefVars[2,] %*% Base.Beta), name="expMeanDZ" ), mxMatrix( type="Full", nrow=nth, ncol=ntv, free=T, values=svTh,labels=lathdz, name="ThDZ"), mxMatrix( type="Lower", nrow=nth, ncol=nth, free=FALSE, values=1, name="inc" ), mxAlgebra( expression= inc %*% ThDZ, name="threDZ" ), mxAlgebra( expression= rbind( cbind(Base.A+Base.C+Base.E , 0.5%x%Base.A+Base.C), cbind(0.5%x%Base.A+Base.C , Base.A+Base.C+Base.E)), name="expCovDZ" ), mxData( observed=dzData, type="raw" ), mxExpectationNormal( covariance="expCovDZ", means="expMeanDZ", dimnames=selVars, thresholds="threDZ" ), mxFitFunctionML() ) Conf1 <- mxCI (c ('Base.h2[1,1]','Base.c2[1,1]','Base.e2[1,1]','Base.h2[2,2]','Base.c2[2,2]','Base.e2[2,2]')) Conf2 <- mxCI (c ('Base.Ra[2,1]', 'Base.Rc[2,1]', 'Base.Re[2,1]','Base.Rph[2,1]') ) Conf3 <- mxCI (c ("Base.pa","Base.pc","Base.pe") ) AceModel <- mxModel( "ACE", baseACE,modelMZ, modelDZ, Conf1, Conf2,Conf3, mxFitFunctionMultigroup(c('MZ.fitfunction','DZ.fitfunction'))) # ------------------------------------------------------------------------------ # 4) RUN AceModel AceFit <- mxTryHardOrdinal(AceModel, intervals=T,extraTries = 15,OKstatuscodes=c(0,1,5)) (AceSumm <- summary(AceFit,verbose = T))

Replied on Thu, 10/10/2019 - 12:22

AdminRobK Joined: 01/24/2014

suggestions

First off, is your installation of OpenMx NPSOL-enabled? You can tell from the output of `mxVersion()`. If it is, and you haven't tried switching optimizers to NSPOL, then try that first.

Another possibility is to use mxBootstrap(), and then use mxBootstrapEval() and mxBootstrapEvalByName() to get confidence intervals for elements of MxAlgebras.

I expect bootstrapping will work OK if you use enough replications, but if not, you could make a custom compute plan to use a derivative-free optimizer for the confidence-limit search:
AceFitCI <- AceFit AceFitCI$compute$steps <- list(CI=omxDefaultComputePlan(intervals=T)$steps$CI) AceFitCI$compute$steps$CI$constraintType <- "ineq" AceFitCI$compute$steps$CI$plan <- mxComputeNelderMead(ineqConstraintMthd="eqMthd",centerIniSimplex=T) # Alternative plan: # AceFitCI$compute$steps$CI$constraintType <- "none" # AceFitCI$compute$steps$CI$plan <- mxComputeNelderMead() AceFitCI <- mxRun(AceFitCI) summary(AceFitCI,verbose=T)

If all else fails, you could use mxSE() to get standard errors for your MxAlgebra elements for which you still don't have valid CIs, and form Wald confidence intervals for them. But, only do that as a last resort.

Log in or register to post comments

Replied on Tue, 10/15/2019 - 10:19

diana Joined: 04/17/2019

error of running mxBootstrap()

I have already tried NSPOL and CI still contained NA.

I am sorry it is my first time to use mxBootstrap(), and I don't know why this error was turned out after running this code:
testBoot <- mxBootstrap(AceFit,10)
and the error was:
*Error in strsplit(e1, imxSeparatorChar, fixed = TRUE) :
non-character argument*
Can you tell what was wrong?

Thanks!

Log in or register to post comments

Replied on Tue, 10/15/2019 - 11:03

AdminRobK Joined: 01/24/2014

traceback()

You may be encountering a bug. Try running `traceback()` right after you get the error message and report what it says.

Log in or register to post comments

Replied on Tue, 10/15/2019 - 22:56

diana Joined: 04/17/2019

hello, this was the result

> testBoot  <-  mxBootstrap(AceFit,10)

Error in strsplit(e1, imxSeparatorChar, fixed = TRUE) :

  non-character argument

> traceback()

6: strsplit(e1, imxSeparatorChar, fixed = TRUE)

5: unlist(strsplit(e1, imxSeparatorChar, fixed = TRUE))

4: FUN(X[[i]], ...)

3: vapply(data, function(e1) {

       path <- unlist(strsplit(e1, imxSeparatorChar, fixed = TRUE))

       if (length(path) == 1) {

           e1 <- paste(path, "data", sep = imxSeparatorChar)

       }

       e1

   }, "")

2: mxComputeBootstrap(data, model@compute)

1: mxBootstrap(AceFit, 10)

Log in or register to post comments

Replied on Wed, 10/16/2019 - 07:07

jpritikin Joined: 05/23/2012

not helpful

Unfortunately this is not helpful. I already figured out the call stack. The question is how to reproduce the problem. See my post below.

Log in or register to post comments

Replied on Wed, 10/16/2019 - 10:17

AdminRobK Joined: 01/24/2014

mxVersion?

diana, what is your `mxVersion()` output?

Log in or register to post comments

Replied on Wed, 10/16/2019 - 10:48

diana Joined: 04/17/2019

> mxVersion()

> mxVersion()
OpenMx version: 2.12.2 [GIT v2.12.2]
R version: R version 3.5.0 (2018-04-23)
Platform: x86_64-w64-mingw32
Default optimizer: NPSOL
NPSOL-enabled?: Yes
OpenMP-enabled?: No

Log in or register to post comments

Replied on Wed, 10/16/2019 - 11:06

AdminRobK Joined: 01/24/2014

update and try again

The latest release of OpenMx is v2.14.11. Try updating, and see if the issue persists afterward.

Log in or register to post comments

Replied on Thu, 10/17/2019 - 01:56

diana Joined: 04/17/2019

I updated OpenMx to the

I updated OpenMx to the latest version but it didn't work.
> mxVersion()
OpenMx version: 2.14.11 [GIT v2.14.11-dirty]
R version: R version 3.5.1 (2018-07-02)
Platform: x86_64-w64-mingw32
Default optimizer: NPSOL
NPSOL-enabled?: Yes
OpenMP-enabled?: No
-----------------------------------------------------------------------------
> testBoot <- mxBootstrap(AceFit,10)
Error in strsplit(e1, imxSeparatorChar, fixed = TRUE) :
non-character argument
-----------------------------------------------------------------------------
> traceback()
6: strsplit(e1, imxSeparatorChar, fixed = TRUE)
5: unlist(strsplit(e1, imxSeparatorChar, fixed = TRUE))
4: FUN(X[[i]], ...)
3: vapply(data, function(e1) {
path <- unlist(strsplit(e1, imxSeparatorChar, fixed = TRUE))
if (length(path) == 1) {
e1 <- paste(path, "data", sep = imxSeparatorChar)
}
e1
}, "")
2: mxComputeBootstrap(data, model@compute)
1: mxBootstrap(AceFit, 10)

Log in or register to post comments

Replied on Thu, 10/17/2019 - 09:14

jpritikin Joined: 05/23/2012

how to reproduce?

Please provide a [minimal working example](https://en.wikipedia.org/wiki/Minimal_working_example)

Log in or register to post comments

Replied on Thu, 10/17/2019 - 21:54

diana Joined: 04/17/2019

Thank you very much. I am

Thank you very much. I am trying to sample the data and make sure the data can reproduce the problem.

Log in or register to post comments

Replied on Fri, 10/18/2019 - 04:58

diana Joined: 04/17/2019

Hello! I am so sorry for

Hello! I am so sorry for asking questions without a minimal working example. My data generated by mxGenerateData() was uploaded and the code was paseted above.

Thank you very much for your help!

Log in or register to post comments

Replied on Thu, 10/17/2019 - 11:32

AdminRobK Joined: 01/24/2014

mxGenerateData()

I second Joshua's request for a minimal working example. mxGenerateData() may be useful for the purpose.

Log in or register to post comments

Replied on Fri, 10/18/2019 - 04:42

diana Joined: 04/17/2019

Thank you

Thank you so much! The data I uploaded was generated by the code:

samdata <- mxGenerateData(AceModel, returnModel=FALSE, use.miss = TRUE)

mzData <- samdata$MZ

dzData <- samdata$DZ

The code entailed to reproduce the problem is above pasted and you can skip the step of generating the mzData and dzData.

So, all in all, the problem I have is the NA in the CI for elements of MxAlgebras, and the error of running testBoot <- mxBootstrap(AceFit,10)

Thank you very much for your help !

File attachments

mzData.csv

dzData.csv

Log in or register to post comments

Replied on Fri, 10/18/2019 - 08:29

jpritikin Joined: 05/23/2012

fixed

Ah, finally. This bug is fixed in v2.14.11-24-gfd7b726f5. I think you can get your code working without an upgrade if you do,

AceBoot <- mxModel(AceFit, mxComputeBootstrap(c("MZ","DZ"), AceFit@compute)) testBoot <- mxBootstrap(AceBoot,10)

Thanks for your help.

Log in or register to post comments

Replied on Sun, 10/20/2019 - 04:48

diana Joined: 04/17/2019

There still exists an error

Thank you for your debugging, but there still exists an error.

> AceBoot <- mxModel(AceFit, mxComputeBootstrap(c("MZ","DZ"), AceFit@compute))

testBoot  <-  mxBootstrap(AceBoot,10)

results:
Running ACE with 20 parameters
MxComputeGradientDescent(NPSOL) evaluations 52732 fit 56009.5 change -8.307e-005Error in runHelper(model, frontendStart, intervals, silent, suppressWarnings, :
Can only compute CIs once

Log in or register to post comments

Replied on Sun, 10/20/2019 - 07:41

jpritikin Joined: 05/23/2012

mxCI

Yes, if you're going to use this work-around then you'll need to omit mxCI from the bootstrapped model.

Log in or register to post comments

Replied on Tue, 10/22/2019 - 07:59

diana Joined: 04/17/2019

Here are something confusing

1) When I use mxCI to compute CI, **Base.Ra is 0.79(NA,NA)**

2) When I use
AceFitCI <- AceFit AceFitCI$compute$steps <- list(CI=omxDefaultComputePlan(intervals=T)$steps$CI) AceFitCI$compute$steps$CI$constraintType <- "ineq" AceFitCI$compute$steps$CI$plan <- mxComputeNelderMead(ineqConstraintMthd="eqMthd",centerIniSimplex=T) AceFitCI <- mxRun(AceFitCI) summary(AceFitCI,verbose=T)
**Base.Ra is 0.79(0.18,NA)**

When I use
AceFitCI <- AceFit AceFitCI$compute$steps <- list(CI=omxDefaultComputePlan(intervals=T)$steps$CI) AceFitCI$compute$steps$CI$constraintType <- "none" AceFitCI$compute$steps$CI$plan <- mxComputeNelderMead() AceFitCI <- mxRun(AceFitCI) summary(AceFitCI,verbose=T)
**Base.Ra is 0.79(0.18,NA)**

3) When I use
AceBoot <- mxModel(AceFit, mxComputeBootstrap(c("MZ","DZ"), AceFit@compute)) testBoot <- mxBootstrap(AceBoot,10) mxBootstrapEvalByName("Base.Ra", testBoot, bq=c(.025,.975))
**result is :**
SE 2.5% 97.5%
[1,] 1.790181e-16 1.0000000 1.00000
[2,] 2.512318e-01 0.5509924 1.20444
[3,] 2.512318e-01 0.5509924 1.20444
[4,] 1.489520e-16 1.0000000 1.00000
**How can I get the value and CI of Ra?**

4)When I use
AceBoot <- mxModel(AceFit, mxComputeBootstrap(c("MZ","DZ"), AceFit@compute)) testBoot <- mxBootstrap(AceBoot,50) mxBootstrapEvalByName("Base.Ra", testBoot, bq=c(.025,.975))

an error and two warnings exist:
Error: The following error occurred while evaluating the subexpression 'solve(sqrt(Base.I * Base.A))' during the evaluation of 'Base.Ra' in model 'ACE' : Lapack routine dgesv: system is exactly singular: U[1,1] = 0
In addition: Warning messages:
1: Only 38% of the bootstrap replications converged acceptably. Accuracy is much less than the 50 replications requested. Examine table(model$compute$output$raw$statusCode)
2: In sqrt(c(-0.137782917394753, 0, 0, 0.97874966037632)) : NaNs produced

**I wonder why the error and warnings exist after I change the number of resampling replications, is there anything wrong?
Is it ok to use the result calculated by replying resample 10 times?**

Log in or register to post comments

Replied on Tue, 10/22/2019 - 12:11

AdminRobK Joined: 01/24/2014

some answers

How can I get the value and CI of Ra?

In your script, 'Base.Ra' is an MxAlgebra that evaluates to a 2x2 matrix. `mxBootstrapEvalByName()` is outputting the standard errors and confidence limits for the 4 elements of that matrix. To see the elements themselves, just do `mxEvalByName("Base.Ra", AceFit, T)` or `mxEval(Ra, AceFit$Base, T)`.

I wonder why the error and warnings exist after I change the number of resampling replications, is there anything wrong?

More replications means more calls to `mxRun()`, and therefore more opportunities for errors or warnings to be raised. That's all there is to it.

Is it ok to use the result calculated by replying resample 10 times?

Absolutely not. If you want 95% CIs, you need to be using at minimum 1000 replications.

Base.Ra is 0.79(0.18,NA)

What does the corresponding row of the CI details table say about that upper confidence limit?

Log in or register to post comments

Replied on Tue, 10/22/2019 - 22:38

diana Joined: 04/17/2019

Thanks for your reply!

Thanks for your reply!
The CI details are in the picture I uploaded.

File attachments

CIdetails.png

Log in or register to post comments

Replied on Wed, 10/23/2019 - 21:13

diana Joined: 04/17/2019

It seems that those two

It seems that those two methods both fails, and I have to use mxSE()to calculate CI, is that right?

Log in or register to post comments

Replied on Thu, 10/24/2019 - 10:34

AdminRobK Joined: 01/24/2014

Any luck with bootstrapping?

Any luck with bootstrapping? If not, then it still looks like you have valid CIs for everything you want except `Base.Ra[2,1]`, right? So, you'd only need to get SEs for that algebra element.

Something I forgot to mention is that, if you're going to calculate CIs from SEs, they ought to be robust SEs. I suggest the following syntax:
AceFit <- imxRobustSE(AceFit) mxSE("Base.Ra[2,1]", AceFit)
The second line will output the SE of the genetic correlation.

BTW, which method did you use to get the results shown in your attachment, CIdetails.png ?

Log in or register to post comments

Replied on Fri, 10/25/2019 - 02:06

diana Joined: 04/17/2019

1)The code getting the

1)The code getting the results shown in that attachment is as follows:

AceFitCI <- AceFit

AceFitCI$compute$steps <- list(CI=omxDefaultComputePlan(intervals=T)$steps$CI)

AceFitCI$compute$steps$CI$constraintType <- "none"

AceFitCI$compute$steps$CI$plan <- mxComputeNelderMead()

AceFitCI <- mxRun(AceFitCI)

summary(AceFitCI,verbose=T)

And I want to make sure if the CI of the Base.Ra is value±1.96*SE, is that right?

2) There exists an error when I run imxRobustSE() right after running the AceModel.
> AceFitSE <- imxRobustSE(AceFit)
Error in if (grep(pattern = model@submodels[[i]]$name, x = contributingModelNames)) { :
argument is of length zero

I also don't know why(T＿T)

Log in or register to post comments

Replied on Fri, 10/25/2019 - 11:22

AdminRobK Joined: 01/24/2014

thanks for reporting another bug

1)The code getting the results shown in that attachment is as follows:
AceFitCI <- AceFit AceFitCI$compute$steps <- list(CI=omxDefaultComputePlan(intervals=T)$steps$CI) AceFitCI$compute$steps$CI$constraintType <- "none" AceFitCI$compute$steps$CI$plan <- mxComputeNelderMead() AceFitCI <- mxRun(AceFitCI) summary(AceFitCI,verbose=T)

OK, a Nelder-Mead success story.

And I want to make sure if the CI of the Base.Ra is value±1.96*SE, is that right?

Correct (assuming you want a 95% CI).

2) There exists an error when I run imxRobustSE() right after running the AceModel.
> AceFitSE <- imxRobustSE(AceFit)
Error in if (grep(pattern = model@submodels[[i]]$name, x = contributingModelNames)) { :
argument is of length zero

That's another bug. But you can work around it by changing the fitfunction object in your MxModel, like this:
AceFit$fitfunction$groups <- c("MZ","DZ")

Log in or register to post comments

Replied on Sat, 10/26/2019 - 10:31

diana Joined: 04/17/2019

I'm sorry, do you mean run

I'm sorry, do you mean run
AceFit$fitfunction$groups <- c("MZ","DZ")
right after running
AceFit <- mxTryHardOrdinal(AceModel, intervals=T,extraTries = 15,OKstatuscodes=c(0,1,5),exhaustive=F)?

If so, running AceFit <- imxRobustSE(AceFit) still reports the same error as before.

Log in or register to post comments

Replied on Mon, 10/28/2019 - 10:13

AdminRobK Joined: 01/24/2014

In that case, when you make

In that case, when you make `AceModel`, do `mxFitFunctionMultigroup(c('MZ','DZ'))` instead. And if that doesn't help, then I must not understand the nature of this bug.

Log in or register to post comments

Replied on Mon, 10/28/2019 - 13:13

AdminRobK Joined: 01/24/2014

It turns out I did

It turns out I did misunderstand the nature of the bug. I pushed a commit that prevents the specific error you're seeing, but it still doesn't make `imxRobustSE()` work with your model. I think one possible workaround might be to put your "Base" MxModel *inside* your MZ and DZ MxModels. But, I'm not sure it's worth the trouble for just one confidence interval. You could just proceed directly to `mxSE("Base.Ra[2,1]", AceFit)`.

Edit: if you try the workaround, you will see warnings about "submodels of submodels". In your case those warnings are ignorable, since "Base" contains no data of its own.

Log in or register to post comments

Replied on Tue, 10/15/2019 - 11:46

jpritikin Joined: 05/23/2012

data?

I can execute summary(mxBootstrap(twinACEFit, 10)) against [univACEP](https://github.com/OpenMx/OpenMx/blob/master/inst/models/passing/univACEP.R). So it would be helpful if you could post [exact instructions to reproduce the problem](https://stackoverflow.com/help/minimal-reproducible-example).

Log in or register to post comments

Replied on Tue, 10/22/2019 - 04:09

ReddyM Joined: 10/02/2019

Some thougths

There are two things i observed with ML based CIs in openmx: first is the number of NAs you get depend on the starting values, even if your model is fitted, the CI estimation depends on the starting values you given. Second: the problem occurs more i think, when there is a possibility of values around zero, or either a negative or positive value can appaear in the algebra you would like to compute the CI for. Maybe it has somthing to do with possible singularities in an algebra, i dont know. For example in my cholesky models i get NA-s for the cross correlation oftne times, which is H/V or U/V etc, which has singularities where the elements of V can be zeros, so estimation around zeros could be problematic? I am not sure, since i dont know how the numeric algorithms work at all, just some guesses, and observations. One advice: try to set the starting values as the exact estimates, sometimes it can reduce, or even eliminate NAs in the CIs.
For the mods: if i said somthing misleading, or plain wrong, feel free to delete my comment, i just wanted to state some observations, maybe they can be helpful, but maybe i am in the wrong.

Log in or register to post comments

Replied on Tue, 10/22/2019 - 22:42

diana Joined: 04/17/2019

Thanks for your suggestion!

Thanks for your suggestion! I will try to change start values and see if it works.

Log in or register to post comments

Replied on Wed, 10/23/2019 - 10:27

AdminRobK Joined: 01/24/2014

please clarify?

**ReddyM**, I'm not sure I follow your post. Could you describe the specific procedure you've used to "set the starting values to the exact estimates" that leads to different results for the confidence limits? Does your procedure throw warnings about the model being changed since it was last run?

Log in or register to post comments

Replied on Wed, 10/23/2019 - 14:14

ReddyM Joined: 10/02/2019

Clarification

I run the model with some values for the starting values for the free parameters, and after it ran, i rerun the model with different starting values, it did not throw warnings, i reran it from the beginning. I did this because i wanted to see if CIs depend on the choise of starting values or not (not the actual values, that would be strange). The computed limits are the same , but this problem with NAs is not, for some starting values i got more NAs for limits, and for some starting values i got less NAs, i thought maybe if i choose my starting values really close to my estimates it will minimize the number of NAs. But to tell you the truth i didn't really continue with this experimetn, i just used bootstrapping after, so i have no idea how NAs and starting values have this relation. When i did this it kind of looked like the number of NAs are generally less, when i choose the starting values closer to the estimates of the free paramters, but i am not sure.

Log in or register to post comments

Replied on Mon, 11/11/2019 - 21:28

AdminNeale Joined: 03/01/2013

Starting values for CI search

Classic Mx - which didn’t have this problem afaik - would start CI searches at the solution, which should not vary between runs if the model is identified. If the model uses ordinal FIML there may be slight differences in the solutions due to limited numerical precision. Does CI search start at the solution?

Log in or register to post comments

Replied on Tue, 11/12/2019 - 09:53

AdminRobK Joined: 01/24/2014

yes

Does CI search start at the solution?

Yes, AFAIK. Which is why I was and remain confused by ReddyM's comment.

Log in or register to post comments

Replied on Sat, 11/30/2019 - 10:14

ReddyM Joined: 10/02/2019

I thought it starts at the

I thought it starts at the starting values. That is how i justified this behaviour, but if it starts at the solution, than i have no idea, why this happened... sorry for the confusion.

Log in or register to post comments

Replied on Fri, 08/06/2021 - 12:27

AdminRobK Joined: 01/24/2014

ReddyM, I believe I owe you

**ReddyM**, I believe I owe you an apology. In a different thread, another user shows that confidence-interval results can indeed be influenced by the choice of start values for the primary optimization, even though the different start values lead to the same point estimates at the solution, out to 6 decimal places. I reproduced that user's report and talked it over with the other developers, and we agreed that the confidence-limit search really can be sensitive to differences in the point estimates beyond 7 significant figures. So, you were indeed correct that using different start values can make the difference between a confidence limit being reported as `NA` or not.

Log in or register to post comments

Log in or register to post comments

News

how to avoid NA in confidence interval

News