Bivariate twin model(ACEorADE)

Posted on
No user picture. Fengxia Joined: 11/02/2017
Hello,

Recently, I have done a bivariate models.BUT I am confused with how to choose ACE or ADE in the bivariate models? I have tried using ICC in each variable apartly, BUT the fitting model one is ACE when another is ADE. NOW the question is that which one bivariate model I should use ? ACE-ACE or ADE-ADE or ACE-ADE? If the the answer is ACE-ADE, how to write the script?

Thanks for some help!

Cheers,
Fengxia

Replied on Wed, 03/21/2018 - 10:20
Picture of user. AdminRobK Joined: 01/24/2014

I'm guessing this is a follow-up to your thread about your triphenotype analysis of sleep duration, diabetes, and HbA1c, right? I really don't think it's wise to divide your sample into two overlapping subsamples and run separate analyses on each. You can incorporate the U-shaped trend into your triphenotype model by creating a dummy-coded covariate (a "definition variable") for sleep duration--make the 6-8h group the reference group, and create one indicator variable for being in the <6h group and another for being in the >8h group. Condition the liability-scale mean of diabetes on these dummy variables, and fix the threshold to zero to identify the model. If the regression coefficients for the two indicator variables are both positive and differ significantly from zero, then you have results consistent with your logistic regression analyses.
Replied on Mon, 03/26/2018 - 08:29
No user picture. Fengxia Joined: 11/02/2017

In reply to by AdminRobK

Thanks for your suggestion!

YES, but my analysis is biphenotype, that's sleep duration and diabetes, sleep duration and HbA1c separately. According to my understanding about "definition variable" for sleep duration, I should define sleep1 with 6-8h=0 and <6h=1, sleep2 with sleep1 with 6-8h=0 and >8h=1? BUT there will have missing data in sleep1 when sleep=>8h, and sleep2 when sleep=>8h, then I can't impute all my sample. How can I solve this problem?

If I have sleep1 and sleep2, the model would be sleep1-sleep2-diabetes ?

Finally, I really hope to receive your suggestion about ACE-ACE or ADE-ADE or ACE-ADE?

Replied on Mon, 04/02/2018 - 13:56
Picture of user. AdminRobK Joined: 01/24/2014

In reply to by Fengxia

Are you not familiar with dummy coding? You'd have two indicator variables: one for <6h ('sleep1', say), and one for >8h ('sleep2', say). So, 6-8h would be scored (0,0), <6h would be scored (1,0), and >8h would be scored (0,1). You would condition either the mean or the threshold of diabetes on the two variables via regression, in addition to how you're currently conditioning on age.

The thing is, though, now that I think about it, I don't think the resulting model would be identified. It would be sort of like regressing diabetes onto sleep duration while simultaneously estimating a covariance between the two phenotypes. Biometrically decomposing a U-shaped relationship between two ordinal phenotypes doesn't seem so easy to me now.

Finally, I really hope to receive your suggestion about ACE-ACE or ADE-ADE or ACE-ADE?

If all three of those models seem plausible to you, why not try all three, and compare their fit to the data? Try adapting a script for your purposes, and post again for help if you get stuck. Note that with the ACE-ADE model, there would be only one free path from C, going to the first phenotype, and only one free path from D, going to the second phenotype. Thus, only A and E would contribute to the phenotypic covariance, and only A would contribute to the cross-twin cross-trait covariance.

Replied on Thu, 10/11/2018 - 10:15
No user picture. Fengxia Joined: 11/02/2017

In reply to by AdminRobK

I tried to adapt ACE-ADE model, but I have no idea about how to set one ACE and another ADE. When run ACE-ACE, we use:
<% # Algebra for expected Variance/Covariance Matrices in MZ & DZ twins
covMZ <- mxAlgebra( expression= rbind( cbind(A+C+E , A+C),
cbind(A+C , A+C+E)) , name="ExpCovMZ" )
covDZ <- mxAlgebra( expression= rbind( cbind(A+C+E , 0.5%x%A+C),
cbind(0.5%x%A+C , A+C+E)), name="ExpCovDZ" ) %>

how to change it when need to run an ACE-ADE?

Replied on Mon, 10/15/2018 - 23:51
No user picture. Fengxia Joined: 11/02/2017

In reply to by AdminRobK

I just run the ACE-ACE model, and because I don't know how to run with dummy coding, I divided the sample into two subsample (one include sleep duration <6h and 6~8h,another >8h and 6~8h). My script in <6h and 6~8h subsample:

<%rm(list=ls())
#ls()

require(OpenMx)
require(psych)

# change optimizer ########
### COMMENT RALF: I faced some issues in finding confidence intervals using CSOLNP, but it worked fine with SLSQP.
#mxOption(NULL, "Default optimizer", "NPSOL")
mxOption( NULL, 'Default optimizer' ,'SLSQP' )
#mxOption(NULL, "Default optimizer", "CSOLNP")

#Bivdata <- read.table ('ADHDdia_bage.csv', header=T, sep=",", na.strings="NA")
Bivdata <- read.csv('b1_2.csv', header=T, sep=',')

Bivdata <- Bivdata[Bivdata$zygsex!=3,]

names(Bivdata)
summary(Bivdata)
str(Bivdata)
describe(Bivdata)

#Make the integer variables ordered factors (Note:because of the 'cuting' above, this is not necessary here)
Bivdata$sleepd <-mxFactor(Bivdata$sleepd, levels=c(0:1) )
Bivdata$sleepd_2 <-mxFactor(Bivdata$sleepd_2, levels=c(0:1) )
Bivdata$dia_b <- mxFactor(Bivdata$dia_b, levels= c(0:1))
Bivdata$dia_b2 <- mxFactor(Bivdata$dia_b2, levels= c(0:1))

mzData <- subset(Bivdata, zygosity==1, c(dia_b,dia_b2,sleepd,sleepd_2,age,age2,SEX,SEX2))
dzData <- subset(Bivdata, zygosity==2, c(dia_b,dia_b2,sleepd,sleepd_2,age,age2,SEX,SEX2))

mzData$sleepd1<-mzData$sleepd
mzData$sleepd2<-mzData$sleepd_2
dzData$sleepd1<-dzData$sleepd
dzData$sleepd2<-dzData$sleepd_2

mzData$dia_b1<-mzData$dia_b
dzData$dia_b1<-dzData$dia_b

vars <-c('dia_b','sleepd')
selVars <-c('dia_b1', 'sleepd1', 'dia_b2', 'sleepd2' )
useVars <-c('dia_b1', 'sleepd1', 'dia_b2', 'sleepd2' , 'age', 'age2','SEX','SEX2')

#Make the integer variables ordered factors (Note:because of the 'cuting' above, this is not necessary here)
sum(is.na(mzData))

sum(is.na(dzData))

mzData <- na.omit(mzData)

describe(mzData)
describe(dzData)

# To get the CTs for dia_b
table(mzData$dia_b1, mzData$dia_b2)
table(dzData$dia_b1, dzData$dia_b2)
# To get the CTs for ADHD
table(mzData$sleepd1, mzData$sleepd2)
table(dzData$sleepd1, dzData$sleepd2)

rMZ <- cor(mzData,use="complete")
rDZ <- cor(dzData,use="complete")
rMZ
rDZ
cMZ <- cov(mzData,use="complete")
cDZ <- cov(dzData,use="complete")
cMZ
cDZ

### COMMENT RALF: You can quickly get (unadjusted) tetrachoric correlations using the following library and function:
library( polycor )
polychor( table(mzData$dia_b1, mzData$dia_b2) , std.err = T , ML=T )
polychor( table(dzData$dia_b1, dzData$dia_b2) , std.err = T , ML=T )

polychor( table(dzData$sleepd, dzData$sleepd2) , std.err = T , ML=T )
polychor( table(dzData$sleepd, dzData$sleepd2) , std.err = T , ML=T )

polychor( table(mzData$dia_b1, mzData$sleepd1) , std.err = T , ML=T )
polychor( table(dzData$dia_b1, dzData$sleepd1) , std.err = T , ML=T )

polychor( table(mzData$dia_b1, mzData$sleepd2) , std.err = T , ML=T )
polychor( table(dzData$dia_b1, dzData$sleepd2) , std.err = T , ML=T )

#sink('sleepdia_bbivsen916.doc')

nv <- 2 # number of variables per twin
ntv <- nv*2 # number of variables per pair
nth <- 1 # number of max thresholds

# 1) Fits a constrained Polychoric correlation model
# TH same across twins but different across zyg groups
# Age effect is different acros variables, but same across thresholds within variables (if c>2)
# There is one overall rPH between var1-2 and the x-trait x-twin correlations are symmetric
# ------------------------------------------------------------------------------------------------------------------------------

# CREATE LABELS & START VALUES as objects(to ease specification)
#LabThMZ <-c('Tmz_11','Tmz_21','imz_21') # THs for var 1 and 2 for a twin individual (mz)
LabCorMZ <-c('r21','rMZ1','MZxtxt','MZxtxt','rMZ2','r21')
#LabThDZ <-c('Tdz_11','idz_11','Tdz_21','idz_21') # THs for var 1 and 2 for a twin individual (mz)
LabCorDZ <-c('r21','rDZ1','DZxtxt','DZxtxt','rDZ2','r21')
#LabCov <-c('BageThdia_b', 'BageThdia_b', 'BageThadhd', 'BageThadhd')
ThPat <-c(T,T)
#StCorMZ <-c(-.3, .7, .7, -.3)
#StCorDZ <-c(-.3, .4, .3, -.3)
StTH <-c(1.1, -1.1 )

# CREATE LABELS & START VALUES as objects(to ease specification)
LabThMZ <-c('Tmz_11','Tmz_21') # THs for var 1 and 2 for a twin individual (mz)
#LabCorMZ <-c('r21','rMZ1','rMZ2','r21')
LabThDZ <-c('Tdz_11','Tdz_21') # THs for var 1 and 2 for a twin individual (mz)
#LabCorDZ <-c('r21','rDZ1','rDZ2','r21')
LabCov <-c('BageThdia', 'BageThdia', 'BageThsleep', 'BageThsleep')
#ThPat <-c(T,F,T,T)
StCorMZ <-c(-.3, .7, -.3,-.3, .7, -.3)
StCorDZ <-c(-.3, .4, -.15,-.15, .3, -.3)
#tTH <-c(1.5, 0, -.8, 1.4 )

# Define definition variables to hold the Covariates
obsAge1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=F, labels=c("data.age"), name="Age1")
obsAge2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=F, labels=c("data.age2"), name="Age2")

obsSex1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=F, labels=c("data.SEX"), name="Sex1")
obsSex2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=F, labels=c("data.SEX2"), name="Sex2")
# Matrix & Algebra for expected means (SND), Thresholds, effect of Age on Th and correlations
Mean <-mxMatrix( type="Zero", nrow=1, ncol=ntv, name="M" )
betaA <-mxMatrix( type="Full", nrow=nth, ncol=nv, free=T, values=.2, labels=LabCov, name="BageTH" )
betaS <-mxMatrix( type="Full", nrow=nth, ncol=nv, free=T, values=.2, labels=LabCov, name="BsexTH" )

Tmz <-mxMatrix( type="Full", nrow=nth, ncol=nv, free=ThPat, values=StTH, lbound=c(-3,.001), ubound=3, labels=LabThMZ, name="ThMZ")
#inc <-mxMatrix( type="Lower",nrow=nth, ncol=nth, free=F, values=1, name="Low")
ThresMZ <-mxAlgebra( expression= cbind(ThMZ + BageTH%x%Age1+BsexTH%x%Sex1, ThMZ + BageTH%x%Age2+BsexTH%x%Sex2), name="expThresMZ")
CorMZ <-mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=StCorMZ, labels=LabCorMZ, lbound=-.99, ubound=.99, name="expCorMZ")

Tdz <-mxMatrix( type="Full", nrow=nth, ncol=nv, free=ThPat, values=StTH, lbound=c(-3,.001), ubound=3, labels=LabThDZ, name="ThDZ")
ThresDZ <-mxAlgebra( expression= cbind(ThDZ + BageTH%x%Age1+BsexTH%x%Sex1, ThDZ + BageTH%x%Age2+BsexTH%x%Sex2), name="expThresDZ")
CorDZ <-mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=StCorDZ, labels=LabCorDZ, lbound=-.99, ubound=.99, name="expCorDZ")

# Data objects for Multiple Groups
dataMZ <- mxData( observed=mzData, type="raw" )
dataDZ <- mxData( observed=dzData, type="raw" )

# Objective objects for Multiple Groups
objMZ <- mxFIMLObjective( covariance="expCorMZ", means="M", dimnames=selVars, thresholds="expThresMZ" )
objDZ <- mxFIMLObjective( covariance="expCorDZ", means="M", dimnames=selVars, thresholds="expThresDZ" )

# Combine Groups
modelMZ <- mxModel( obsAge1, obsAge2, obsSex1,obsSex2,Mean,betaA, betaS, Tmz, ThresMZ, CorMZ, dataMZ, objMZ, name="MZ" )
modelDZ <- mxModel( obsAge1, obsAge2,obsSex1,obsSex2, Mean,betaA, betaS, Tdz, ThresDZ, CorDZ, dataDZ, objDZ, name="DZ" )
minus2ll <- mxAlgebra( expression=MZ.objective + DZ.objective, name="m2LL" )
obj <- mxAlgebraObjective( "m2LL" )
Conf <- mxCI (c ('MZ.expCorMZ[2,1]', 'MZ.BageTH[1,1]', 'MZ.BageTH[1,2]') )
SatModel <- mxModel( "Sat", modelMZ, modelDZ, minus2ll, obj, Conf )

# -------------------------------------------------------------------------------------------------------------------------------
# 1) RUN Saturated Model

SatFit <- mxRun(SatModel, intervals=TRUE)
(SatSumm <- summary(SatFit))

# Generate SatModel Output

round(SatFit@output$estimate,4)

SatFit$MZ.expCorMZ
SatFit$DZ.expCorDZ

(ExpTetraCorMZ <-mxEval(MZ.expCorMZ, SatFit))
(ExpTetraCorDZ <-mxEval(DZ.expCorDZ, SatFit))

# ---------------------------------------------------------------------------------------------------
# 2) Specify and Run Sub1Model: equating Thresholds across MZ and DZ group

Sub1Model <- mxModel(SatModel, name="sub1")
Sub1Model <- omxSetParameters(Sub1Model, labels=LabThDZ, newlabels=LabThMZ, free=ThPat, values=StTH, lbound=c(-3,.001), ubound=3 )
Sub1Fit <- mxRun(Sub1Model, intervals=TRUE)
(Sub1Summ <- summary(Sub1Fit))
round(Sub1Fit@output$estimate,4)

# Generate Sub1Model Output
(ExpTetraCorMZ <-mxEval(MZ.expCorMZ, Sub1Fit))
(ExpTetraCorDZ <-mxEval(DZ.expCorDZ, Sub1Fit))
mxCompare(SatFit,Sub1Fit)

# CREATE LABELS for Cholesky decomposition with a function
laLower <- function(la,nv) { paste(la,rev(nv+1-sequence(1:nv)),rep(1:nv,nv:1),sep="_") }
LabTh <-c('T_11','i_11','T_12','i_12')

# 3) Specify ACE Model, with ONE overall set of Thresholds with Age effects on thresholds
# ------------------------------------------------------------------------------------------

# CREATE LABELS for Cholesky decomposition with a function
laLower <- function(la,nv) { paste(la,rev(nv+1-sequence(1:nv)),rep(1:nv,nv:1),sep="_") }
LabTh <-c('T_11','i_11','T_12','i_12')

# Define definition variables to hold the Covariates
obsAge1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=F, labels=c("data.age"), name="Age1")
obsAge2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=F, labels=c("data.age2"), name="Age2")
obsSex1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=F, labels=c("data.SEX"), name="Sex1")
obsSex2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=F, labels=c("data.SEX2"), name="Sex2")
# Matrix & Algebra for expected means (SND), Thresholds, effect of Age on Th and correlations
Mean <-mxMatrix( type="Zero", nrow=1, ncol=ntv, name="M" )
betaA <-mxMatrix( type="Full", nrow=nth, ncol=nv, free=T, values=.2, labels=LabCov, name="BageTH" )
betaS <-mxMatrix( type="Full", nrow=nth, ncol=nv, free=T, values=.2, labels=LabCov, name="BsexTH" )

Tr <-mxMatrix( type="Full", nrow=nth, ncol=nv, free=ThPat, values=StTH, labels=LabTh, lbound=c(-3,.001), ubound=3, name="Th")
#inc <-mxMatrix( type="Lower",nrow=nth, ncol=nth, free=F, values=1, name="Low")
Thres <-mxAlgebra( expression= cbind(Th + BageTH%x%Age1+BsexTH%x%Sex1, Th + BageTH%x%Age2+BsexTH%x%Sex2), name="expThres")

# Matrices declared to store a, c, and e Path Coefficients
pathA <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=c(.8,-.3,.6), label=laLower("a",nv), name="a" )
pathC <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=c(.3,-.1,.3), label=laLower("c",nv), name="c" )
pathE <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=c(.8,-.3,.6), label=laLower("e",nv), name="e" )

# Matrices generated to hold A, C, and E components and total Variance
covA <- mxAlgebra( expression=a %*% t(a), name="A" )
covC <- mxAlgebra( expression=c %*% t(c), name="C" )
covE <- mxAlgebra( expression=e %*% t(e), name="E" )
covP <- mxAlgebra( expression=A+C+E, name="V" )

# Algebra to compute Phenotypic, A, C & E correlations
# Algebra to compute total variances and standard deviations (dia_bgonal only)

matI <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I")
invSD <- mxAlgebra( expression=solve(sqrt(I*V)), name="iSD")
rph <- mxAlgebra( expression= solve(sqrt(I*V)) %*% V %*% solve(sqrt(I*V)), name="Rph")
rA <- mxAlgebra( expression= solve(sqrt(I*A)) %*% A %*% solve(sqrt(I*A)), name="Ra" )
rC <- mxAlgebra( expression= solve(sqrt(I*C)) %*% C %*% solve(sqrt(I*C)), name="Rc" )
rE <- mxAlgebra( expression= solve(sqrt(I*E)) %*% E %*% solve(sqrt(I*E)), name="Re" )

h2 <- mxAlgebra( expression=A/V, name="h2")
c2 <- mxAlgebra( expression=C/V, name="c2")
e2 <- mxAlgebra( expression=E/V, name="e2")

# Algebra to compute standardized variance components
covP <- mxAlgebra( expression=A+C+E, name="V" )
StA <- mxAlgebra( expression=A/V, name="h2")
StC <- mxAlgebra( expression=C/V, name="c2")
StE <- mxAlgebra( expression=E/V, name="e2")

#std <- mxAlgebra(expression=cbind(iSD %*% a,iSD %*% c,iSD %*% e),name="std")
stdA <- mxAlgebra(expression=iSD %*% a,name="stda")
stdC <- mxAlgebra(expression=iSD %*% c,name="stdc")
stdE <- mxAlgebra(expression=iSD %*% e,name="stde")

#std <- mxAlgebra(expression=cbind(iSD %*% a,iSD %*% c,iSD %*% e),name="std")

# Algebra for expected Variance/Covariance Matrices in MZ & DZ twins
covMZ <- mxAlgebra( expression= rbind( cbind(A+C+E , A+C),
cbind(A+C , A+C+E)), name="expCovMZ" )
covDZ <- mxAlgebra( expression= rbind( cbind(A+C+E , 0.5%x%A+C),
cbind(0.5%x%A+C , A+C+E)), name="expCovDZ" )

# Constraint on total variance of Ordinal variables (A+C+E=1)
mUnv <- mxMatrix( type="Unit", nrow=nv, ncol=1, name="Unv" )
varL <- mxConstraint( expression=diag2vec(V)==Unv, name="VarL" )

# Data objects for Multiple Groups
dataMZ <- mxData( observed=mzData, type="raw" )
dataDZ <- mxData( observed=dzData, type="raw" )

# Objective objects for Multiple Groups
objMZ <- mxFIMLObjective( covariance="expCovMZ", means="M", dimnames=selVars, thresholds="expThres" )
objDZ <- mxFIMLObjective( covariance="expCovDZ", means="M", dimnames=selVars, thresholds="expThres" )

# Combine Groups

pars1 <-list(obsAge1, obsAge2, obsSex1,obsSex2,Mean, betaA,betaS, Tr, Thres,matI, invSD,mUnv, varL)
pars2 <-list(pathA, pathC, pathE, covA, covC, covE, covP, StA, StC, StE, matI, rph, rA, rC, rE,invSD,stdA, stdC,stdE )

modelMZ <- mxModel( pars1, pars2, covMZ, dataMZ, objMZ, name="MZ" )
modelDZ <- mxModel( pars1,pars2, covDZ, dataDZ, objDZ, name="DZ" )

minus2ll <- mxAlgebra( expression=MZ.objective + DZ.objective, name="m2LL" )
obj <- mxAlgebraObjective( "m2LL" )
Conf <- mxCI (c ('h2[1,1]', 'h2[2,2]', 'c2[1,1]', 'c2[2,2]', 'e2[1,1]', 'e2[2,2]','stda','stdc','stde'))
Conf2 <- mxCI (c ('Rph[2,1]', 'Ra[2,1]', 'Rc[2,1]', 'Re[2,1]') )
AceModel <- mxModel( "ACE", pars2, modelMZ, modelDZ, minus2ll, obj, Conf, Conf2)

# ------------------------------------------------------------------------------
# 4) RUN AceModel

AceFit <- mxRun(AceModel, intervals=TRUE)
AceSumm <- summary(AceFit)
AceSumm

round(AceFit@output$estimate,4)

(h2 <-mxEval(MZ.A, AceFit))
(c2 <-mxEval(MZ.C, AceFit))
(e2 <-mxEval(MZ.E, AceFit))

AceFit$MZ.V

mxCompare(SatFit,AceFit)%>

And I have a problem in the ACE-ACE result, I have no idea why there were ACE.stda[1,2], ACE.stdc[1,2], ACE.stde[1,2] in confidence intervals.
Hope for your help, THANKS!

Summary of ACE

free parameters:
name matrix row col Estimate lbound ubound
1 a_1_1 a 1 1 6.445126e-01
2 a_2_1 a 2 1 4.813791e-01
3 a_2_2 a 2 2 2.526926e-01
4 c_1_1 c 1 1 6.662922e-01
5 c_2_1 c 2 1 -7.729928e-01
6 c_2_2 c 2 2 2.637099e-05
7 e_1_1 e 1 1 3.750442e-01
8 e_2_1 e 2 1 8.722846e-02
9 e_2_2 e 2 2 3.151094e-01
10 T_11 MZ.Th 1 1 9.148539e-01 -3 3

confidence intervals:
lbound estimate ubound note
ACE.h2[1,1] 7.745476e-22 4.153965e-01 0.9240819
ACE.h2[2,2] 8.359753e-23 2.955794e-01 0.8889913
ACE.c2[1,1] 7.761656e-42 4.439454e-01 0.9134585
ACE.c2[2,2] 1.753807e-02 5.975179e-01 0.9138442
ACE.e2[1,1] 4.492750e-02 1.406581e-01 0.3351444
ACE.e2[2,2] 4.188246e-02 1.069028e-01 0.2239186
ACE.stda[1,1] -9.612906e-01 6.445126e-01 0.9612907
ACE.stda[2,1] -9.231730e-01 4.813791e-01 0.9231730
ACE.stda[1,2] NA 0.000000e+00 NA !!!
ACE.stda[2,2] -9.416668e-01 2.526926e-01 0.9416665
ACE.stdc[1,1] -9.408472e-01 6.662922e-01 0.9557514
ACE.stdc[2,1] -9.559521e-01 -7.729928e-01 0.9559520
ACE.stdc[1,2] NA 0.000000e+00 NA !!!
ACE.stdc[2,2] -9.478683e-01 2.637099e-05 0.9478681
ACE.stde[1,1] -5.789146e-01 3.750442e-01 0.5789149
ACE.stde[2,1] -1.748776e-01 8.722846e-02 0.3333841
ACE.stde[1,2] NA 0.000000e+00 NA !!!
ACE.stde[2,2] -4.665202e-01 3.151094e-01 0.4665116
ACE.Rph[2,1] -5.036850e-01 -1.720697e-01 0.1994960
ACE.Ra[2,1] -1.000000e+00 8.854216e-01 NA !!!
ACE.Rc[2,1] NA -1.000000e+00 1.0000000 !!!
ACE.Re[2,1] -5.090314e-01 2.667864e-01 0.8844167

Replied on Tue, 10/16/2018 - 11:52
Picture of user. AdminRobK Joined: 01/24/2014

In reply to by Fengxia

And I have a problem in the ACE-ACE result, I have no idea why there were ACE.stda[1,2], ACE.stdc[1,2], ACE.stde[1,2] in confidence intervals.

You requested CIs for 'stda', 'stdc', and 'stde' in line 290:

Conf <- mxCI (c ('h2[1,1]', 'h2[2,2]', 'c2[1,1]', 'c2[2,2]', 'e2[1,1]', 'e2[2,2]','stda','stdc','stde'))

Since matrices 'a', 'c', and 'e' are lower-triangular, their [1,2] elements are fixed to zero; that is also true of 'stda', 'stdc', and 'stde'.

To fit the ACE-ADE model, the simplest change to make is probably to change this,

# Matrices declared to store a, c, and e Path Coefficients
pathA <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=c(.8,-.3,.6), label=laLower("a",nv), name="a" )
pathC <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=c(.3,-.1,.3), label=laLower("c",nv), name="c" )
pathE <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=c(.8,-.3,.6), label=laLower("e",nv), name="e" )

# Matrices generated to hold A, C, and E components and total Variance
covA <- mxAlgebra( expression=a %*% t(a), name="A" )
covC <- mxAlgebra( expression=c %*% t(c), name="C" )
covE <- mxAlgebra( expression=e %*% t(e), name="E" )
covP <- mxAlgebra( expression=A+C+E, name="V" )

, to this,

# Matrices declared to store a, c, d, and e Path Coefficients
pathA <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=c(.8,-.3,.6), label=laLower("a",nv), name="a" )
pathC <- mxMatrix( type="Full", nrow=nv, ncol=nv, free=c(T,F,F,F), values=c(.3,0,0,0), labels=c("c11",NA,NA,NA), name="c" )
pathD <- mxMatrix( type="Full", nrow=nv, ncol=nv, free=c(F,F,F,T), values=c(0,0,0,0.3), labels=c(NA,NA,NA,"d22"), name="d" )
pathE <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=c(.8,-.3,.6), label=laLower("e",nv), name="e" )

# Matrices generated to hold A, C, D, and E components and total Variance
covA <- mxAlgebra( expression=a %*% t(a), name="A" )
covC <- mxAlgebra( expression=c %*% t(c), name="C" )
covD <- mxAlgebra( expression=d %*% t(d), name="D" )
covE <- mxAlgebra( expression=e %*% t(e), name="E" )
covP <- mxAlgebra( expression=A+C+D+E, name="V" )

. Then, redefine the model-expected covariance matrices to be:

covMZ <- mxAlgebra(
expression= rbind(
cbind(A+C+D+E , A+C+D),
cbind(A+C+D , A+C+D+E)), name="expCovMZ" )
covDZ <- mxAlgebra(
expression= rbind(
cbind(A+C+D+E , (0.5%x%A)+C+(0.25%x%D)),
cbind((0.5%x%A)+C+(0.25%x%D) , A+C+D+E)), name="expCovDZ" )

Feel free to add additional algebras involving 'd' and 'D' to your script (e.g., a dominance-genetic correlation matrix, or standardized 'd' paths).

Finally, I note the comment in your script about CI difficulties using CSOLNP, and that you place the identifying equality MxConstraint, varL, into both group models as well as the container model. With CSOLNP, you were almost certainly running into a known issue with redundant equality constraints.