CI is confusing
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?
I agree that those results
Log in or register to post comments
I am anxious to know why
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))
Log in or register to post comments
In reply to I am anxious to know why by qingwen
need more info
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`.
Log in or register to post comments
In reply to need more info by AdminRobK
Look forward to your reply
#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.**
Log in or register to post comments
I think the model might not
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?
Log in or register to post comments
In reply to I think the model might not by AdminRobK
is it acceptable that I change some path to zero to make submode
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.
Log in or register to post comments
In reply to is it acceptable that I change some path to zero to make submode by qingwen
some answers
Yes, the model is at least locally identified, at the current values of its free parameters.
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.
That was probably an artifact of model unidentification.
Log in or register to post comments
In reply to some answers by AdminRobK
Problems with the result
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.
Log in or register to post comments
In reply to Problems with the result by qingwen
result
Why is it a problem if they're negative? You haven't mentioned what trait you're studying.
I don't see that there's anything wrong in doing that.
Log in or register to post comments