#_________________________ # Following I run each ADHD scale item through a saturated, ACE, # ADE, AE and E model to look for variations in genetic and # environmental aetiology across symptoms while controlling # for gender and age effects #__________________________ require(OpenMx) require(psych) source('GenEpiHelperFunctions.R') nu <-read.csv('nu/nu.csv') names(nu) dim(nu) # Dealing with missing definition variables (no NA's allowed) table(twins$age_1) table(twins$age_2) twins$age_1[is.na(twins$age_1)]<-twins$age_2[is.na(twins$age_1)] twins$age_2[is.na(twins$age_2)]<-twins$age_1[is.na(twins$age_2)] table(twins$age_1) table(twins$age_2) twins$sex_1[is.na(twins$sex_1)]<-.5 twins$sex_2[is.na(twins$sex_2)]<-.5 #__________________________ # fitting univariate saturated #__________________________ mzData1 <-subset(twins, zygosity==1 | zygosity==2, c(sex_1,age_1,sex_2, age_2,ia1_1,ia1_2)) dzData1 <-subset(twins, zygosity==3 | zygosity==4, c(sex_1,age_1,sex_2,age_2,ia1_1,ia1_2)) selVars <-c('ia1_1','ia1_2') nv <-1 ntv <-nv*2 nth <-1 thLabels <-'th1' thValues <-1.1 mzData1[,5:6] <-mxFactor(mzData1[,5:6], levels = c(0:1)) dzData1[,5:6] <-mxFactor(dzData1[,5:6], levels = c(0:1)) satOrdModel1 <-mxModel('satOrd1', mxModel('MZ', mxMatrix(type='Zero', nrow=nv, ncol=ntv, name='expMeanMZ'), # setting up matrices for covariate betas mxMatrix(type='Full', nrow=nth, ncol=2, free=TRUE, values=.1, label=c('bAge','bSex'), name='beta'), # matrix for observed age and sex covariates mxMatrix(type='Full', nrow=2, ncol=2, free=FALSE, label=c('data.age_1','data.sex_1','data.age_2','data.sex_2'), name='mzDefVars'), mxMatrix(type='Full', nrow=nth, ncol=ntv, free=TRUE, values=thValues, lbound=.001, name='thre'), mxMatrix(type='Lower', nrow=nth, ncol=nth, free=FALSE, values=1, name='inc'), mxAlgebra(expression= inc %*% thre, name='threInc'), mxAlgebra(expression=threInc, dimnames=list('th1', selVars), name='exthreMZ'), mxMatrix(type='Stand', nrow=ntv, ncol=ntv, free=TRUE, values=.5, name='expCovMZ'), # applying regression to the thresholds mxAlgebra(expression=exthreMZ+beta%*%mzDefVars, dimnames=list('th1', selVars), name='expThreshMZ'), mxData(observed=mzData1, type='raw'), mxFIMLObjective(covariance='expCovMZ', means='expMeanMZ', dimnames=selVars, thresholds='expThreshMZ') ), mxModel('DZ', mxMatrix(type='Zero', nrow=1, ncol=ntv, name='expMeanDZ'), mxMatrix(type='Stand', nrow=ntv, ncol=ntv, values=.3, free=TRUE, name='expCovDZ'), mxMatrix(type='Full', nrow=nth, ncol=ntv, free=TRUE, values=thValues, lbound=.001, name='thre'), mxMatrix(type='Lower', nrow=nth, ncol=nth, free=FALSE, values=1, name='inc'), mxAlgebra(expression= inc %*% thre, name='threInc'), mxAlgebra(expression=threInc, dimnames=list('th1', selVars), name='exthreDZ'), mxMatrix(type='Full', nrow=nth, ncol=2,free=TRUE, values=.1, label=c('bAge','bSex'), name='beta'), mxMatrix(type='Full', nrow=2, ncol=2, free=FALSE, label=c('data.age_1','data.sex_1','data.age_2','data.sex_2'), name='dzDefVars'), # applying regression to the thresholds mxAlgebra(expression=exthreDZ+beta%*%dzDefVars, dimnames=list('th1', selVars), name='expThreshDZ'), mxData(observed=dzData1, type='raw'), mxFIMLObjective(covariance='expCovDZ', means='expMeanDZ', dimnames=selVars, thresholds='expThreshDZ') ), mxAlgebra(MZ.objective + DZ.objective, name='minus2sumloglikelihood'), mxAlgebraObjective('minus2sumloglikelihood') ) satOrdModel1 <-mxOption(satOrdModel1, 'Function precision', 1e-9) satOrdFit1 <-mxRun(satOrdModel1) satOrdSumm1 <-summary(satOrdFit1) satOrdSumm1 #________________________ # running ACE threshold model code green # (for item 1 - makes careless mistakes) #________________________ aceOrdModel1 <-mxModel('aceOrd1', mxModel('ACE', # matrices a, c and e to store path coefficients mxMatrix(type='Lower', nrow=nv, ncol=nv, free=TRUE, values=.8, labels='a11', name='a'), mxMatrix(type='Lower', nrow=nv, ncol=nv, free=TRUE, values=.6, labels='c11', name='c'), mxMatrix(type='Lower', nrow=nv, ncol=nv, free=FALSE, values=1, labels='e11', name='e'), # matrices A, C and E to compute variance components mxAlgebra(expression=a %*% t(a), name='A'), mxAlgebra(expression=c %*% t(c), name='C'), mxAlgebra(expression=e %*% t(e), name='E'), # algebra to compute total variance and std. dev. (diagonal) mxAlgebra(expression=A+C+E, name='V'), mxMatrix(type='Iden', nrow=nv, ncol=nv, name='I'), mxAlgebra(expression=solve(sqrt(I*V)), name='iSD'), mxAlgebra(expression=A/V, name='heritability'), mxAlgebra(expression=C/V, name='comm_env'), mxAlgebra(expression=E/V, name='unique_env'), mxCI(c('heritability','comm_env','unique_env')), mxMatrix(type='Zero', nrow=1, ncol=ntv, name='expMean'), mxMatrix(type='Full', nrow=nth, ncol=nv, free=TRUE, values=thValues, lbound=.001, labels=thLabels, name='Thre'), mxMatrix(type='Lower', nrow=nth, ncol=nth, free=FALSE, values=1, name='Inc'), mxAlgebra(expression= Inc %*% Thre, name='ThreInc'), mxAlgebra(expression=cbind(ThreInc,ThreInc), dimnames=list('th1', selVars ), name='exthre'), # setting up matrices for covariate betas mxMatrix(type='Full', nrow=nth, ncol=2, free=TRUE, values=.1, label=c('bAge','bSex'), name='beta'), # algebra for expected variance /covariance matrix in MZ mxAlgebra(expression=rbind( cbind(A+C+E, A+C), cbind(A+C, A+C+E)), name='expCovMZ'), # algebra for exprected variance/covariance in DZ twins mxAlgebra(expression=rbind( cbind(A+C+E, 0.5%x%A+C), cbind(0.5%x%A+C, A+C+E)), name='expCovDZ') ), mxModel('MZ', mxData(observed=mzData1, type='raw'), # matrix for observed age and sex covariates mxMatrix(type='Full', nrow=2, ncol=2, free=FALSE, label=c('data.age_1','data.sex_1','data.age_2','data.sex_2'), name='mzDefVars'), # applying regression to the thresholds mxAlgebra(expression=ACE.exthre + ACE.beta%*%mzDefVars, dimnames=list('th1', selVars), name='expThreMZ'), mxFIMLObjective(covariance='ACE.expCovMZ', means='ACE.expMean', dimnames=selVars, thresholds='expThreMZ') ), mxModel('DZ', mxData(observed=dzData1, type='raw'), # matrix for observed age and sex covariates mxMatrix(type='Full', nrow=2, ncol=2, free=FALSE, label=c('data.age_1','data.sex_1','data.age_2','data.sex_2'), name='dzDefVars'), # applying regression to the thresholds mxAlgebra(expression=ACE.exthre + ACE.beta%*%dzDefVars, dimnames=list('th1', selVars), name='expThreDZ'), mxFIMLObjective(covariance='ACE.expCovDZ', means='ACE.expMean', dimnames=selVars, thresholds='expThreDZ') ), mxAlgebra(expression=MZ.objective + DZ.objective, name='minus2sumloglikelihood'), mxAlgebraObjective('minus2sumloglikelihood') ) aceOrdModel1 <-mxOption(aceOrdModel1, 'Function precision', 1e-9) aceOrdFit1 <-mxRun(aceOrdModel1, intervals=TRUE) aceOrdSum1 <-summary(aceOrdFit1) aceOrdSum1 ACEresults1 <-c(aceOrdFit1$ACE@algebras[c(6,7,8)]) ACEresults1 #_______________________________ # Now I fit a twin ADE model #_______________________________ adeOrdModel1 <-mxModel(aceOrdModel1, name='adeOrd1', mxModel(aceOrdFit1$ACE, mxMatrix(type='Full', nrow=1, ncol=1, free=FALSE, values=0, label='c11', name='c'), mxMatrix(type='Full', nrow=1, ncol=1, free=TRUE, values=.9, label='d11', name='d'), mxAlgebra(expression=d%*%t(d), name='D'), mxAlgebra(expression=A+D+E, name='V'), mxAlgebra(expression=rbind( cbind(A+D+E, A+D), cbind(A+D, A+D+E)), name='expCovMZ'), mxAlgebra(expression=rbind( cbind(A+D+E, 0.5%x%A+0.25%x%D), cbind(0.5%x%A+0.25%x%D, A+D+E)), name='expCovDZ'), mxAlgebra(expression=D/V, name='dominance'), mxCI(c('heritability','dominance','unique_env')) ) ) adeOrdFit1 <-mxRun(adeOrdModel1, intervals=TRUE) adeOrdSumm1 <-summary(adeOrdFit1) adeOrdSumm1 ADEresults1 <-c(adeOrdFit1$ACE@algebras[c(6,7,8,14)]) ADEresults1 #_______________________________ # Now I fit a twin AE model code green #_______________________________ aeOrdModel1 <-mxModel(aceOrdModel1, name='aeOrd1', mxModel(aceOrdFit1$ACE, mxMatrix(type='Full', nrow=1, ncol=1, free=FALSE, values=0, label='c11', name='c'), mxCI(c('heritability','comm_env','unique_env')) ) ) aeOrdFit1 <-mxRun(aeOrdModel1, intervals=TRUE) aeOrdSumm1 <-summary(aeOrdFit1) aeOrdSumm1 AEresults1 <-c(aeOrdFit1$ACE@algebras[c(6,7,8)]) AEresults1 #_______________________________ # Now I fit a twin E model code red #_______________________________ eOrdModel1 <-mxModel(aeOrdModel1, name='eOrd1', mxModel(aeOrdFit1$ACE, mxMatrix(type='Full', nrow=1, ncol=1, free=FALSE, values=0, label='a11', name='a'), mxCI(c('heritability','comm_env','unique_env')) ) ) eOrdFit1 <-mxRun(eOrdModel1, intervals=TRUE) eOrdSumm1 <-summary(eOrdFit1) eOrdSumm1 Eresults1 <-c(eOrdFit1$ACE@algebras[c(6,7,8)]) Eresults1 #_________________ # comparing models #_________________ nested.fit1 <- rbind( mxCompare(satOrdFit1, aceOrdFit1), mxCompare(satOrdFit1, adeOrdFit1)[2,], mxCompare(aceOrdFit1, aeOrdFit1), mxCompare(aeOrdFit1, eOrdFit1)[2,] ) nested.fit1