# ------------------------------------------------------------------------------ # Program: BivACE_OrdCon_covars.R (adapted from Hermine Mae's twoACEvj.R) # Author: Yi Ting Tan # Date: 28.01.2020 # # Bivariate ACE model with one continuous variable and one ordinal variable, as well as covariates. # -------|---------|---------|---------|---------|---------|---------|---------| #CHANGE THIS to the file destination for your output file # Load Libraries & Options #clearing the workspace and retaining only the dataset rm(list=setdiff(ls(),c("twinData2019","Vars","VarList","PVars","PVarList"))) 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") source('/Users/yitingtan/OneDrive - The University of Melbourne/Project 2017/Singtwins 2018 analysis/R codes 2019/GenEpiHelperFunctions.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/20200130bivACE_ConOrd.txt", "a") sink(outputFile, append=TRUE, split=TRUE) sink(outputFile, append=TRUE, type="message") mxOption(NULL,"Default optimizer","NPSOL") # ------------------------------------------------------------------ # PREPARE DATA #================================ # Select Variables for Analysis #CHANGE THESE to the two variables of interest to compare <- Vars and PVars Vars <- 'lnSingcomb.' PVars <- 'SIntSelf.' combVars <- c('Age.','Sex.','YrsEdNew.','AgeRelatedHearingLoss.','BilateralHearingLoss.',Vars,PVars) nvInit <- 7 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 oc <- c(0,1) # 0 for continuous, 1 for ordinal variables # 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$lnSingcomb.1)&!is.na(twinData2019$lnSingcomb.2)&!is.na(twinData2019$SIntSelf.1)&!is.na(twinData2019$SIntSelf.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 #filtering out twins likely to have MTNDeviceProb #twinData <- twinData[twinData$MTNDeviceProb.1!=1 & twinData$MTNDeviceProb.2!=1,] 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','BHearing1','Vars1','PVars1','Age2','Sex2','yrsEd2','AHearing2','BHearing2','Vars2','PVars2') mzDataF <- cbind(mzData[,c('Age1','Sex1','yrsEd1','AHearing1','BHearing1','Vars1','Age2','Sex2','yrsEd2','AHearing2','BHearing2','Vars2')],mxFactor( x=mzData[,c('PVars1','PVars2')], levels=c(1:7)) ) dzDataF <- cbind(dzData[,c('Age1','Sex1','yrsEd1','AHearing1','BHearing1','Vars1','Age2','Sex2','yrsEd2','AHearing2','BHearing2','Vars2')],mxFactor( x=dzData[,c('PVars1','PVars2')], 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 the ordinal variable up for OpenMx # mzData$PVars1 <- mxFactor(mzData$PVars1, levels=c(1:7) ) # mzData$PVars2 <- mxFactor(mzData$PVars2, levels=c(1:7) ) # dzData$PVars1 <- mxFactor(dzData$PVars1, levels=c(1:7) ) # dzData$PVars2 <- mxFactor(dzData$PVars2, levels=c(1:7) ) # Set Starting Values frMV <- c(TRUE,FALSE) # free status for variables svMe <- c(3,3) # start value for means svPa <- .2 # start value for path coefficient #svPe <- .5 # start value for path coefficient for e svLTh <- -.5 # start value for first threshold svITh <- 1 # 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 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_SexOrd <- mxMatrix( type="Full", nrow=nth, ncol=1, free=TRUE, values= .01, labels="betaSexOrd", name="bSexOrd") B_SexCon <- mxMatrix( type="Full", nrow=1, ncol=2, free=c(T,F), values= c(.01, 0), labels=c('bSexV1','bSexV2'), name="bSexCon") meanSexOrd <- mxAlgebra( bSexOrd%x%Sex, name="SexROrd") meanSexCon <- mxAlgebra( bSexCon%x%Sex, name="SexRCon") #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_AgeOrd <- mxMatrix( type="Full", nrow=nth, ncol=1, free=FALSE, values= 0, labels="betaAgeOrd", name="bAgeOrd") B_AgeCon <- mxMatrix( type="Full", nrow=1, ncol=2, free=c(T,F), values= c(.01,0), labels=c('bAgeV1','bAgeV2'), name="bAgeCon") meanAgeOrd <- mxAlgebra( bAgeOrd%x%Age, name="AgeROrd") meanAgeCon <- mxAlgebra( bAgeCon%x%Age, name="AgeRCon") #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_YEdOrd <- mxMatrix( type="Full", nrow=nth, ncol=1, free=FALSE, values= 0, labels="betaYEdOrd", name="bYEdOrd") B_YEdCon <- mxMatrix( type="Full", nrow=1, ncol=2, free=c(T,F), values= c(.01, 0), labels=c('bYEdV1','bYEdV2'), name="bYEdCon") meanYEdOrd <- mxAlgebra( bYEdOrd%x%YEd, name="YEdROrd") meanYEdCon <- mxAlgebra( bYEdCon%x%YEd, name="YEdRCon") #Age-related hearing condition defAHearing <- mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, labels=c("data.AHearing1","data.AHearing2"), name="AHearing") # Matrices declared to store linear Coefficients for covariate B_AHearingOrd <- mxMatrix( type="Full", nrow=nth, ncol=1, free=TRUE, values= .01, labels="betaAHearOrd", name="bAHearingOrd") B_AHearingCon <- mxMatrix( type="Full", nrow=1, ncol=2, free=c(T,F), values= c(.01, 0), labels=c('bAHearV1','bAHearV2'), name="bAHearingCon") meanAHearingOrd <- mxAlgebra( bAHearingOrd%x%AHearing, name="AHearingROrd") meanAHearingCon <- mxAlgebra( bAHearingCon%x%AHearing, name="AHearingRCon") #Bilateral hearing condition defBHearing <- mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, labels=c("data.BHearing1","data.BHearing2"), name="BHearing") # Matrices declared to store linear Coefficients for covariate B_BHearingOrd <- mxMatrix( type="Full", nrow=nth, ncol=1, free=FALSE, values= 0, labels="betaBHearOrd", name="bBHearingOrd") B_BHearingCon <- mxMatrix( type="Full", nrow=1, ncol=2, free=c(T,F), values= c(.01, 0), labels=c('bBHearV1','bBHearV2'), name="bBHearingCon") meanBHearingOrd <- mxAlgebra( bBHearingOrd%x%BHearing, name="BHearingROrd") meanBHearingCon <- mxAlgebra( bBHearingCon%x%BHearing, name="BHearingRCon") # Matrix & Algebra for expected means vector and expected thresholds intercept <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=c(T,F), values=c(3,0), labels=c("meanP","binary"), name="intercept" ) threG <- mxMatrix( type="Full", nrow=nth, ncol=nv, 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 + AgeROrd + SexROrd + YEdROrd + AHearingROrd + BHearingROrd, name = "expThreC") #with covariates expMean <- mxAlgebra( intercept + AgeRCon + SexRCon + YEdRCon + AHearingRCon + BHearingRCon, name="expMean") inclusions <- list (defSex, B_SexOrd, B_SexCon, meanSexOrd, meanSexCon, defAge, B_AgeOrd, B_AgeCon, meanAgeOrd, meanAgeCon, defYEd, B_YEdOrd, B_YEdCon, meanYEdOrd, meanYEdCon, defAHearing, B_AHearingOrd, B_AHearingCon, meanAHearingOrd, meanAHearingCon, defBHearing, B_BHearingOrd, B_BHearingCon, meanBHearingOrd, meanBHearingCon, expMean, intercept, threG, threT, threC) # ------------------------------------------------------------------------------ # 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" ) covD <- mxAlgebra( expression=d %*% t(d), name="D" ) covC <- mxAlgebra( expression=c %*% t(c), name="C" ) covE <- mxAlgebra( expression=e %*% t(e), name="E" ) # Algebra to compute total variances and standard deviations (diagonal only) covP <- mxAlgebra( expression=A+D+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 # Constrain Variance of Binary Variables matUnv <- mxMatrix( type="Unit", nrow=1, ncol=1, name="Unv1" ) matOc <- mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, values=c(0,1), name="Oc" ) var1 <- mxConstraint( expression=diag2vec(Oc%&%V)==Unv1, name="Var1" ) # Algebras generated to hold Parameter Estimates and Derived Variance Components rowVars <- rep('vars',nv) colVars <- rep(c('A','D','C','E','SA','SD','SC','SE'),each=nv) #each variable is repeated, A of Var1, A of var2, D of Var1, D of Var2... estVars <- mxAlgebra( expression=cbind(A,D,C,E,A/V,D/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") # Data objects for Multiple Groups dataMZ <- mxData( observed=mzDataF, type="raw" ) dataDZ <- mxData( observed=dzDataF, type="raw" ) covMZ <- mxAlgebra( expression= rbind( cbind(V , A+D+C), cbind(A+D+C , V)), name="expCovMZ" ) covDZ <- mxAlgebra( expression= rbind( cbind(V, 0.5%x%A+0.25%x%D+C), cbind(0.5%x%A+0.25%x%D+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")) # Objective objects for Multiple Groups expMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=c('Vars1','PVars1','Vars2','PVars2'), thresholds="expThreC", threshnames=c('PVars1','PVars2') ) expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=c('Vars1','PVars1','Vars2','PVars2'), thresholds="expThreC", threshnames=c('PVars1','PVars2') ) funML <- mxFitFunctionML() # Combine Groups pars <- list( pathA, pathD, pathC, pathE, SA,SD,SC,SE,covA, covD, covC, covE, covP, matUnv, matOc, matI, invSD, estVars, corA, corC, corE, corrA, corrC, corrE, inc) modelMZ <- mxModel( pars, inclusions, covMZ, dataMZ, expMZ, funML, name="MZ" ) modelDZ <- mxModel( pars, inclusions, 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) } ) #to avoid the optimality prob. #BivAceFitciAgain <- mxRun(BivAceFitci,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=""))