# ------------------------------------------------------------------------------ # Program: BivACE_Ord_ASEAh.R # Author: Hermine Maes # Date: 09 08 2015 # # Bivariate Twin Analysis model to estimate causes of variation # Matrix style model - Raw data - Ordinal data # -------|---------|---------|---------|---------|---------|---------|---------| # Load Libraries & Options #clearing the workspace and retaining only the dataset rm(list=setdiff(ls(),c("twinData2019","Vars","VarList"))) library(OpenMx) library(psych); library(polycor) #mac source for mi function source("/Users/yitingtan/OneDrive - The University of Melbourne/Project 2017/Singtwins 2018 analysis/R codes 2019/miFunctions.R") #PC source for mi function #source("C:/Users/YiTing/OneDrive - The University of Melbourne/Project 2017/Singtwins 2018 analysis/R codes 2019/miFunctions.R") #Mac source for output file outputFile <- file("/Users/yitingtan/OneDrive - The University of Melbourne/Project 2017/Singtwins 2018 analysis/R output/20200130_BivACEord_misimasa.txt", "a") #PC source for output file #outputFile <- file("C:/Users/YiTing/OneDrive - The University of Melbourne/Project 2017/Singtwins 2018 analysis/R output/20191210Sat_ASE_lnVarsHet5WithMsg.txt") sink(outputFile, append=TRUE, split=TRUE) sink(outputFile, append=TRUE, type="message") mxOption(NULL,"Default optimizer","NPSOL") #mxOption(NULL,"Default optimizer","SLSQP") # ------------------------------------------------------------------ # PREPARE DATA #================================ # Select Variables for Analysis #CHANGE THESE to the two variables of interest to compare <- Vars and PVars Vars <- 'SAbSelf.' PVars <- 'SIntSelf.' combVars <- c('Age.','Sex.','YrsEdNew.','AgeRelatedHearingLoss.',Vars,PVars) nvInit <- 6 nv <- 2 # number of variables ntv <- nv*2 # number of total variables selVars <- paste(combVars,c(rep(1,nvInit),rep(2,nvInit)),sep="") #eg c('ht1','wt1,'ht2','wt2') nth <- 6 # number of thresholds (for likert scale of 7, there will be 7-1=6 thresholds; binary it will be 1 threshold) ; change accordingly # Select Covariates for Analysis twinData <- subset(twinData2019, (!is.na(twinData2019$Age.1) & !is.na(twinData2019$Age.2) &!is.na(twinData2019$Sex.1)& !is.na(twinData2019$Sex.2) &!is.na(twinData2019$YrsEdNew.1)& !is.na(twinData2019$YrsEdNew.2) &!is.na(twinData2019$AgeRelatedHearingLoss.1)& !is.na(twinData2019$AgeRelatedHearingLoss.2) #&!is.na(twinData2019$BilateralHearingLoss.1)& !is.na(twinData2019$BilateralHearingLoss.2) #&!is.na(twinData2019$touchscreenDummy.1)& !is.na(twinData2019$touchscreenDummy.2) #&!is.na(twinData2019$MTNearphone.1)& !is.na(twinData2019$MTNearphone.2) )) twinData[,'Age.1'] <- twinData[,'Age.1']/100 twinData[,'Age.2'] <- twinData[,'Age.2']/100 twinData[,'YrsEdNew.1'] <- twinData[,'YrsEdNew.1']/100 twinData[,'YrsEdNew.2'] <- twinData[,'YrsEdNew.2']/100 cat("\n") cat("----------------------\n") cat("MZ v DZall\n") Vars PVars cat("----------------------\n") mzData <- subset(twinData, DeterminedZyggroup==1|DeterminedZyggroup==3, selVars) #change accordingly. can be ZygosityGrp too dzData <- subset(twinData, DeterminedZyggroup==2|DeterminedZyggroup==4|DeterminedZyggroup==5, selVars) colnames(mzData)=colnames(dzData)=c('Age1','Sex1','yrsEd1','AHearing1','Vars1','PVars1','Age2','Sex2','yrsEd2','AHearing2','Vars2','PVars2') # mzDataF <- cbind(mzData[,c('Age1','Sex1','yrsEd1','AHearing1','Age2','Sex2','yrsEd2','AHearing2')],mxFactor( x=mzData[,c('Vars1','PVars1','Vars2','PVars2')], levels=c(1:7)) ) # dzDataF <- cbind(dzData[,c('Age1','Sex1','yrsEd1','AHearing1','Age2','Sex2','yrsEd2','AHearing2')],mxFactor( x=dzData[,c('Vars1','PVars1','Vars2','PVars2')], levels=c(1:7)) ) #change accordingly. The column numbers of t1 and t2; as well as the level values. mzData[,5] <- mxFactor( x=mzData[,5], levels=c(1:7)) #misimasa c(1:7) mzData[,11] <- mxFactor( x=mzData[,11], levels=c(1:7)) mzData[,6] <- mxFactor( x=mzData[,6], levels=c(1:7) ) #misimasa c(1:7) mzData[,12] <- mxFactor( x=mzData[,12], levels=c(1:7) ) dzData[,5] <- mxFactor( x=dzData[,5], levels=c(1:7) ) dzData[,11] <- mxFactor( x=dzData[,11], levels=c(1:7) ) dzData[,6] <- mxFactor( x=dzData[,6], levels=c(1:7) ) dzData[,12] <- mxFactor( x=dzData[,12], levels=c(1:7) ) cat("number of valid MZ cases with non-missing Vars and PVars:") sum(!is.na(mzData$Vars1)&!is.na(mzData$Vars2)&!is.na(mzData$PVars1)&!is.na(mzData$PVars2)) cat("number of valid DZ cases with non-missing Vars and PVars:") sum(!is.na(dzData$Vars1)&!is.na(dzData$Vars2)&!is.na(dzData$PVars1)&!is.na(dzData$PVars2)) cat("\n") describe(mzData) describe(dzData) cat("\n") # Set Starting Values svLTh <- -2 # start value for first threshold svITh <- .5 # start value for increments svTh <- matrix(rep(c(svLTh,(rep(svITh,nth-1)))),nrow=nth,ncol=nv) # start value for thresholds lbTh <- matrix(rep(c(-3,(rep(0.001,nth-1))),nv),nrow=nth,ncol=nv) # lower bounds for thresholds svPa <- .2 # start value for path coefficient #svPe <- .5 # start value for path coefficient for e labThZ <- c(paste("t",1:nth,"Z",sep=""),paste("t",1:nth,"Z",sep="")) laMe <- c('Vars1_Mean','PVars1_Mean','Vars2_Mean','PVars2_Mean') svPas <- diag(svPa,nv,nv) laA <- paste("a",rev(nv+1-sequence(1:nv)),rep(1:nv,nv:1),sep="") # c("a11","a21","a22") laD <- paste("d",rev(nv+1-sequence(1:nv)),rep(1:nv,nv:1),sep="") laC <- paste("c",rev(nv+1-sequence(1:nv)),rep(1:nv,nv:1),sep="") laE <- paste("e",rev(nv+1-sequence(1:nv)),rep(1:nv,nv:1),sep="") # ------------------------------------------------------------------------------ # PREPARE MODEL # Matrix for moderating/interacting variable defSex <- mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, labels=c("data.Sex1","data.Sex2"), name="Sex") # Matrices declared to store linear Coefficients for covariate B_Sex <- mxMatrix( type="Full", nrow=nth, ncol=2, free=c(F,T), values= c(0,.1), label=c('bSexV1','bSexV2'), name="bSex", byrow = TRUE) meanSex <- mxAlgebra( bSex%x%Sex, name="SexR") #age defAge <- mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, labels=c("data.Age1","data.Age2"), name="Age") # Matrices declared to store linear Coefficients for covariate B_Age <- mxMatrix( type="Full", nrow=nth, ncol=2, free=c(T,F), values= c(.01,0), label=c('bAgeV1','bAgeV2'), name="bAge", byrow = TRUE) meanAge <- mxAlgebra( bAge%x%Age, name="AgeR") #YrsEd defYEd <- mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, labels=c("data.yrsEd1","data.yrsEd2"), name="YEd") # Matrices declared to store linear Coefficients for covariate B_YEd <- mxMatrix( type="Full", nrow=nth, ncol=2, free=c(T,F), values= c(.01,0), labels=c('bYEdV1','bYEdV2'), name="bYEd", byrow = TRUE) meanYEd <- mxAlgebra( bYEd%x%YEd, name="YEdR") #Age-related hearing condition defHearing <- mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, labels=c("data.AHearing1","data.AHearing2"), name="Hearing") # Matrices declared to store linear Coefficients for covariate B_Hearing <- mxMatrix( type="Full", nrow=nth, ncol=2, free=TRUE, values= .01, labels=c('bAHearV1','bAHearV2'), name="bHearing", byrow = TRUE) meanHearing <- mxAlgebra( bHearing%x%Hearing, name="HearingR") #*************************************************************** defs <- list( defSex, B_Sex, meanSex, defAge, B_Age, meanAge, defYEd, B_YEd, meanYEd, defHearing, B_Hearing, meanHearing) #setting up the regression # Matrix & Algebra for expected means vector and expected thresholds meanG <- mxMatrix( type="Zero", nrow=1, ncol=ntv, name="expMean" ) threG <- mxMatrix( type="Full", nrow=nth, ncol=ntv, free=TRUE, values=svTh, lbound=lbTh, labels=labThZ, name="Thre" ) inc <- mxMatrix( type="Lower", nrow=nth, ncol=nth, free=FALSE, values=1, name="Inc" ) threT <- mxAlgebra( expression = Inc %*% Thre, name="expThre" ) threC <- mxAlgebra( expression = expThre + AgeR + SexR + YEdR + HearingR, name = "expThreC") #with covariates # ACE Model # Matrices declared to store a, c, and e Path Coefficients pathA <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPas, label=laA, name="a" ) #pathD <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=FALSE, values=0, label=laD, name="d" ) pathC <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPas, label=laC, name="c" ) pathE <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPas, label=laE, 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" ) covP <- mxAlgebra( expression=A+C+E, name="V" ) matI <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I") invSD <- mxAlgebra( expression=solve(sqrt(I*V)), name="iSD") #this is to help to get stdized paths # Algebras generated to hold Parameter Estimates and Derived Variance Components rowVars <- rep('vars',nv) colVars <- rep(c('A','C','E','SA','SC','SE'),each=nv) estVars <- mxAlgebra( expression=cbind(A,C,E,A/V,C/V,E/V), name="Vars", dimnames=list(rowVars,colVars)) SA <- mxAlgebra(expression=A/V, name="SA") SD <- mxAlgebra(expression=D/V, name="SD") SC <- mxAlgebra(expression=C/V, name="SC") SE <- mxAlgebra(expression=E/V, name="SE") # Algebra for expected Mean and Variance/Covariance Matrices in MZ & DZ twins covMZ <- mxAlgebra( expression= rbind( cbind(V, A+C), cbind(A+C, V)), name="expCovMZ" ) covDZ <- mxAlgebra( expression= rbind( cbind(V, 0.5%x%A+C), cbind(0.5%x%A+C, V)), name="expCovDZ" ) # Calculate genetic and environmental correlations corA <- mxAlgebra( expression=solve(sqrt(I*A))%&%A, name ="rA" ) corC <- mxAlgebra( expression=solve(sqrt(I*C))%&%C, name ="rC" ) corE <- mxAlgebra( expression=solve(sqrt(I*E))%&%E, name ="rE" ) #calculate genetic and environmental correlations due to A, D, C and E corrA <- mxAlgebra(expression = sqrt(SA[1,1])*rA[2,1]*sqrt(SA[2,2]), name="corrA") corrC <- mxAlgebra(expression = sqrt(SC[1,1])*rC[2,1]*sqrt(SC[2,2]), name="corrC") corrE <- mxAlgebra(expression = sqrt(SE[1,1])*rE[2,1]*sqrt(SE[2,2]), name="corrE") ci <- mxCI (c("SA[1,1]","SA[2,2]","SD[1,1]","SD[2,2]","SC[1,1]","SC[2,2]","SE[1,1]","SE[2,2]","rA[2,1]","rC[2,1]","rE[2,1]","corrA","corrC","corrE")) # Constraint on variance of Ordinal variables matUnv <- mxMatrix( type="Unit", nrow=nv, ncol=1, name="Unv1" ) var1 <- mxConstraint( expression=diag2vec(V)==Unv1, name="Var1" ) # Create Data Objects for Multiple Groups dataMZ <- mxData( observed=mzData, type="raw" ) dataDZ <- mxData( observed=dzData, type="raw" ) # Create Expectation Objects for Multiple Groups expMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=c('Vars1','PVars1','Vars2','PVars2'), thresholds="expThreC" ) expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=c('Vars1','PVars1','Vars2','PVars2'), thresholds="expThreC" ) funML <- mxFitFunctionML() # Combine Groups pars <- list( pathA, pathC, pathE, SA,SC,SE,covA, covC, covE, covP, matUnv, matI, invSD, estVars, corA, corC, corE, corrA, corrC, corrE, inc) modelMZ <- mxModel( pars, defs, meanG, threG, threT, threC, covMZ, dataMZ, expMZ, funML, name="MZ" ) modelDZ <- mxModel( pars, defs, meanG, threG, threT, threC, covDZ, dataDZ, expDZ, funML, name="DZ" ) fitML <- mxFitFunctionMultigroup(c("MZ.fitfunction","DZ.fitfunction")) BivAceModel <- mxModel( "BivACE", pars, var1, modelMZ, modelDZ, fitML, ci) # ------------------------------------------------------------------------------ # RUN MODEL # Run Bivariate ACE model BivAceFit <- tryCatch( expr = { mxRun( BivAceModel, intervals=T ) }, warning = function(w){ mxTryHard(BivAceModel, extraTries = 20,bestInitsOutput=T, intervals=T) } ) BivAceSumm <- summary(BivAceFit) BivAceSumm # Generate List of Parameter Estimates and Derived Quantities using formatOutputMatrices ACEpathMatrices <- c("a","c","e","iSD","iSD %*% a","iSD %*% c","iSD %*% e","(iSD %*% a)*(iSD %*% a)","(iSD %*% c)*(iSD %*% c)","(iSD %*% e)*(iSD %*% e)") ACEpathLabels <- c("pathEst_a","pathEst_d","pathEst_e","isd","stPathEst_a","stPathEst_c","stPathEst_e", "st2PathEst_a","st2PathEst_c","st2PathEst_e") formatOutputMatrices(BivAceFit,ACEpathMatrices,ACEpathLabels,Vars,4) ACEcovMatrices <- c("A","C","E","V","A/V","C/V","E/V") ACEcovLabels <- c("covComp_A","covComp_C","covComp_E","Var","stCovComp_A","stCovComp_C","stCovComp_E") formatOutputMatrices(BivAceFit,ACEcovMatrices,ACEcovLabels,Vars,4) ACEcorMatrices <- c("rA","rC","rE","corrA","corrC","corrE") ACEcorLabels <- c("rA","rC","rE","corrA","corrC","corrE") formatOutputMatrices(BivAceFit,ACEcorMatrices,ACEcorLabels,Vars,4) #drop a21, c21, and e21 respectively BivAceModelna21 <- mxModel( BivAceFit, name="ACEna21" ) BivAceModelna21 <- omxSetParameters( BivAceModelna21, labels=c("a21"), free=FALSE, values=0 ) #BivAcena21Fit <- mxRun(BivAceModelna21, intervals=T) BivAcena21Fit <- tryCatch( expr = { mxRun( BivAceModelna21, intervals=T ) }, warning = function(w){ mxTryHard(BivAceModelna21, extraTries = 20,bestInitsOutput=T, intervals=T) } ) BivAcena21Summ <- summary(BivAcena21Fit) BivAceModelnc21 <- mxModel( BivAceFit, name="ACEnc21" ) BivAceModelnc21 <- omxSetParameters( BivAceModelnc21, labels=c("c21"), free=FALSE, values=0 ) #BivAcenc21Fit <- mxRun(BivAceModelnc21, intervals=T) BivAcenc21Fit <- tryCatch( expr = { mxRun( BivAceModelnc21, intervals=T ) }, warning = function(w){ mxTryHard(BivAceModelnc21, extraTries = 20,bestInitsOutput=T, intervals=T) } ) BivAcenc21Summ <- summary(BivAcenc21Fit) BivAceModelne21 <- mxModel( BivAceFit, name="ACEne21" ) BivAceModelne21 <- omxSetParameters( BivAceModelne21, labels=c("e21"), free=FALSE, values=0 ) #BivAcene21Fit <- mxRun(BivAceModelne21, intervals=T) BivAcene21Fit <- tryCatch( expr = { mxRun( BivAceModelne21, intervals=T ) }, warning = function(w){ mxTryHard(BivAceModelne21, extraTries = 20,bestInitsOutput=T, intervals=T) } ) BivAcene21Summ <- summary(BivAcene21Fit) BivAceNested <- list(BivAcena21Fit, BivAcenc21Fit, BivAcene21Fit) tableFitStatistics(BivAceFit,BivAceNested) #find the most parsimonious model based on lowest AIC BivAcena21Summ # Generate List of Parameter Estimates and Derived Quantities using formatOutputMatrices ACEpathMatrices <- c("a","c","e","iSD","iSD %*% a","iSD %*% c","iSD %*% e","(iSD %*% a)*(iSD %*% a)","(iSD %*% c)*(iSD %*% c)","(iSD %*% e)*(iSD %*% e)") ACEpathLabels <- c("pathEst_a","pathEst_d","pathEst_e","isd","stPathEst_a","stPathEst_c","stPathEst_e", "st2PathEst_a","st2PathEst_c","st2PathEst_e") formatOutputMatrices(BivAcena21Fit,ACEpathMatrices,ACEpathLabels,Vars,4) ACEcovMatrices <- c("A","C","E","V","A/V","C/V","E/V") ACEcovLabels <- c("covComp_A","covComp_C","covComp_E","Var","stCovComp_A","stCovComp_C","stCovComp_E") formatOutputMatrices(BivAcena21Fit,ACEcovMatrices,ACEcovLabels,Vars,4) ACEcorMatrices <- c("rA","rC","rE","corrA","corrC","corrE") ACEcorLabels <- c("rA","rC","rE","corrA","corrC","corrE") formatOutputMatrices(BivAcena21Fit,ACEcorMatrices,ACEcorLabels,Vars,4) BivAcenc21Summ # Generate List of Parameter Estimates and Derived Quantities using formatOutputMatrices ACEpathMatrices <- c("a","c","e","iSD","iSD %*% a","iSD %*% c","iSD %*% e","(iSD %*% a)*(iSD %*% a)","(iSD %*% c)*(iSD %*% c)","(iSD %*% e)*(iSD %*% e)") ACEpathLabels <- c("pathEst_a","pathEst_d","pathEst_e","isd","stPathEst_a","stPathEst_c","stPathEst_e", "st2PathEst_a","st2PathEst_c","st2PathEst_e") formatOutputMatrices(BivAcenc21Fit,ACEpathMatrices,ACEpathLabels,Vars,4) ACEcovMatrices <- c("A","C","E","V","A/V","C/V","E/V") ACEcovLabels <- c("covComp_A","covComp_C","covComp_E","Var","stCovComp_A","stCovComp_C","stCovComp_E") formatOutputMatrices(BivAcenc21Fit,ACEcovMatrices,ACEcovLabels,Vars,4) ACEcorMatrices <- c("rA","rC","rE","corrA","corrC","corrE") ACEcorLabels <- c("rA","rC","rE","corrA","corrC","corrE") formatOutputMatrices(BivAcenc21Fit,ACEcorMatrices,ACEcorLabels,Vars,4) BivAcene21Summ # Generate List of Parameter Estimates and Derived Quantities using formatOutputMatrices ACEpathMatrices <- c("a","c","e","iSD","iSD %*% a","iSD %*% c","iSD %*% e","(iSD %*% a)*(iSD %*% a)","(iSD %*% c)*(iSD %*% c)","(iSD %*% e)*(iSD %*% e)") ACEpathLabels <- c("pathEst_a","pathEst_d","pathEst_e","isd","stPathEst_a","stPathEst_c","stPathEst_e", "st2PathEst_a","st2PathEst_c","st2PathEst_e") formatOutputMatrices(BivAcene21Fit,ACEpathMatrices,ACEpathLabels,Vars,4) ACEcovMatrices <- c("A","C","E","V","A/V","C/V","E/V") ACEcovLabels <- c("covComp_A","covComp_C","covComp_E","Var","stCovComp_A","stCovComp_C","stCovComp_E") formatOutputMatrices(BivAcene21Fit,ACEcovMatrices,ACEcovLabels,Vars,4) ACEcorMatrices <- c("rA","rC","rE","corrA","corrC","corrE") ACEcorLabels <- c("rA","rC","rE","corrA","corrC","corrE") formatOutputMatrices(BivAcene21Fit,ACEcorMatrices,ACEcorLabels,Vars,4) # ---------------------------------------------------------------------------------------------------------------------- sink(type="message") sink() #save.image(paste(filename,".Ri",sep=""))