Bivariate twin model(ACEorADE)
Posted on
Fengxia
Joined: 11/02/2017
Forums
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
don't split up sample
Log in or register to post comments
In reply to don't split up sample by AdminRobK
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 ?
Finally, I really hope to receive your suggestion about ACE-ACE or ADE-ADE or ACE-ADE?
Log in or register to post comments
In reply to definition variable? by Fengxia
dummy coding
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.
Log in or register to post comments
In reply to dummy coding by AdminRobK
need some help about ACE-ADE
<% # 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?
Log in or register to post comments
In reply to need some help about ACE-ADE by Fengxia
script?
Log in or register to post comments
In reply to script? by AdminRobK
I just run the ACE-ACE model,
<%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
Log in or register to post comments
In reply to I just run the ACE-ACE model, by Fengxia
some hints
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.Log in or register to post comments