# # Copyright 2007-2015 The OpenMx Project # # Licensed under the Apache License, Version 2.0 (the "License"); # you may not use this file except in compliance with the License. # You may obtain a copy of the License at # # http://www.apache.org/licenses/LICENSE-2.0 # # Unless required by applicable law or agreed to in writing, software # distributed under the License is distributed on an "AS IS" BASIS, # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. # See the License for the specific language governing permissions and # limitations under the License. # ----------------------------------------------------------------------------- # "Bivariate" biometric moderation example # Adapted from code provided by Conor V. Dolan # Model C is the fixed regressor model, Model B is the extended fixed regressor model, and Model A is # the full "bivariate" moderation model, a la van der Sluis, Posthuma, & Dolan (2012) # Simulate according to full Purcell ACE model # # # library(MASS) # mvrnorm library(OpenMx) # # data simulation # power good Nmz=Ndz=5000 # sample sizes in MZ and DZ number of pairs # # get the moderator data # components of moderator M am=sqrt(.4) # the variance components of the moderator Model A cm=sqrt(.2) # em=sqrt(.4) # # trait T at0=sqrt(.15) # A variance components of the trait model A (see above) at1=.025 # A moderation ct0=sqrt(.10) # C ct1=.025 # C et0=sqrt(.25) # E et1=.025 # E # # cross paths M->T in cholesky atm0=sqrt(.1) # Acovariances in tems of cholesky paths atm1=.025 # A paths are moderated ctm0=sqrt(.1) # C ctm1=.025 # C mod etm0=sqrt(.3) # E etm1=.025 # E mod # # we simulate data by first generating scores on the latent variables A,C, and E of the moderator # and A, C, and E of the trait conditional on the moderator. # These two sets are uncorrelated as the latter is trait | moderators # This is consistent with the cholesky decomp as depicted in model A # define the expected ccorelation matrices for MZ and DZ Smztmp=Sdztmp=Zerotmp=matrix(0,6,6) diag(Smztmp)=diag(Sdztmp)=1 Smztmp[4,1]=Smztmp[1,4]=1 # A Sdztmp[4,1]=Sdztmp[1,4]=.5 # A Smztmp[5,2]=Smztmp[2,5]=Sdztmp[5,2]=Sdztmp[2,5]=1 # C Smztmp2=rbind(cbind(Smztmp,Zerotmp),cbind(Zerotmp,Smztmp)) Sdztmp2=rbind(cbind(Sdztmp,Zerotmp),cbind(Zerotmp,Sdztmp)) # Latent scores: A C E (m) A C E (t|m) tmpmz=mvrnorm(Nmz,mu=rep(0,12),Sigma=Smztmp2,emp=F) tmpdz=mvrnorm(Ndz,mu=rep(0,12),Sigma=Sdztmp2,emp=F) # trait | M – just a reorganization ttmpmz=tmpmz[,7:12] ttmpdz=tmpdz[,7:12] # M– just a reorganization mtmpmz=tmpmz[,1:6] mtmpdz=tmpdz[,1:6] # – data matrices really for context mdatmz=matrix(0,Nmz,2) mdatdz=matrix(0,Ndz,2) tdatmz=matrix(0,Nmz,2) tdatdz=matrix(0,Ndz,2) # for (i in 1:Nmz) { j=1 # J=1 twin 1 mz. mod=am*mtmpmz[i,1]+cm*mtmpmz[i,2]+em*mtmpmz[i,3] # created the phenotypic moderation score # create the phenotypic trait score, depending on M and on T|M # T|M M atmp1=(at0+at1*mod)*ttmpmz[i,1] + (atm0+atm1*mod)*mtmpmz[i,1] ctmp1=(ct0+ct1*mod)*ttmpmz[i,2] + (ctm0+ctm1*mod)*mtmpmz[i,2] etmp1=(et0+et1*mod)*ttmpmz[i,3] + (etm0+etm1*mod)*mtmpmz[i,3] mdatmz[i,j]=mod # moderator tdatmz[i,j]=atmp1+ctmp1+etmp1 # trait j=2 # twin 2 mod=am*mtmpmz[i,4]+cm*mtmpmz[i,5]+em*mtmpmz[i,6] atmp1=(at0+at1*mod)*ttmpmz[i,4] + (atm0+atm1*mod)*mtmpmz[i,4] ctmp1=(ct0+ct1*mod)*ttmpmz[i,5] + (ctm0+ctm1*mod)*mtmpmz[i,5] etmp1=(et0+et1*mod)*ttmpmz[i,6] + (etm0+etm1*mod)*mtmpmz[i,6] mdatmz[i,j]=mod tdatmz[i,j]=atmp1+ctmp1+etmp1 } # same for the dz twins for (i in 1:Ndz) { j=1 mod=am*mtmpdz[i,1]+cm*mtmpdz[i,2]+em*mtmpdz[i,3] atmp1=(at0+at1*mod)*ttmpdz[i,1] + (atm0+atm1*mod)*mtmpdz[i,1] ctmp1=(ct0+ct1*mod)*ttmpdz[i,2] + (ctm0+ctm1*mod)*mtmpdz[i,2] etmp1=(et0+et1*mod)*ttmpdz[i,3] + (etm0+etm1*mod)*mtmpdz[i,3] mdatdz[i,j]=mod tdatdz[i,j]=atmp1+ctmp1+etmp1 j=2 mod=am*mtmpdz[i,4]+cm*mtmpdz[i,5]+em*mtmpdz[i,6] atmp1=(at0+at1*mod)*ttmpdz[i,4] + (atm0+atm1*mod)*mtmpdz[i,4] ctmp1=(ct0+ct1*mod)*ttmpdz[i,5] + (ctm0+ctm1*mod)*mtmpdz[i,5] etmp1=(et0+et1*mod)*ttmpdz[i,6] + (etm0+etm1*mod)*mtmpdz[i,6] mdatdz[i,j]=mod tdatdz[i,j]=atmp1+ctmp1+etmp1 } # into data frames add names and reorder the data. datmz=as.data.frame(cbind(mdatmz,mdatmz,tdatmz)) datdz=as.data.frame(cbind(mdatdz,mdatdz,tdatdz)) datmz=datmz[,c(1,2,3,5,4,6)] datdz=datdz[,c(1,2,3,5,4,6)] colnames(datmz)=colnames(datdz)=c('defm1','defm2','m1','t1','m2','t2') # --------------- Data are ready for openMx # OpenMx # ------------------------------------ Full Model Model A varnames=colnames(datmz)[3:6] # the dependent phenotypes m1 t1 m2 t2 # # big model MZ and DZ # modmodelMZ=mxModel("modmodMZ", # def …. Definition variables which we need to create the moderated paths (from M->T) mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.defm1"), name="m1" ) , mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.defm2"), name="m2" ) , # ------------------------------------------------------------ the means …. Not starting values are good. mxMatrix( type="Full", nrow=1, ncol=4, free=TRUE, values=0, labels=c("meanm","meant","meanm","meant"), name="expMean" ), # ------------------------------------------------------------------------ # ACE model moderator – cholesky elements… now defined as scalars (1x1), called lower for no good reason. # moderator parameters mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="a0m", name="A0M" ) , # path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="c0m", name="C0M" ) , # path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="e0m", name="E0M" ), # ACE model trait 0 trait parameters: main effects mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="a0t", name="A0T" ) , # path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="c0t", name="C0T" ) , # path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="e0t", name="E0T" ), # ACE model trait 1 moderator effects mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06, label="a1t", name="A1T" ) , # path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06, label="c1t", name="C1T" ) , # path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06, label="e1t", name="E1T" ), # # ACE model moderator-trait 0 moderated corvariances. First main effects mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="a0mt", name="A0MT" ) , # path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="c0mt", name="C0MT" ) , # path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="e0mt", name="E0MT" ), # ACE model moderator-trait 1 moderated covariances…. moderation mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06, label="a1mt", name="A1MT" ) , # path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06, label="c1mt", name="C1MT" ) , # path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06, label="e1mt", name="E1MT" ), # assumble parameters in to a cholesky decomp for twin 1 and twin 2 # do this in three parts : A, C, and E. # mxMatrix(type="Symm", nrow=4, ncol=4, free=F, values= c(1,0,1,0,1,0,1,1,0,1),name="PsAmz"), mxMatrix(type="Symm", nrow=4, ncol=4, free=F, values= c(1,0,1,0,1,0,1,1,0,1),name="PsC"), # # Matrices generated to hold A and E computed Variance Components # A This is a Cholesky decomposition of A for twin 1 and twin 2 (mz and dz) # not that m1 appear in the top left part (tinw 1) and m2 in the bottom right part (twin 2) # A…. mxAlgebra( rbind( cbind(A0M,0,0,0), cbind(A0MT+m1%x%A1MT,A0T+m1%x%A1T,0,0), cbind(0,0,A0M,0), cbind(0,0,A0MT+m2%x%A1MT,A0T+m2%x%A1T) ), name="chA"), # C chol mxAlgebra( rbind( cbind(C0M,0,0,0), cbind(C0MT+m1%x%C1MT,C0T+m1%x%C1T,0,0), cbind(0,0,C0M,0), cbind(0,0,C0MT+m2%x%C1MT,C0T+m2%x%C1T) ), name="chC"), # E chol mxAlgebra( rbind( cbind(E0M,0,0,0), cbind(E0MT+m1%x%E1MT,E0T+m1%x%E1T,0,0), cbind(0,0,E0M,0), cbind(0,0,E0MT+m2%x%E1MT,E0T+m2%x%E1T) ), name="chE"), # # get the expected covariance matrices mxAlgebra(chA%*%PsAmz%*%t(chA),name="Amz" ), # variance component A mxAlgebra(chC%*%PsC%*%t(chC),name="C" ), # variance component A mxAlgebra(chE%*%t(chE),name="E" ), # variance component A # # and assemble in to phenotypic 4x4 covariance matrix of MZ and DZ mxAlgebra(Amz+C+E, name="expCovMZ" ), mxData( observed=datmz, type="raw" ), mxFIMLObjective( covariance="expCovMZ", means="expMean", dimnames=varnames) ) # Same as above modmodelDZ=mxModel("modmodDZ", # def mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.defm1"), name="m1" ) , mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.defm2"), name="m2" ) , # ------------------------------------------------------------ means mxMatrix( type="Full", nrow=1, ncol=4, free=TRUE, values=0, labels=c("meanm","meant","meanm","meant"), name="expMean" ), # ------------------------------------------------------------------------ # ACE model moderator mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="a0m", name="A0M" ) , # path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="c0m", name="C0M" ) , # path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="e0m", name="E0M" ), # ACE model trait 0 mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="a0t", name="A0T" ) , # path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="c0t", name="C0T" ) , # path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="e0t", name="E0T" ), # ACE model trait 1 mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06, label="a1t", name="A1T" ) , # path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06, label="c1t", name="C1T" ) , # path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06, label="e1t", name="E1T" ), # ACE model moderator-trait 0 mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="a0mt", name="A0MT" ) , # path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="c0mt", name="C0MT" ) , # path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="e0mt", name="E0MT" ), # ACE model moderator-trait 1 mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06, label="a1mt", name="A1MT" ) , # path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06, label="c1mt", name="C1MT" ) , # path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06, label="e1mt", name="E1MT" ), # , mxMatrix(type="Symm", nrow=4, ncol=4, free=F, values= c(1,0,.5,0,1,0,.5,1,0,1),name="PsAdz"), mxMatrix(type="Symm", nrow=4, ncol=4, free=F, values= c(1,0,1,0,1,0,1,1,0,1),name="PsC"), # Matrices generated to hold A and E computed Variance Components # A mxAlgebra( rbind( cbind(A0M,0,0,0), cbind(A0MT+m1%x%A1MT,A0T+m1%x%A1T,0,0), cbind(0,0,A0M,0), cbind(0,0,A0MT+m2%x%A1MT,A0T+m2%x%A1T) ), name ="chA"), # C chol mxAlgebra( rbind( cbind(C0M,0,0,0), cbind(C0MT+m1%x%C1MT,C0T+m1%x%C1T,0,0), cbind(0,0,C0M,0), cbind(0,0,C0MT+m2%x%C1MT,C0T+m2%x%C1T) ), name ="chC"), # E chol mxAlgebra( rbind( cbind(E0M,0,0,0), cbind(E0MT+m1%x%E1MT,E0T+m1%x%E1T,0,0), cbind(0,0,E0M,0), cbind(0,0,E0MT+m2%x%E1MT,E0T+m2%x%E1T) ), name ="chE"), # # mxAlgebra(chA%*%PsAdz%*%t(chA),name="Adz" ), # variance component A mxAlgebra(chC%*%PsC%*%t(chC),name="C" ), # variance component A mxAlgebra(chE%*%t(chE),name="E" ), # variance component A # # Algebra for expected Variance/Covariance Matrices in MZ & DZ twins # this should look familiar! expected MZ covariance matrix as function of A,E mxAlgebra(Adz+C+E, name="expCovDZ" ), mxData( observed=datdz, type="raw" ), mxFIMLObjective( covariance="expCovDZ", means="expMean", dimnames=varnames) ) modmodelA=mxModel("modmodel",modmodelMZ, modmodelDZ, mxAlgebra( modmodMZ.objective + modmodDZ.objective, name="m2LL" ), mxAlgebraObjective( "m2LL" )) modmodelAfit=mxRun(modmodelA) # # a check consisting of fixing the parameters to their true values and refitting the model. # the likelihood ratio tests should not be significant (barring type I errors). # check check check ...................... check # check fix to true values ................... # 3 df modmodelAf1=omxSetParameters(modmodelA, label=c("a0m","c0m","e0m"), free=F, values=c(sqrt(.4),sqrt(.2),sqrt(.4))) modmodelAfitf1=mxRun(modmodelAf1) mxCompare(modmodelAfit,modmodelAfitf1) # 6 df modmodelAf2=omxSetParameters(modmodelA, label=c("a0mt","c0mt","e0mt","a1mt","c1mt","e1mt"), free=F, values=c(sqrt(.1),sqrt(.1),sqrt(.3),.025,.025,.025)) modmodelAfitf2=mxRun(modmodelAf2) mxCompare(modmodelAfit,modmodelAfitf2) # 6 df modmodelAf3=omxSetParameters(modmodelA, label=c("a0t","c0t","e0t","a1t","c1t","e1t"), free=F, values=c(sqrt(.15),sqrt(.1),sqrt(.25),.025,.025,.025)) modmodelAfitf3=mxRun(modmodelAf3) mxCompare(modmodelAfit,modmodelAfitf3) # check all 15df modmodelAfall=modmodelAf1 # a0m c0m e0m fixed modmodelAfall=omxSetParameters(modmodelAfall, label=c("a0t","c0t","e0t","a1t","c1t","e1t"), free=F, values=c(sqrt(.15),sqrt(.1),sqrt(.25),.025,.025,.025)) modmodelAfall=omxSetParameters(modmodelAfall, label=c("a0mt","c0mt","e0mt","a1mt","c1mt","e1mt"), free=F, values=c(sqrt(.10),sqrt(.1),sqrt(.3),.025,.025,.025)) modmodelAfitfall=mxRun(modmodelAfall) mxCompare(modmodelAfit,modmodelAfitfall) # # end check................................................................. # # now the alternative model which treats m as a fixed regressor. # there are two versions. both include moderation of A, C, and E effects on the trait # in the first version the moderator is used as a covariate as follows (model C) # y_i1 = b0 + b1*modi1 + residual_i1 # y_i2 = b0 + b1*modi2 + residual_i2 # in the second we have the extended model (model B): # y_i1 = b0 + b1*modi1 + b2*modi2 + residual_i1 # y_i2 = b0 + b1*modi2 + b2*modi1 + residual_i2 # # these can be used in some, but not all circumstances (see the moderation talk, boulder 2014). # # -------------------------------------------------------- Extended model MODEL B varnames2=c('t1','t2') datmz2=datmz[,c(1,2,4,6)] # defm1 defm2 t1 t2 datdz2=datdz[,c(1,2,4,6)] # defm1 defm2 t1 t2 # # modmodelMZ2=mxModel("modmodMZ2", # def mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.defm1"), name="m1" ) , mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.defm2"), name="m2" ) , # ------------------------------------------------------------ means mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=0, labels=c("b0"), name="B0" ), mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=0, labels=c("b1"), name="B1" ), mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=0, labels=c("b2"), name="B2" ), mxAlgebra( cbind(B0+B1*m1+B2*m2, B0+B1*m2+B2*m1),name="expMean"), # ------------------------------------------------------------------------ # ACE model trait 0 mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="a0t", name="A0T" ) , # path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="c0t", name="C0T" ) , # path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="e0t", name="E0T" ), # ACE model trait 1 mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06, label="a1t", name="A1T" ) , # path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06, label="c1t", name="C1T" ) , # path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06, label="e1t", name="E1T" ), # mxMatrix( type="Symm", nrow=2, ncol=2, free=F, values=c(1,1,1),name='PsAmz'), mxMatrix( type="Symm", nrow=2, ncol=2, free=F, values=c(1,1,1),name='PsC'), mxAlgebra(rbind(cbind(A0T+A1T*m1,0),cbind(0,A0T+A1T*m2)), name="A2"), mxAlgebra(rbind(cbind(C0T+C1T*m1,0),cbind(0,C0T+C1T*m2)), name="C2"), mxAlgebra(rbind(cbind(E0T+E1T*m1,0),cbind(0,E0T+E1T*m2)), name="E2"), # mxAlgebra(A2%*%PsAmz%*%t(A2),name="Amz" ), # variance component A mxAlgebra(C2%*%PsC%*%t(C2),name="C" ), # variance component A mxAlgebra(E2%*%t(E2),name="E" ), # variance component A # mxAlgebra(Amz+C+E, name="expCovMZ" ), mxData( observed=datmz2, type="raw" ), mxFIMLObjective( covariance="expCovMZ", means="expMean", dimnames=varnames2) ) # modmodelDZ2=mxModel("modmodDZ2", # def # def mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.defm1"), name="m1" ) , mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.defm2"), name="m2" ) , # ------------------------------------------------------------ means mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=0, labels=c("b0"), name="B0" ), mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=0, labels=c("b1"), name="B1" ), mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=0, labels=c("b2"), name="B2" ), mxAlgebra(cbind(B0+B1*m1+B2*m2, B0+B1*m2+B2*m1),name="expMean"), # ------------------------------------------------------------------------ # ACE model trait 0 mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="a0t", name="A0T" ) , # path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="c0t", name="C0T" ) , # path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6, label="e0t", name="E0T" ), # ACE model trait 1 mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06, label="a1t", name="A1T" ) , # path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06, label="c1t", name="C1T" ) , # path coefficient A mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06, label="e1t", name="E1T" ), # mxMatrix( type="Symm", nrow=2, ncol=2, free=F, values=c(1,.5,1),name='PsAdz'), mxMatrix( type="Symm", nrow=2, ncol=2, free=F, values=c(1,1,1),name='PsC'), # mxAlgebra(rbind(cbind(A0T+A1T*m1,0),cbind(0,A0T+A1T*m2)),name="A2"), mxAlgebra(rbind(cbind(C0T+C1T*m1,0),cbind(0,C0T+C1T*m2)),name="C2"), mxAlgebra(rbind(cbind(E0T+E1T*m1,0),cbind(0,E0T+E1T*m2)),name="E2"), # mxAlgebra(A2%*%PsAdz%*%t(A2),name="Adz" ), # variance component A mxAlgebra(C2%*%PsC%*%t(C2),name="C" ), # variance component A mxAlgebra(E2%*%t(E2),name="E" ), # variance component A # # mxAlgebra(Adz+C+E, name="expCovMZ" ), mxData( observed=datdz2, type="raw" ), mxFIMLObjective( covariance="expCovMZ", means="expMean", dimnames=varnames2) ) modmodelB=mxModel("modmodelB",modmodelMZ2, modmodelDZ2, mxAlgebra( modmodMZ2.objective + modmodDZ2.objective, name="m2LL" ), mxAlgebraObjective( "m2LL" )) modmodelBfit=mxRun(modmodelB) # model C...... modmodelC=omxSetParameters(modmodelB,label=c("b2"),value=0,free=F) modmodelCfit=mxRun(modmodelC)