## You are here

8 posts / 0 new Offline
Joined: 11/02/2017 - 03:33

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 Offline
Joined: 01/24/2014 - 12:15
don't split up sample

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. Offline
Joined: 11/02/2017 - 03:33
definition variable?

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 ? Offline
Joined: 01/24/2014 - 12:15
dummy coding

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.

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. Offline
Joined: 11/02/2017 - 03:33

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? Offline
Joined: 01/24/2014 - 12:15
script?

You need to redefine the A, C, and E matrices, and appropriately define the D matrix. Post your current script, and I can make some suggestions. Offline
Joined: 11/02/2017 - 03:33
I just run the ACE-ACE model,

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

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 Offline
Joined: 01/24/2014 - 12:15
some hints
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.