rm(list=ls(all=TRUE)) library(OpenMx) library(umx) library(MASS) # # ntime = t = 4 # number of time points ny = ntime # number of phenotype ny2 = ny*2 # 2x number of phenotypes ny for each twin # # DATA #--------------------------------------------------------------- # iqdat <- read.table(file="IQdat2018",header=TRUE) # head(iqdat) dim(iqdat) # nmz <- sum(iqdat$zyg==1) # number of MZ twins in data ndz <- sum(iqdat$zyg==2) # number of MZ twins in data ntot <- ndz + nmz # total iqdatmz=iqdat[iqdat$zyg==1,2:9] iqdatdz=iqdat[iqdat$zyg==2,2:9] varnames=c( c('iq5_T1','iq7_T1','iq10_T1','iq12_T1'), c('iq5_T2','iq7_T2','iq10_T2','iq12_T2')) colnames(iqdatmz)=colnames(iqdatdz)=varnames # check print(nmz) print(dim(iqdatmz)) print(dim(iqdatdz)) # # Get starting values for SA SC SE unconstrained model # Smz=cov(iqdatmz,use="pairwise.complete.obs") # MZ covariance matrix Smz1=(Smz[1:4,1:4]/3 + Smz[5:8,5:8]/3)/2 # divide by three (A=C=E effects) # # starting values # stSA=Smz1 stSC=Smz1 stSE=Smz1 meanMZ=apply(iqdatmz,2,mean,na.rm=T) # means meanDZ=apply(iqdatdz,2,mean,na.rm=T) # means stmean=(meanMZ[1:t]+meanDZ[1:t])/2 # starting values for the means # # model: saturated model selvars=varnames doSimplex=T doACESAT=T doSAT=T # # mxOption(NULL,"Default optimizer","NPSOL") if (doSAT) { # Satmz=mxModel('Satmz', mxMatrix(type='Stand',nrow=ny2,ncol=ny2,free=TRUE,value=.5, lbound=-.9,ubound=.9,name='cormz'), mxMatrix(type='Diag',nrow=ny2,ncol=ny2,free=TRUE,value=15,name='sdmz'), # IQ data variance of about 225 stdev of 15 mxMatrix(type='Full',nrow=1,ncol=ny2,free=TRUE, labels=c('m1','m2','m3','m4','m1','m2','m3','m4'), value=stmean,name='memz'), mxAlgebra(expression=sdmz%*%cormz%*%t(sdmz),name='Ssatmz') ) Satdz=mxModel('Satdz', mxMatrix(type='Stand',nrow=ny2,ncol=ny2,free=TRUE,value=.4, lbound=-.9,ubound=.9,name='cordz'), mxMatrix(type='Diag',nrow=ny2,ncol=ny2,free=TRUE,value=15,name='sddz'), mxMatrix(type='Full',nrow=1,ncol=ny2,free=TRUE, labels=c('m1','m2','m3','m4','m1','m2','m3','m4'), value=stmean,name='medz'), mxAlgebra(expression=sddz%*%cordz%*%t(sddz),name='Ssatdz') ) # # a model the data, the fit function (MZ) MZmodel0=mxModel("MZ", mxData( observed=iqdatmz, type="raw"), # the data mxExpectationNormal( covariance="Satmz.Ssatmz", means="Satmz.memz", dimnames=selvars), # the fit function mxFitFunctionML() ) # a model the data, the fit function (DZ) DZmodel0=mxModel("DZ", mxData( observed=iqdatdz, type="raw"), mxExpectationNormal( covariance="Satdz.Ssatdz", means="Satdz.medz", dimnames=selvars), mxFitFunctionML() ) twinModel0 <- mxModel("twinmodel0", Satmz, Satdz, MZmodel0, DZmodel0, mxAlgebra(MZ.objective + DZ.objective, name="minus2loglikelihood"), mxFitFunctionAlgebra("minus2loglikelihood") ) # fit the model twinModel0_out <- mxRun(twinModel0) # Rmz_est=round(mxEval(Satmz.cormz,twinModel0_out),3) Rdz_est=round(mxEval(Satdz.cordz,twinModel0_out),3) Means_est=round(mxEval(Satdz.medz,twinModel0_out),3) # } # end doSAT # # selDVs=colnames(iqdatmz)[1:8] dzData=iqdatdz mzData=iqdatmz twinModel1unx_out=umxACEv(name = "ACE", selDVs, selCovs=NULL, dzData, mzData, covMethod=NULL) # if (doACESAT) { # ACE simplex model SATACEs=mxModel("SATace", mxMatrix(type='Symm',nrow=ny,ncol=ny,free=TRUE, value=stSA,name='SA'), mxMatrix(type='Symm',nrow=ny,ncol=ny,free=TRUE, value=stSC,name='SC'), mxMatrix(type='Symm',nrow=ny,ncol=ny,free=TRUE, value=stSE,name='SE') ) # SumstatMZ=mxModel("sumstatmz", mxAlgebra(expression=SATace.SA+SATace.SC+SATace.SE, name='SigmaPh11'), mxAlgebra(expression=SATace.SA+SATace.SC, name='SigmaPh21'), mxAlgebra(expression=cbind(rbind(SigmaPh11,SigmaPh21),rbind(SigmaPh21,SigmaPh11)),name='SigmaMZ'), mxMatrix(type='Full',nrow=1,ncol=ny2,free=TRUE, labels=c('m1','m2','m3','m4','m1','m2','m3','m4'), value=stmean,name='MuMZ')) # SumstatDZ=mxModel("sumstatdz", mxAlgebra(expression=SATace.SA+SATace.SC+SATace.SE, name='SigmaPh11'), mxAlgebra(expression=.5%x%SATace.SA+SATace.SC, name='SigmaPh21'), mxAlgebra(expression=cbind(rbind(SigmaPh11,SigmaPh21),rbind(SigmaPh21,SigmaPh11)),name='SigmaDZ'), mxMatrix(type='Full',nrow=1,ncol=ny2,free=TRUE, labels=c('m1','m2','m3','m4','m1','m2','m3','m4'), value=stmean,name='MuDZ')) # a model the data, the fit function (MZ) MZmodel1=mxModel("MZ", mxData( observed=iqdatmz, type="raw"), # the data mxExpectationNormal( covariance="sumstatmz.SigmaMZ", means="sumstatmz.MuMZ", dimnames=selvars), # the fit function mxFitFunctionML() ) # a model the data, the fit function (DZ) DZmodel1=mxModel("DZ", mxData( observed=iqdatdz, type="raw"), mxExpectationNormal(covariance="sumstatdz.SigmaDZ", means="sumstatdz.MuDZ", dimnames=selvars), mxFitFunctionML() ) # twinModel1 <- mxModel("twinmodel1", SATACEs, SumstatMZ,SumstatDZ, MZmodel1, DZmodel1, mxAlgebra(MZ.objective + DZ.objective, name="minus2loglikelihood"), mxFitFunctionAlgebra("minus2loglikelihood") ) twinModel1_out <- mxRun(twinModel1) # mxCompare(twinModel0_out,twinModel1_out) # fit of ACE cholesky SA_est=round(mxEval(SATace.SA,twinModel1_out),3) # unconstrained SA and RA RA_est=round(cov2cor(SA_est),3) SC_est=round(mxEval(SATace.SC,twinModel1_out),3) # unconstrained SC and RC RC_est=round(cov2cor(SC_est),3) SE_est=round(mxEval(SATace.SE,twinModel1_out),3) # unconstrained SE and RE RE_est=round(cov2cor(SE_est),3) } # end do cholesky if (doSimplex) { # # ACE simplex model # SAmodel=mxModel("SA", mxMatrix(type='Diag',nrow=ny,ncol=ny,free=c(T,T,T,T), labels=c("PsA11","PsA22","PsA33","PsA44"), value=c(100,5,5,5),name='PsA'), mxMatrix(type='Diag',nrow=ny,ncol=ny,free=c(T,T,T,T), labels=c("TeA11","TeA11","TeA11","TeA11"), value=c(2,2,2,2),name='TeA'), mxMatrix(type='Iden',nrow=ny,ncol=ny,name='I'), mxMatrix(type='Full',nrow=ny,ncol=ny, free=matrix(c(F,F,F,F, T,F,F,F, F,T,F,F, F,F,T,F),ny,ny,byrow=T), labels=matrix(c('NA','NA','NA','NA', 'BeA21','NA','NA','NA', 'NA','BeA32','NA','NA', 'NA','NA','BeA43','NA'),ny,ny,byrow=T), value=matrix(c(0,0,0,0, .9,0,0,0, 0,.9,0,0, 0,0,.9,0),ny,ny,byrow=T), name='BeA'), mxAlgebra(expression=solve(I-BeA),name='iBeA'), mxAlgebra(expression=iBeA%*%(PsA)%*%t(iBeA)+(TeA),name='SigmaA') ) SCmodel=mxModel("SC", mxMatrix(type='Diag',nrow=ny,ncol=ny,free=c(T,T,T,T), labels=c("PsC11","PsC22","PsC33","PsC44"), value=c(70,10,10,10),name='PsC'), mxMatrix(type='Diag',nrow=ny,ncol=ny,free=c(T,T,T,T), labels=c("TeC11","TeC11","TeC11","TeC11"), value=c(5,5,5,5),name='TeC'), mxMatrix(type='Iden',nrow=ny,ncol=ny,name='I'), mxMatrix(type='Full',nrow=ny,ncol=ny, free=matrix(c(F,F,F,F, T,F,F,F, F,T,F,F, F,F,T,F),ny,ny,byrow=T), labels=matrix(c('NA','NA','NA','NA', 'BeC21','NA','NA','NA', 'NA','BeC32','NA','NA', 'NA','NA','BeC43','NA'),ny,ny,byrow=T), value=matrix(c(0,0,0,0, .8,0,0,0, 0,.8,0,0, 0,0,.8,0),ny,ny,byrow=T), name='BeC'), mxAlgebra(expression=solve(I-BeC),name='iBeC'), mxAlgebra(expression=iBeC%*%(PsC)%*%t(iBeC)+(TeC),name='SigmaC') ) SEmodel=mxModel("SE", mxMatrix(type='Diag',nrow=ny,ncol=ny,free=c(F,F,F,F), labels=c("PsE11","PsE22","PsE33","PsE44"), value=c(0,0,0,0),name='PsE'), mxMatrix(type='Diag',nrow=ny,ncol=ny,free=c(T,T,T,T), labels=c("TeE11","TeE22","TeE33","TeE44"), value=c(50,50,50,50),name='TeE'), mxMatrix(type='Iden',nrow=ny,ncol=ny,name='I'), mxMatrix(type='Full',nrow=ny,ncol=ny, free=matrix(c(F,F,F,F, F,F,F,F, F,F,F,F, F,F,F,F),ny,ny,byrow=T), labels=matrix(c('NA','NA','NA','NA', 'BeE21','NA','NA','NA', 'NA','BeE32','NA','NA', 'NA','NA','BeE43','NA'),ny,ny,byrow=T), value=matrix(c(0,0,0,0, .0,0,0,0, 0,.0,0,0, 0,0,.0,0),ny,ny,byrow=T), name='BeE'), mxAlgebra(expression=solve(I-BeE),name='iBeE'), mxAlgebra(expression=iBeE%*%(PsE)%*%t(iBeE)+(TeE),name='SigmaE') ) SumstatMZ=mxModel("sumstatmz", mxAlgebra(expression=SA.SigmaA+SC.SigmaC+SE.SigmaE, name='SigmaPh11'), mxAlgebra(expression=SA.SigmaA+SC.SigmaC, name='SigmaPh21'), mxAlgebra(expression=cbind(rbind(SigmaPh11,SigmaPh21),rbind(SigmaPh21,SigmaPh11)),name='SigmaMZ'), mxMatrix(type='Full',nrow=1,ncol=ny2,free=TRUE, labels=c('m1','m2','m3','m4','m1','m2','m3','m4'), value=stmean,name='MuMZ')) SumstatDZ=mxModel("sumstatdz", mxAlgebra(expression=SA.SigmaA+SC.SigmaC+SE.SigmaE, name='SigmaPh11'), mxAlgebra(expression=.5%x%(SA.SigmaA)+SC.SigmaC, name='SigmaPh21'), mxAlgebra(expression=cbind(rbind(SigmaPh11,SigmaPh21),rbind(SigmaPh21,SigmaPh11)),name='SigmaDZ'), mxMatrix(type='Full',nrow=1,ncol=ny2,free=TRUE, labels=c('m1','m2','m3','m4','m1','m2','m3','m4'), value=stmean,name='MuDZ')) # a model the data, the fit function (MZ) MZmodel1=mxModel("MZ", mxData( observed=iqdatmz, type="raw"), # the data mxExpectationNormal( covariance="sumstatmz.SigmaMZ", means="sumstatmz.MuMZ", dimnames=selvars), # the fit function mxFitFunctionML() ) # a model the data, the fit function (DZ) DZmodel1=mxModel("DZ", mxData( observed=iqdatdz, type="raw"), mxExpectationNormal(covariance="sumstatdz.SigmaDZ", means="sumstatdz.MuDZ", dimnames=selvars), mxFitFunctionML() ) twinModel2 <- mxModel("twinmodel2", SAmodel,SCmodel,SEmodel, SumstatMZ,SumstatDZ, MZmodel1, DZmodel1, mxAlgebra(MZ.objective + DZ.objective, name="minus2loglikelihood"), mxFitFunctionAlgebra("minus2loglikelihood") ) twinModel2_out <- mxRun(twinModel2) } # end do simplex mxCompare(twinModel1_out,twinModel2_out) # fit of ACE cholesky 3df test # #PsE_est=round(mxEval(SE.PsE,twinModel2_out)^2,3) #BeE_est=round(mxEval(SE.BeE,twinModel2_out),3) TeE_est=round(mxEval(SE.TeE,twinModel2_out),3) # PsC_est=round(mxEval(SC.PsC,twinModel2_out),3) BeC_est=round(mxEval(SC.BeC,twinModel2_out),3) TeC_est=round(mxEval(SC.TeC,twinModel2_out),3) # PsA_est=round(mxEval(SA.PsA,twinModel2_out),3) BeA_est=round(mxEval(SA.BeA,twinModel2_out),3) TeA_est=round(mxEval(SA.TeA,twinModel2_out),3) # # twinModel3 = omxSetParameters(twinModel2,labels=c("TeA11"), values=c(0),free=c(F)) # twinModel3_out <- mxRun(twinModel3) mxCompare(twinModel2_out,twinModel3_out) # TeE_est=round(mxEval(SE.TeE,twinModel3_out),3) # PsC_est=round(mxEval(SC.PsC,twinModel3_out),3) BeC_est=round(mxEval(SC.BeC,twinModel3_out),3) TeC_est=round(mxEval(SC.TeC,twinModel3_out),3) # PsA_est=round(mxEval(SA.PsA,twinModel3_out),3) BeA_est=round(mxEval(SA.BeA,twinModel3_out),3) # # twinModel4 = omxSetParameters(twinModel3,labels=c("PsA22","PsA33","PsA44"), values=c(0),free=c(F)) twinModel4_out <- mxRun(twinModel4) mxCompare(twinModel3_out,twinModel4_out) # TeE_est=round(mxEval(SE.TeE,twinModel4_out),3) # PsC_est=round(mxEval(SC.PsC,twinModel4_out),3) BeC_est=round(mxEval(SC.BeC,twinModel4_out),3) TeC_est=round(mxEval(SC.TeC,twinModel4_out),3) # PsA_est=round(mxEval(SA.PsA,twinModel4_out),3) BeA_est=round(mxEval(SA.BeA,twinModel4_out),3) #TeA_est=round(mxEval(SA.TeA,twinModel4_out),3) # SC_est=round(mxEval(SC.SigmaC,twinModel4_out),3) SA_est=round(mxEval(SA.SigmaA,twinModel4_out),3) SE_est=round(mxEval(SE.SigmaE,twinModel4_out),3) Sph_est=SC_est+SA_est+SE_est h2=round(diag(SA_est) / diag(Sph_est),3) c2=round(diag(SC_est) / diag(Sph_est),3) e2=round(diag(SE_est) / diag(Sph_est),3) Rph_est=cov2cor(Sph_est) # diag(SA_est) diag(SC_est) diag(SE_est) round(SA_est/Sph_est,3) round(SC_est/Sph_est,3) round(SE_est/Sph_est,3)