#------------------------------------------------------------------------------------- # Script: BivOrdCont.R # Bivariate Twin model for Combined Continuous and Ordered Categorical variables # Data File: ADHDiqage.csv # Phenotypes # Variable 1: IQ continuous variable # Variable 2: ADHD categorised in 3 classes: 0 (unafected), # 1 (above TH for a clinical diagnosis), 2 (sub-clinical group) # Type of data: Raw ordinal data # ZYG: 1=MZ 2=DZ # ------------------------------------------------------------------------------------- # Note: Continuous variable(s) are always selected first, then Ordinal variable(s) rm(list=ls()) #ls() require(OpenMx) require(psych) Bivdata <- read.csv ('bivtwindata.csv', header=T, sep=",", na.strings="NA") names(Bivdata) summary(Bivdata) str(Bivdata) describe(Bivdata) BivdataOrd <- data.frame(Bivdata) #Make the integer variables ordered factors (Note:because of the 'cuting' above, this is not necessary here) BivdataOrd$sleepda1 <-mxFactor(BivdataOrd$sleepda1, levels=c(0:2) ) BivdataOrd$sleepda2 <-mxFactor(BivdataOrd$sleepda2, levels=c(0:2) ) describe(BivdataOrd) vars <-c('Hb','sleepda') selVars <-c('Hb1', 'sleepda1', 'Hb2', 'sleepda2' ) useVars <-c('Hb1', 'sleepda1', 'Hb2', 'sleepda2', 'age1', 'age2') # Select Data for Analysis mzData <- subset(BivdataOrd, zygosity==1, useVars) dzData <- subset(BivdataOrd, zygosity==2, useVars) describe(mzData) describe(dzData) # To get the CTs for ADHD table(mzData$sleepda1, mzData$sleepda2) table(dzData$sleepda1, dzData$sleepda2) nv <- 2 # number of variables per twin nvo <- 1 # number of ordinal variables per twin nvc <- nv-nvo # number of continuous variables per twin poso <- nvc+1 # position where ordinal variables start ntv <- nv*2 # number of variables per pair nth <- 2 # number of max thresholds ninc <- nth-1 # number of max increments ## ----------------------------------------------------------------------------------------- # Part 1: # Constrained Polyserial correlation model # MZ and DZ Thresholds, but constrained across twins # Age effect is different acros variables, but same across thresholds within ordinal variable(s)if c>2 # There is one overall rPH between var1-2 and the x-trait x-twin correlations are symmetric ## ------------------------------------------------------------------------------------------------------------------------------ # CREATE LABELS in objects(to ease specification) LabThMZ <-c('Tmz_11','imz_11') # THs for ordinal var for a twin individual (mz) LabThDZ <-c('Tdz_11','idz_11') # THs for ordinal var for a twin individual (dz) LabCorMZ <-c('r21','rMZ1','MZxtxt','MZxtxt','rMZ2','r21') LabCorDZ <-c('r21','rDZ1','DZxtxt','DZxtxt','rDZ2','r21') (LabCovM <-c( paste("Bcv",1:nvc,sep=""), rep(NA, nvo))) (LabCovTH <-rep( paste("Bov",1:nvo,sep=""), nth)) (LabM <- c( paste("m",1:nv,sep=""))) (LabSD <- c( paste("sd",1:nv,sep=""))) (PatSD <- c( rep(T,nvc), rep(F, nvo))) (PatM <- c( rep(T,nvc), rep(F, nvo))) (PatBetaM <- c( rep(T,nvc), rep(F, nvo))) # CREATE START VALUES in objects #(StM <-c( colMeans(mzData[,1:nvc],na.rm=TRUE), rep(0,nvo))) StM <-c( 100, 0) (StSD <-c( sd(mzData[,1:nvc],na.rm=TRUE), rep(1,nvo)) ) (StBetaM <-c( rep(.2,nvc), rep(0,nvo))) StCorMZ <-c(-.3, .7, -.3,-.3, .7, -.3) StCorDZ <-c(-.3, .4, -.15,-.15, .3, -.3) StTH <-c(-1.4, 1.7 ) # Define definition variables to hold the Covariates obsAge1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=F, labels=c("data.age1"), name="Age1") obsAge2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=F, labels=c("data.age2"), name="Age2") # Matrix & Algebra for expected observed means, Thresholds, effect of Age on Th and correlations Mean <-mxMatrix( type="Full", nrow=1, ncol=nv, free=PatM, values=StM, labels=LabM, name="M" ) betaAm <-mxMatrix( type="Full", nrow=nv, ncol=1, free=PatBetaM, values=StBetaM, labels=LabCovM, name="BageM" ) betaAth <-mxMatrix( type="Full", nrow=nth, ncol=1, free=T, values=.2, labels=LabCovTH, name="BageTH" ) Mu1 <-mxAlgebra( expression= M + t(BageM%*%Age1), name="MeanCtw1" ) Mu2 <-mxAlgebra( expression= M + t(BageM%*%Age2), name="MeanCtw2" ) Mu12 <-mxAlgebra( expression= cbind(MeanCtw1,MeanCtw2), name="ExpMean" ) inc <-mxMatrix( type="Lower",nrow=nth, ncol=nth, free=F, values=1, name="Low") Tmz <-mxMatrix( type="Full", nrow=nth, ncol=nvo, free=T, values=StTH, lbound=c(-3,rep(.001,ninc)), ubound=3, labels=LabThMZ, name="ThMZ") ThresMZ <-mxAlgebra( expression= cbind(Low%*%ThMZ + BageTH%x%Age1, Low%*%ThMZ + BageTH%x%Age2), name="ExpThresMZ") Tdz <-mxMatrix( type="Full", nrow=nth, ncol=nvo, free=T, values=StTH, lbound=c(-3,rep(.001,ninc)), ubound=3, labels=LabThDZ, name="ThDZ") ThresDZ <-mxAlgebra( expression= cbind(Low%*%ThDZ + BageTH%x%Age1, Low%*%ThDZ + BageTH%x%Age2), name="ExpThresDZ") SD <-mxMatrix( type="Diag", nrow=ntv, ncol=ntv, free=c(PatSD,PatSD), values=c(StSD,StSD), labels=c(LabSD,LabSD), name="sd") rMZ <-mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=StCorMZ, labels=LabCorMZ, lbound=-.99, ubound=.99, name="CorMZ") rDZ <-mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=StCorDZ, labels=LabCorDZ, lbound=-.99, ubound=.99, name="CorDZ") # Matrices for the Covariance model MZCov <-mxAlgebra( expression=sd %&% CorMZ, name="ExpCovMZ" ) DZCov <-mxAlgebra( expression=sd %&% CorDZ, name="ExpCovDZ" ) # Data objects DataMZ <-mxData(observed=mzData, type="raw") DataDZ <-mxData(observed=dzData, type="raw") # Objective objects for Multiple Groups objmz <-mxFIMLObjective( covariance="ExpCovMZ", means="ExpMean", dimnames=selVars, thresholds="ExpThresMZ", threshnames=c("sleepda1","sleepda2")) objdz <-mxFIMLObjective( covariance="ExpCovDZ", means="ExpMean", dimnames=selVars, thresholds="ExpThresDZ", threshnames=c("sleepda1","sleepda2")) # Combine Groups pars <-list (obsAge1,obsAge2,Mean,betaAm,betaAth,Mu1,Mu2,Mu12,inc,SD) modelMZ <-mxModel(pars, Tmz, ThresMZ, rMZ, MZCov, DataMZ, objmz, name="MZ" ) modelDZ <-mxModel(pars, Tdz, ThresDZ, rDZ, DZCov, DataDZ, objdz, name="DZ" ) minus2ll <-mxAlgebra( expression=MZ.objective + DZ.objective, name="m2LL" ) obj <-mxAlgebraObjective("m2LL" ) Conf1 <-mxCI (c ('MZ.rMZ[2,1]' ) ) # IQ-ADHD phenotypic correlations Conf2 <-mxCI (c ('MZ.rMZ[3,1]', 'MZ.rMZ[4,2]','MZ.rMZ[3,2]') ) # MZ twin cor var1, var 2 and xtr-xtwin Conf3 <-mxCI (c ('DZ.rDZ[3,1]', 'DZ.rDZ[4,2]','DZ.rDZ[3,2]') ) # DZ twin cor var1, var 2 and xtr-xtwin CorModel <-mxModel( "Cor", modelMZ, modelDZ, minus2ll, obj, Conf1, Conf2, Conf3) # RUN Correlation MODEL CorFit <-mxRun(CorModel, intervals=F) (CorSumm <-summary(CorFit)) round(CorFit@output$estimate,4) CorFit$Cor$MZ.Rph CorFit$Cor$MZ.Rbmz CorFit$Cor$DZ.Rbdz # ******************************************************************************************************** # ------------------------------------------------------------------------------------------------------------------- # (2) Specify Bivariate ACE Model using Cholesky Decomposition, # One set of THRESHOLDS (same across twins and zyg group) # ------------------------------------------------------------------------------------------------------------------- LabTh <-c('T_11','i_11') # THs and inc for ordinal var for a twin individual # CREATE LABELS FOR Lower Triangular Matrices (aLabs <- paste("a", do.call(c, sapply(seq(1, nv), function(x){ paste(x:nv, x,sep="") })), sep="")) (cLabs <- paste("c", do.call(c, sapply(seq(1, nv), function(x){ paste(x:nv, x,sep="") })), sep="")) (eLabs <- paste("e", do.call(c, sapply(seq(1, nv), function(x){ paste(x:nv, x,sep="") })), sep="")) # CREATE START VALUES FOR Lower Triangular Matrices (SDcon <-(sqrt(vech(cov(mzData[,1:nvc],mzData[,1:nvc],use="complete")))))/3 Stpaths <- c(5,-0.2,3) # Define definition variables to hold the Covariates obsAge1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=F, labels=c("data.age1"), name="Age1") obsAge2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=F, labels=c("data.age2"), name="Age2") # Matrix & Algebra for expected observed means, Thresholds, effect of Age on Th and correlations Mean <-mxMatrix( type="Full", nrow=1, ncol=nv, free=PatM, values=StM, labels=LabM, name="M" ) betaAm <-mxMatrix( type="Full", nrow=nv, ncol=1, free=PatBetaM, values=StBetaM, labels=LabCovM, name="BageM" ) betaAth <-mxMatrix( type="Full", nrow=nth, ncol=1, free=T, values=.2, labels=LabCovTH, name="BageTH" ) Mu1 <-mxAlgebra( expression= M + t(BageM%*%Age1), name="MeanCtw1" ) Mu2 <-mxAlgebra( expression= M + t(BageM%*%Age2), name="MeanCtw2" ) Mu12 <-mxAlgebra( expression= cbind(MeanCtw1,MeanCtw2), name="ExpMean" ) inc <-mxMatrix( type="Lower",nrow=nth, ncol=nth, free=F, values=1, name="Low") T <-mxMatrix( type="Full", nrow=nth, ncol=nvo, free=T, values=StTH, lbound=c(-3,rep(.001,ninc)), ubound=3, labels=LabTh, name="Th") Thres <-mxAlgebra( expression= cbind(Low%*%Th + BageTH%x%Age1, Low%*%Th + BageTH%x%Age2), name="ExpThres") # MATRICES FOR THE ACE COV MODEL # Matrices to store a, c, and e Path Coefficients pathA <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=Stpaths, labels=aLabs, name="a" ) pathC <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=Stpaths, labels=cLabs, name="c" ) pathE <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=Stpaths, labels=eLabs, name="e" ) # Matrices generated to hold A, C, and E computed Variance Components covA <- mxAlgebra( expression=a %*% t(a), name="A" ) covC <- mxAlgebra( expression=c %*% t(c), name="C" ) covE <- mxAlgebra( expression=e %*% t(e), name="E" ) # 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") # Algebra to compute Phenotypic, A, C & E correlations matI <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I") 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" ) # 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" ) # Unity constraint on total variance of Ordinal variable(s) mUnv <- mxMatrix( type="Unit", nrow=nvo, ncol=1, name="Unv" ) varL <- mxConstraint( expression=diag2vec(V[(poso:nv),(poso:nv)])==Unv, name="VarL" ) # Data objects DataMZ <-mxData(observed=mzData, type="raw") DataDZ <-mxData(observed=dzData, type="raw") # Objective objects for Multiple Groups objmz <-mxFIMLObjective( covariance="ExpCovMZ", means="ExpMean", dimnames=selVars, thresholds="ExpThres", threshnames=c("sleepda1","sleepda2")) objdz <-mxFIMLObjective( covariance="ExpCovDZ", means="ExpMean", dimnames=selVars, thresholds="ExpThres", threshnames=c("sleepda1","sleepda2")) # Combine Groups pars1 <-list(obsAge1,obsAge2,Mean,betaAm,betaAth,Mu1,Mu2,Mu12,inc,T,Thres,mUnv,varL) pars2 <-list(pathA, pathC, pathE, covA, covC, covE, covP, StA, StC, StE, matI, rph, rA, rC, rE ) 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" ) Conf1 <- mxCI (c ('h2[1,1]', 'h2[2,2]', 'c2[1,1]', 'c2[2,2]', 'e2[1,1]', 'e2[2,2]') ) Conf2 <- mxCI (c ('Rph[2,1]', 'Ra[2,1]', 'Rc[2,1]', 'Re[2,1]') ) AceModel <- mxModel( "ACE", pars2, modelMZ, modelDZ, minus2ll, obj, Conf1, Conf2) # ------------------------------------------------------------------------------ # RUN Constrained Bivariate ACE Model AceFit <- mxRun(AceModel, intervals=TRUE) (AceSumm <- summary(AceFit)) round(AceFit@output$estimate,4) round(AceFit$Est@result,4) AceFit$ACE.h2 AceFit$ACE.c2 AceFit$ACE.e2 AceFit$ACE.Rph AceFit$ACE.Ra AceFit$ACE.Rc AceFit$ACE.Re # Generate AceModel Output # Print Comparative Fit Statistics # ----------------------------------------------------------------------- mxCompare(CorFit,AceFit) # Fit ACE_noa1 model # ----------------------------------------------------------------------- ACE_noa1_Cholesky_Model <- omxSetParameters(model=AceFit,labels=c('a11','a21'),free=F,values=0,name="ACE_noa1_Cholesky") ACE_noa1_Cholesky_Fit <- mxRun(ACE_noa1_Cholesky_Model, intervals=TRUE) ACE_noa1_Cholesky_Summ <- summary(ACE_noa1_Cholesky_Fit) ACE_noa1_Cholesky_Summ mxCompare(AceFit,ACE_noa1_Cholesky_Fit) tableFitStatistics(AceFit,ACE_noa1_Cholesky_Fit) # Fit ACE_noa2 model # ----------------------------------------------------------------------- ACE_noa2_Cholesky_Model <- omxSetParameters(model=AceFit,labels="a22",free=F,values=0,name="ACE_noa2_Cholesky") ACE_noa2_Cholesky_Fit <- mxRun(ACE_noa2_Cholesky_Model, intervals=TRUE) ACE_noa2_Cholesky_Summ <- summary(ACE_noa2_Cholesky_Fit) ACE_noa2_Cholesky_Summ mxCompare(AceFit,ACE_noa2_Cholesky_Fit) tableFitStatistics(AceFit,ACE_noa2_Cholesky_Fit) # Fit ACE_noa model # ----------------------------------------------------------------------- ACE_noa_Cholesky_Model <- omxSetParameters(model=AceFit,labels=c('a11','a21','a22'),free=F,values=0,name="ACE_noa_Cholesky") ACE_noa_Cholesky_Fit <- mxRun(ACE_noa_Cholesky_Model, intervals=TRUE) ACE_noa_Cholesky_Summ <- summary(ACE_noa_Cholesky_Fit) ACE_noa_Cholesky_Summ mxCompare(AceFit,ACE_noa_Cholesky_Fit) tableFitStatistics(AceFit,ACE_noa_Cholesky_Fit) CholeskyACENesteda <- list(ACE_noa1_Cholesky_Fit, ACE_noa1_Cholesky_Fit, ACE_noa_Cholesky_Fit) mxCompare(AceFit,CholeskyACENesteda) # Fit ACE_noc1 model # ----------------------------------------------------------------------- ACE_noc1_Cholesky_Model <- omxSetParameters(model=AceFit,labels=c('c11','c21'),free=F,values=0,name="ACE_noc1_Cholesky") ACE_noc1_Cholesky_Fit <- mxRun(ACE_noc1_Cholesky_Model, intervals=TRUE) ACE_noc1_Cholesky_Summ <- summary(ACE_noc1_Cholesky_Fit) ACE_noc1_Cholesky_Summ mxCompare(AceFit,ACE_noc1_Cholesky_Fit) tableFitStatistics(AceFit,ACE_noc1_Cholesky_Fit) # Fit ACE_noc2 model # ----------------------------------------------------------------------- ACE_noc2_Cholesky_Model <- omxSetParameters(model=AceFit,labels="c22",free=F,values=0,name="ACE_noc2_Cholesky") ACE_noc2_Cholesky_Fit <- mxRun(ACE_noc2_Cholesky_Model, intervals=TRUE) ACE_noc2_Cholesky_Summ <- summary(ACE_noc2_Cholesky_Fit) ACE_noc2_Cholesky_Summ mxCompare(AceFit,ACE_noc2_Cholesky_Fit) tableFitStatistics(AceFit,ACE_noc2_Cholesky_Fit) # Fit ACE_noc model # ----------------------------------------------------------------------- ACE_noc_Cholesky_Model <- omxSetParameters(model=AceFit,labels=c('c11','c21','c22'),free=F,values=0,name="ACE_noc_Cholesky") ACE_noc_Cholesky_Fit <- mxRun(ACE_noc_Cholesky_Model, intervals=TRUE) ACE_noc_Cholesky_Summ <- summary(ACE_noc_Cholesky_Fit) ACE_noc_Cholesky_Summ mxCompare(AceFit,ACE_noc_Cholesky_Fit) tableFitStatistics(AceFit,ACE_noc_Cholesky_Fit) CholeskyACENestedc <- list(ACE_noc1_Cholesky_Fit, ACE_noc1_Cholesky_Fit, ACE_noc_Cholesky_Fit) mxCompare(AceFit,CholeskyACENestedc) mxCompare(AceFit,CholeskyACENesteda,CholeskyACENestedc) # Fit ACE_noa21 model # ----------------------------------------------------------------------- #ACE_noa21_Cholesky_Model <- mxModel(AceFit, name="ACE_noa21_Cholesky", # mxModel(AceFit$ACE, # mxMatrix( type="Lower", nrow=nv, ncol=nv, free=F, values=0, name="a" ) # drop a21 at 0 # )) #CholAeModel <- omxSetParameters( CholAeModel, labels=c('c11','c21','c31','c22','c32','c33'), free=FALSE, values=0 ) ACE_noa21_Cholesky_Model <- omxSetParameters(model=AceFit,labels="a21",free=F,values=0,name="ACE_noa21_Cholesky") ACE_noa21_Cholesky_Fit <- mxRun(ACE_noa21_Cholesky_Model, intervals=TRUE) reACE_noa21_Cholesky_Fit <- mxRun(ACE_noa21_Cholesky_Fit) ACE_noa21_Cholesky_Summ <- summary(ACE_noa21_Cholesky_Fit) ACE_noa21_Cholesky_Summ mxCompare(AceFit,ACE_noa21_Cholesky_Fit) tableFitStatistics(AceFit,ACE_noa21_Cholesky_Fit)