# how to avoid NA in confidence interval

Attachment | Size |
---|---|

3.png | 33.17 KB |

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

## suggestions

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

In reply to suggestions by AdminRobK

## error of running mxBootstrap()

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

In reply to error of running mxBootstrap() by diana

## traceback()

Log in or register to post comments

In reply to traceback() by AdminRobK

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

In reply to hello, this was the result by diana

## not helpful

Log in or register to post comments

In reply to hello, this was the result by diana

## mxVersion?

diana, what is your `mxVersion()` output?Log in or register to post comments

In reply to mxVersion? by AdminRobK

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

In reply to > mxVersion() by diana

## update and try again

Log in or register to post comments

In reply to update and try again by AdminRobK

## I updated OpenMx to the

> 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

In reply to I updated OpenMx to the by diana

## how to reproduce?

Log in or register to post comments

In reply to how to reproduce? by jpritikin

## Thank you very much. I am

Log in or register to post comments

In reply to how to reproduce? by jpritikin

## Hello! I am so sorry for

`mxGenerateData()`

was uploaded and the code was paseted above.Thank you very much for your help!

Log in or register to post comments

In reply to I updated OpenMx to the by diana

## mxGenerateData()

`mxGenerateData()`

may be useful for the purpose.Log in or register to post comments

In reply to mxGenerateData() by AdminRobK

## Thank you

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

Log in or register to post comments

In reply to Thank you by diana

## fixed

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

testBoot <- mxBootstrap(AceBoot,10)

Thanks for your help.

Log in or register to post comments

In reply to fixed by jpritikin

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

In reply to There still exists an error by diana

## mxCI

Log in or register to post comments

In reply to mxCI by jpritikin

## Here are something confusing

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

In reply to Here are something confusing by diana

## some answers

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

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.

Absolutely not. If you want 95% CIs, you need to be using

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

Log in or register to post comments

In reply to some answers by AdminRobK

## Thanks for your reply!

The CI details are in the picture I uploaded.

Log in or register to post comments

In reply to some answers by AdminRobK

## It seems that those two

`mxSE()`

to calculate CI, is that right?Log in or register to post comments

In reply to It seems that those two by diana

## Any luck with bootstrapping?

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

In reply to Any luck with bootstrapping? by AdminRobK

## 1)The code getting the

`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

In reply to 1)The code getting the by diana

## thanks for reporting another bug

OK, a Nelder-Mead success story.

Correct (assuming you want a 95% CI).

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

In reply to thanks for reporting another bug by AdminRobK

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

In reply to I'm sorry, do you mean run by diana

## In that case, when you make

Log in or register to post comments

In reply to I'm sorry, do you mean run by diana

## It turns out I did

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

## data?

`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

## Some thougths

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

In reply to Some thougths by ReddyM

## Thanks for your suggestion!

Log in or register to post comments

In reply to Some thougths by ReddyM

## please clarify?

Log in or register to post comments

In reply to please clarify? by AdminRobK

## Clarification

Log in or register to post comments

In reply to Clarification by ReddyM

## Starting values for CI search

Log in or register to post comments

In reply to Starting values for CI search by AdminNeale

## yes

Yes, AFAIK. Which is why I was and remain confused by

ReddyM's comment.Log in or register to post comments

In reply to yes by AdminRobK

## I thought it starts at the

Log in or register to post comments

In reply to I thought it starts at the by ReddyM

## ReddyM, I believe I owe you

Log in or register to post comments