CI is confusing

Posted on
No user picture. qingwen Joined: 08/15/2019
I fit a nuclear twin family model. AE model is the best. However, the CI of e makes me confused. Here are the output in free parameters and CI details section. The CI details looks strange, ranging from 0 to 2.053589e+06.

name matrix row col Estimate Std.Error A lbound ubound
1 AddGenPath MZntf.a 1 1 0.708885615 1.004657e-02
2 EnvPath MZntf.e 1 1 0.668286218 1.224575e-02
3 AMCopath MZntf.mu 1 1 0.590740014 1.836063e-02

CI details:
parameter side value fit diagnostic statusCode method
AddGenPath lower 5.859355e-01 11404.18 success OK neale-miller-1997
AddGenPath upper 8.259182e-01 11404.18 success OK neale-miller-1997
EnvPath lower -5.363198e+00 11403.31 alpha level not reached infeasible non-linear constraint neale-miller-1997
EnvPath upper 2.053589e+06 11402.99 alpha level not reached infeasible non-linear constraint neale-miller-1997

Then I add mxBootstrap(AE.NTF.Fit, 1000) to calculate CIs. Here are the results. It looks confusing that the std.error of e is 0.000000 and the CI is so close to the estimate.

name matrix. row col Estimate Std.Error 2.5% 97.5%
1 AddGenPath MZntf.a 1 1 0.708885615 0.06318805 0.56058213 0.81236401
2 EnvPath MZntf.e 1 1 0.668286218 0.00000000 0.66828622 0.66828622
3 AMCopath MZntf.mu 1 1 0.590740014 0.02719992 0.54368092 0.64904098

I would like to know what caused this. Which one is right? Is it acceptable to use CIs in bootstrap method? If not, what can I do to get the accurate CIs?

Replied on Thu, 11/21/2019 - 06:07
No user picture. qingwen Joined: 08/15/2019

Thanks for your reply. I have been working on it for many days but still could not find the solutions. Here is my script.

a <- mxMatrix(type="Lower",nrow=nv,ncol=nv,free=T,values=.6,label="AddGenPath",name="a")
m <- mxMatrix(type="Lower", nrow=nv, ncol=nv, free=T, values=.2, label="FamilialPath", name="m",)
d <- mxMatrix(type="Lower",nrow=nv,ncol=nv,free=T,values=.6,label="DomPath",name="d",)
s <- mxMatrix(type="Lower", nrow=nv, ncol=nv, free=T,values=.6,label="SibPath", name="s")
e <- mxMatrix(type="Lower",nrow=nv,ncol=nv,free=T,values=.6,label="EnvPath",name="e")
mu <- mxMatrix(type="Lower",nrow=nv,ncol=nv,free=T,values=.3,label="AMCopath",name="mu")

Vp1 <- mxMatrix(type="Full",nrow=nv,ncol=nv,free=T,values=1.5,label="VarPhen",name="Vp1")
delta1 <- mxMatrix(type="Full",nrow=nv,ncol=nv,free=T,values=.5,label="CovPhenGen",name="delta1")
q1 <- mxMatrix(type="Full",nrow=nv,ncol=nv,free=T,values=1.2,label="LatentVarAddGen",name="q1")
w1 <- mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=0, label="CovFA",name="w1",lbound=0)
x1 <- mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=0, label="VarF",name="x1")

# mxAlgebra - nonlinear constraints
Vp2 <- mxAlgebra(a %&% q1 + x1 + 2 %x% a %*% w1 + e %*% t(e) + d %*% t(d) + s %*% t(s),name="Vp2")
delta2 <- mxAlgebra(q1 %*% a + w1,name="delta2")
q2 <- mxAlgebra(1 + delta1 %&% mu,name="q2")
x2 <- mxAlgebra(2 %x% m %&% Vp1 + 2 %x% m %*% m %*% Vp1 %*% Vp1 %*% mu, name="x2")
w2 <- mxAlgebra(delta1 %*% m + delta1 %*% mu %*% Vp1 %*% t(m), name="w2")

# Equating nonlinear constraints and parameters
VpCon <- mxConstraint(Vp1==Vp2,name='VpCon')
deltaCon <- mxConstraint(delta1==delta2,name='deltaCon')
qCon <- mxConstraint(q1==q2,name='qCon')
wCon <- mxConstraint(w1==w2,name='wCon')
xCon <- mxConstraint(x1==x2,name='xCon')

# mxAlgebra - relative covariances
CvMz <- mxAlgebra(a %&% q1 + x1 + 2 %x% a %*% w1 + d %*% t(d) + s %*% t(s),name="CvMz")
CvDz <- mxAlgebra(a %&% (q1-.5) + .25 %x% d %*% t(d) + x1 + 2 %x% a %*% w1 + s %*% t(s),name="CvDz")
ParChild <- mxAlgebra(.5 %x% a %*% delta1 + .5 %x% a %*% delta1 %*% mu %*% Vp1 + m %*% Vp1 + m %*% Vp1 %&% mu, name="ParChild")
CvSps <- mxAlgebra(Vp1 %&% mu,name="CvSps")

#Put these relative covariances together into MZ relatives and DZ relatives matrices
MZ.rel.cv <- mxAlgebra(rbind(
cbind(Vp1, CvMz, ParChild, ParChild),
cbind(CvMz, Vp1, ParChild, ParChild),
cbind(ParChild, ParChild, Vp1, CvSps),
cbind(ParChild, ParChild, CvSps, Vp1)),
dimnames=list(colnames(mz.ntf),colnames(mz.ntf)),name="expCovMZRels")

DZ.rel.cv <- mxAlgebra(rbind(
cbind(Vp1, CvDz, ParChild, ParChild),
cbind(CvDz, Vp1, ParChild, ParChild),
cbind(ParChild, ParChild, Vp1, CvSps),
cbind(ParChild, ParChild, CvSps, Vp1)),
dimnames=list(colnames(dz.ntf),colnames(dz.ntf)),name="expCovDZRels")

#Put the relatives means together into means matrix (same for MZ and DZ rels)
Both.means <- mxMatrix(type="Full", nrow=1, ncol=4*nv, free=TRUE, values=.05, label="mean",dimnames=list(NULL, colnames(mz.ntf)), name="expMeanBoth")

# Put in the data
dataMZ <- mxData( observed=mz.ntf, type="raw" )
dataDZ <- mxData( observed=dz.ntf, type="raw" )

#Objective objects for two groups
objMZRels <- mxExpectationNormal(covariance="expCovMZRels",means="expMeanBoth", dimnames=selVars)
objDZRels <- mxExpectationNormal(covariance="expCovDZRels",means="expMeanBoth", dimnames=selVars)

#Combine groups
myfitfun <- mxFitFunctionML()
params <- list(a,d,s,e,mu,m,Vp1,delta1,q1,x1,w1,Vp2,delta2,q2,x2,w2,VpCon,deltaCon,qCon,xCon,wCon,CvMz,CvDz,ParChild,CvSps,MZ.rel.cv,Both.means,DZ.rel.cv,myfitfun)

modelMZ <- mxModel("MZntf",params,dataMZ,objMZRels)
modelDZ <- mxModel("DZntf",params,dataDZ,objDZRels)
obj <- mxFitFunctionMultigroup(c("MZntf","DZntf"))
ci1 <- mxCI(c("AddGenPath","SibPath","DomPath","EnvPath","AMCopath","FamilialPath","VarPhen","CovPhenGen","LatentVarAddGen", "VarF","CovFA","SibPath"))
ADSFE.NTF.Model <- mxModel("ntfASDFE",modelMZ,modelDZ,obj,ci1)

#ADFE model
ADSFE.NTF.Fit <- mxRun(ADSFE.NTF.Model,intervals=F)
ADSFE.NTF.Fit.Sum <- summary(ADSFE.NTF.Fit,verbose = F)
ADSFE.NTF.Fit.Sum

#ADFE model
ADFE.NTF.Model <- mxModel( ADSFE.NTF.Fit, name="ntfADFE")
ADFE.NTF.Model <- omxSetParameters( ADFE.NTF.Model,labels="SibPath",free=FALSE,values=0)
ADFE.NTF.Fit <- mxRun(ADFE.NTF.Model,intervals=F)
ADFE.NTF.Fit.Sum <- summary(ADFE.NTF.Fit,verbose = F)
ADFE.NTF.Fit.Sum

#ADSE model
ADSE.NTF.Model <- mxModel( ADSFE.NTF.Fit, name="ntfADSE")
ADSE.NTF.Model <- omxSetParameters( ADSE.NTF.Model,labels="FamilialPath",free=FALSE,values=0)
ADSE.NTF.Fit <- mxRun(ADSE.NTF.Model,intervals=F)
ADSE.NTF.Fit.Sum <- summary(ADSE.NTF.Fit,verbose = F)
ADSE.NTF.Fit.Sum

#ASFE model
ASFE.NTF.Model <- mxModel( ADSFE.NTF.Fit, name="ntfASFE")
ASFE.NTF.Model <- omxSetParameters( ASFE.NTF.Model,labels=c("DomPath"),free=FALSE,values=0)
ASFE.NTF.Fit <- mxRun(ASFE.NTF.Model,intervals=F)
ASFE.NTF.Fit.Sum <- summary(ASFE.NTF.Fit,verbose = F)
ASFE.NTF.Fit.Sum

#AFE model
AFE.NTF.Model <- mxModel( ADSFE.NTF.Fit, name="ntfAFE")
AFE.NTF.Model <- omxSetParameters( AFE.NTF.Model,labels=c("DomPath","SibPath"),free=FALSE,values=0)
AFE.NTF.Fit <- mxRun(AFE.NTF.Model,intervals=F)
AFE.NTF.Fit.Sum <- summary(AFE.NTF.Fit,verbose = F)
AFE.NTF.Fit.Sum

#ASE nodel
ASE.NTF.Model <- mxModel( ADSFE.NTF.Fit, name="ntfASE")
ASE.NTF.Model <- omxSetParameters( ASE.NTF.Model,labels=c("DomPath","FamilialPath"),free=FALSE,values=0)
ASE.NTF.Fit <- mxRun(ASE.NTF.Model,intervals=F)
ASE.NTF.Fit.Sum <- summary(ASE.NTF.Fit,verbose = F)
ASE.NTF.Fit.Sum

#ADE nodel
ADE.NTF.Model <- mxModel( ADSFE.NTF.Fit, name="ntfADE")
ADE.NTF.Model <- omxSetParameters( ADE.NTF.Model,labels=c("SibPath","FamilialPath"),free=FALSE,values=0)
ADE.NTF.Fit <- mxRun(ADE.NTF.Model,intervals=F)
ADE.NTF.Fit.Sum <- summary(ADE.NTF.Fit,verbose = F)
ADE.NTF.Fit.Sum

#AE model
AE.NTF.Model <- mxModel( ADSFE.NTF.Fit, name="ntfAE")
AE.NTF.Model <- omxSetParameters( AE.NTF.Model,labels=c("DomPath","SibPath","FamilialPath"),free=FALSE,values=0)
AE.NTF.Fit <- mxRun(AE.NTF.Model,intervals=T)
AE.NTF.Fit.Sum <- summary(AE.NTF.Fit,verbose = T)
AE.NTF.Fit.Sum
boot <- mxBootstrap(AE.NTF.Fit, 1000)# start with 10
summary(boot,boot.quantile=c(.025,.975))

Replied on Thu, 11/21/2019 - 10:33
Picture of user. AdminRobK Joined: 01/24/2014

In reply to by qingwen

Which optimizer are you using? If you're using CSOLNP (the on-load default), then you should try either (1) switching optimizers or (2) placing the MxConstraints in either the MZ or DZ MxModel, but not both. You can switch optimizers with mxOption(), e.g. mxOption(NULL,"Default optimizer","SLSQP") switches optimizers to SLSQP.

What do you get from qr(ADSFE.NTF.Fit$output$constraintJacobian)?

What do you get from mxVersion()?

Aside from the confidence intervals...do the point estimates you get look reasonable?

How many phenotypes are you analyzing, anyhow? Your posted script doesn't show what value is assigned to `nv`.

Replied on Fri, 11/22/2019 - 00:16
No user picture. qingwen Joined: 08/15/2019

In reply to by AdminRobK

__Thanks for your quick reply. My openMx version is 2.14.11. I am analyzing one phenotype and nv<-1 in my script. Aside from the CIs, the point estimate looks reasonable. The estimate of E is 0.662516155 and the std.error is 1.503494e-02.__

#The output I get from `qr(ADSFE.NTF.Fit$output$constraintJacobian)`

$qr
AddGenPath DomPath AMCopath FamilialPath VarPhen SibPath
MZntf.VpCon[1,1] 2.0110369 1.0050667 0.000000000 0.00000000 -9.241115e-01 5.415512e-01
MZntf.deltaCon[1,1] 0.5352611 -1.1643035 0.000000000 0.00000000 1.070522e+00 -6.273514e-01
MZntf.qCon[1,1] 0.0000000 0.0000000 0.305701843 0.04914830 3.202266e-04 0.000000e+00
MZntf.xCon[1,1] 0.0000000 0.0000000 0.003955549 0.95232837 6.771227e-03 0.000000e+00
MZntf.wCon[1,1] 0.0000000 0.0000000 0.037370059 0.66355074 -2.248940e-03 1.233176e-29
DZntf.VpCon[1,1] 0.4620558 -0.3661017 0.000000000 0.00000000 7.404974e-14 1.755417e-16
DZntf.deltaCon[1,1] 0.5352611 0.6580158 0.000000000 0.00000000 -4.936649e-14 -9.486833e-01
DZntf.qCon[1,1] 0.0000000 0.0000000 0.706107522 -0.02135926 -1.677573e-02 7.375195e-16
DZntf.xCon[1,1] 0.0000000 0.0000000 0.003955549 0.24386946 -5.335934e-01 2.345862e-14
DZntf.wCon[1,1] 0.0000000 0.0000000 0.037370059 0.66355074 5.979116e-01 -2.628627e-14
EnvPath CovPhenGen LatentVarAddGen VarF CovFA mean
MZntf.VpCon[1,1] 1.224478 -1.070522e+00 6.080053e-01 9.241115e-01 1.842139e+00 0
MZntf.deltaCon[1,1] -1.418477 -9.241115e-01 1.992176e-01 -1.070522e+00 3.024421e-02 0
MZntf.qCon[1,1] 0.000000 4.670787e-01 -1.412215e+00 -7.911097e-03 -7.474012e-02 0
MZntf.xCon[1,1] 0.000000 1.940195e-02 7.288240e-02 -4.875699e-01 -1.325505e+00 0
MZntf.wCon[1,1] 0.000000 9.909658e-03 1.835284e-02 1.327484e+00 -4.872886e-01 0
DZntf.VpCon[1,1] 0.000000 -5.206094e-16 -1.337556e-15 -9.300964e-14 3.421490e-14 0
DZntf.deltaCon[1,1] 0.000000 3.659091e-16 1.059826e-15 7.249687e-14 -2.655185e-14 0
DZntf.qCon[1,1] 0.000000 -1.479618e-17 -1.024771e-17 -2.017033e-16 1.122347e-16 0
DZntf.xCon[1,1] 0.000000 3.517239e-01 -9.219073e-18 -4.955415e-16 2.676370e-16 0
DZntf.wCon[1,1] 0.000000 5.862064e-02 1.338068e-01 -1.855310e-16 -7.219559e-17 0

$rank
[1] 5

$qraux
[1] 1.462056e+00 1.658016e+00 1.706108e+00 1.243869e+00 1.597912e+00 1.316228e+00 0.000000e+00
[8] 1.934267e+00 1.991007e+00 7.249711e-14 2.655195e-14 0.000000e+00

$pivot
[1] 1 2 5 6 7 3 4 8 9 10 11 12

attr(,"class")
[1] "qr"

#The output of free parameters in AE model

free parameters:
name matrix row col Estimate Std.Error A lbound ubound
1 AddGenPath MZntf.a 1 1 0.786053523 1.821280e-02
2 EnvPath MZntf.e 1 1 0.662516155 1.503494e-02 !
3 AMCopath MZntf.mu 1 1 0.354074786 2.500254e-02
4 VarPhen MZntf.Vp1 1 1 1.023611563 2.421320e-02
5 CovPhenGen MZntf.delta1 1 1 0.293082083 1.894200e-02
6 LatentVarAddGen MZntf.q1 1 1 0.875558820 4.474109e-03
7 VarF MZntf.x1 1 1 0.001234334 1.503001e-18
8 CovFA MZntf.w1 1 1 0.027007610 8.958187e-18 0!
9 mean MZntf.expMeanBoth 1 1 0.015469679 2.003104e-02

#The beginning of the mentioned script

setwd(".....")
getwd()
mxOption(NULL,"Default optimizer","SLSQP")
data <- read.spss("data all.sav", to.data.frame = TRUE)
selVars = c("t1self","t2self","faself","moself")
mz.ntf <- subset(data, zyg2012=="MZ", selVars)
dz.ntf <- subset(data, zyg2012=="DZ", selVars)
nv<-1

#Here is the output for mxVersion().

OpenMx version: 2.14.11 [GIT v2.14.11]
R version: R version 3.5.3 (2019-03-11)
Platform: x86_64-apple-darwin15.6.0
MacOS: 10.14.6
Default optimizer: SLSQP
NPSOL-enabled?: No
OpenMP-enabled?: Yes

**Finally, I tried to place the MxConstraints only in MZ model and in DZ model by the following correction. These two results are consistent and the CIs could be estimated. However the result changes a lot compared to placing the MxConstraints both MZ and DZ MxModel. The best model is ADSE model and the CIs of a and s are including zero.**
#The code I changed

params <-list(a,d,s,e,mu,m,Vp1,delta1,q1,x1,w1,Vp2,delta2,q2,x2,w2,VpCon,deltaCon,qCon,xCon,wCon,CvMz,CvDz,ParChild,CvSps,MZ.rel.cv,Both.means,DZ.rel.cv,myfitfun)
params1 <- list(a,d,s,e,mu,m,Vp1,delta1,q1,x1,w1,Vp2,delta2,q2,x2,w2,CvMz,CvDz,ParChild,CvSps,MZ.rel.cv,Both.means,DZ.rel.cv,myfitfun)
modelMZ <- mxModel("MZntf",params,dataMZ,objMZRels)
modelDZ <- mxModel("DZntf",params1,dataDZ,objDZRels)

#the results of ADSE model

ree parameters:
name matrix row col Estimate Std.Error A lbound ubound
1 AddGenPath MZntf.a 1 1 0.430265872 2.218452e-03 !
2 DomPath MZntf.d 1 1 0.534466963 2.250038e-02 !
3 SibPath MZntf.s 1 1 0.286951652 2.521712e-02 !
4 EnvPath MZntf.e 1 1 0.656287839 1.330449e-02
5 AMCopath MZntf.mu 1 1 0.353079233 1.434644e-02 !
6 VarPhen MZntf.Vp1 1 1 1.014386192 NA !
7 CovPhenGen MZntf.delta1 1 1 0.480431662 NA !
8 LatentVarAddGen MZntf.q1 1 1 1.081422035 NA !
9 VarF MZntf.x1 1 1 0.002789132 1.241243e-18 !
10 CovFA MZntf.w1 1 1 0.015196440 2.338501e-18 ! 0!
11 mean MZntf.expMeanBoth 1 1 0.014373885 1.837673e-02 !

CI details:
parameter side value fit diagnostic
1 AddGenPath lower -5.198694e-01 11930.47 active box constraint
2 AddGenPath upper 4.701160e-01 11946.02 alpha level not reached
3 SibPath lower -4.603577e-01 11930.43 active box constraint
4 SibPath upper 3.198604e-01 11940.42 alpha level not reached
5 DomPath lower 4.848444e-01 11937.81 alpha level not reached
6 DomPath upper 6.464772e-01 11930.45 active box constraint
7 EnvPath lower 6.285536e-01 11930.44 active box constraint
8 EnvPath upper 6.994059e-01 11930.43 active box constraint
9 AMCopath lower 3.858827e-01 11938.80 alpha level not reached
10 AMCopath upper 4.055172e-01 11930.43 active box constraint
11 FamilialPath lower 0.000000e+00 11963.10 alpha level not reached
12 FamilialPath upper 0.000000e+00 11930.43 active box constraint
13 VarPhen lower 9.809736e-01 11945.73 alpha level not reached
14 VarPhen upper 1.127305e+00 11941.67 alpha level not reached
15 CovPhenGen lower -5.814936e-01 11930.43 active box constraint
16 CovPhenGen upper 7.634307e-01 11989.53 alpha level not reached
17 LatentVarAddGen lower 1.024691e+00 11954.81 alpha level not reached
18 LatentVarAddGen upper 1.079775e+00 11941.99 alpha level not reached
19 VarF lower 5.837377e-22 11930.43 active box constraint
20 VarF upper -4.168990e-22 11930.43 active box constraint
21 CovFA upper 7.623297e-21 11930.43 active box constraint
22 CovFA lower 0.000000e+00 11925.85 success

**Why are these two results so different? I am looking forward to your reply. Thanks forward.**

Replied on Fri, 11/22/2019 - 11:13
Picture of user. AdminRobK Joined: 01/24/2014

I think the model might not be identified. What do you get if you run mxCheckIdentification() on your MxModel objects?

Also, I notice that you're not adjusting for age of parents and twins, nor for sex of twins. Doing so is generally advised. But, worry about that after trying `mxCheckIdentification()`.

For that matter, your script also models parents and twins as having the same model-expected mean:

Both.means <- mxMatrix(type="Full", nrow=1, ncol=4*nv, free=TRUE, values=.05, label="mean",dimnames=list(NULL, colnames(mz.ntf)), name="expMeanBoth")

Is that really realistic for the trait you're studying?

Replied on Sat, 11/23/2019 - 01:29
No user picture. qingwen Joined: 08/15/2019

In reply to by AdminRobK

You are right. The ADSFE model could not be identified. However, the submodels seemed to be identified. The output for the submodels is just like this. The status is “TRUE” and $non_identified_parameters is "None". Could it be seen as identified? If it is identified, is it acceptable that I set some path to zero to make submodels on a non-identified model? If not, should I model ADFE, ADSE and ASFE separately?

Also, the question I mentioned earlier still make me confused. Why I move constraints to MZmodel make such a difference?

You are right that it seems not feasible to model parents and twins as having the same model-expected mean. However it is always to do so in nuclear twin family study. I am new in this area and also have some confusion. I will consider it seriously and control sex and age factors after dealing with the identification problem.

I am looking forward to your reply.

Replied on Tue, 11/26/2019 - 16:54
Picture of user. AdminRobK Joined: 01/24/2014

In reply to by qingwen

The status is “TRUE” and $non_identified_parameters is "None". Could it be seen as identified?

Yes, the model is at least locally identified, at the current values of its free parameters.

If it is identified, is it acceptable that I set some path to zero to make submodels on a non-identified model? If not, should I model ADFE, ADSE and ASFE separately?

Yes, that sounds reasonable to me. It's analogous to how, in the classical twin design, the *ACDE* model is unidentified, and therefore one must fix either the *C* or *D* component to zero in any given model.

Also, the question I mentioned earlier still make me confused. Why I move constraints to MZmodel make such a difference?

That was probably an artifact of model unidentification.

Replied on Thu, 11/28/2019 - 02:24
No user picture. qingwen Joined: 08/15/2019

In reply to by AdminRobK

Thanks for your reply. I still have some problems with the result. The path diagram of nuclear twin model is in the attachment.

Firstly, the estimates of m and w are significantly, but NEGATIVE. It seems difficult to explain it reasonably. I don’t know if I should report these negative values as the output shows or I should set these factors between 0 and 1 ahead of running the script.

Secondly, I tried to make ADFE, ASFE and ADSE model separately. However, the ADSE model do not have an equal degree of freedom as ASFE and ADFE model because ADSE model has no estimates of m, w and some nonlinear constraints (i.e., w1=w2, x1=x2). Therefore, I am not sure if it is acceptable for me to do this.

Your reply is highly appreciated.

Replied on Tue, 12/03/2019 - 15:03
Picture of user. AdminRobK Joined: 01/24/2014

In reply to by qingwen

Firstly, the estimates of m and w are significantly, but NEGATIVE. It seems difficult to explain it reasonably. I don’t know if I should report these negative values as the output shows or I should set these factors between 0 and 1 ahead of running the script.

Why is it a problem if they're negative? You haven't mentioned what trait you're studying.

Secondly, I tried to make ADFE, ASFE and ADSE model separately. However, the ADSE model do not have an equal degree of freedom as ASFE and ADFE model because ADSE model has no estimates of m, w and some nonlinear constraints (i.e., w1=w2, x1=x2). Therefore, I am not sure if it is acceptable for me to do this.

I don't see that there's anything wrong in doing that.