############################################################################################ #------------------------------------------------------------------------------------- # smoke and frequency of tea drinking # Bivariate Twin model for binary data # Phenotypes # Variable 1: tea_2g categorised in 2 classes: 0,1 # Variable 2: smk_cur categorised in 2 classes:0,1 # Type of data: Raw binary data # zygosity: 0=MZ 1=DZ # ------------------------------------------------------------------------------------- # Note: in reality you would rather analyse the continuous distributions rm(list=ls()) require(OpenMx) require(psych) library(foreign) setwd("C:/Users/dell/Desktop/data") # change optimizer ######## # mxOption(NULL, "Default optimizer", "NPSOL") mxOption( NULL, 'Default optimizer' ,'SLSQP' ) #mxOption(NULL, "Default optimizer", "CSOLNP") # I would suggest increasing mvnRelEps only if the model is taking too long to evaluate. # We generally want to make mvnRelEps closer to zero because that improves integration accuracy. mxOption( NULL, "mvnRelEps", value= mxOption(NULL, key="mvnRelEps")/5) getOption('mxOptions')$mvnRelEps data <- read.csv("data.csv") table(data$tea_2g1,data$tea_2g2) ########################################################### var( c( data$age1 , data$age2 ) ) me <- mean(c( data$age1 , data$age2 )) sde <- sd(c( data$age1 , data$age2 )) data$age1 <- (data$age1 - me)/sde data$age2 <- (data$age2 - me)/sde var( c( data$age1 , data$age2 ) ) ########################################################### vars <-c('tea_2g','smk_cur') covars <- c("age","region_type") nv <- 2 # number of variables per twin ncv <- 2 # number of covariates ntv <- nv*2 # number of variables per pair nth <- 1 # number of max thresholds = the number of category - 1 selVars <-paste(vars,c(rep(1,nv),rep(2,nv)),sep="") #c('tea2g1', 'smk_cur1', 'tea2g2', 'smk_cur2' )# covVars <- paste(covars,c(rep(1,ncv),rep(2,ncv)),sep="") mzData <- subset(data,zygosity1==0,c(selVars,covVars)) dzData <- subset(data,zygosity1==1,c(selVars,covVars)) sum(is.na(mzData)) sum(is.na(dzData)) mzData <- na.omit(mzData) dzData <- na.omit(dzData) dim(mzData); dim(dzData) #1333,605 # Generate Descriptive Statistics mzData[,selVars] <-mxFactor(mzData[,selVars], levels=c(0:nth) ) dzData[,selVars] <-mxFactor(dzData[,selVars], levels=c(0:nth) ) table( c(mzData[,1],mzData[,3]) , c(mzData[,4],mzData[,2])) table( c(dzData[,1],dzData[,3]) , c(dzData[,4],dzData[,2])) # Set Starting Values svTh <- c(1.1,1.5) # start value for thresholds,,get from tryhardordinal svCor <- .5 # start value for correlations svb <- 0.1 # sv for beta,get from tryhardordinal # Create Labels laThMZ <- paste("thMZ",selVars,sep="_") laThDZ <- paste("thDZ",selVars,sep="_") laThZ <- paste("thZ",selVars,sep="_") labdata <- cbind(paste("data.",covars[1],c(1,2),sep = ""),paste("data.",covars[2],c(1,2),sep = "")) #c("data.age1","data.age2","data.region_type1","data.region_type2") labeta <-cbind(paste("B",covars,"_",vars[1],sep = ""),paste("B",covars,"_",vars[2],sep = "")) #c('Bage_tea2g','Bregion_type_tea2g','Bage_smkcur', 'Bregion_type_smkcur') ########################################## # tea2g1 smkcur1 | tea2g2 smkcur2 # tea2g1 VP1T1 | # smkcur1 CVT1* VP2T1 | #----------------------------------------- # tea2g2 WP1* CTT2T1 | VP1T2 # smkcur2 CTT1T2* WP2* | CVT2 VP2T2 ########################################## # VP1T1(variance of pheno1 of twin1) VP2T1(variance of pheno2 of twin1) VP1T2(variance of pheno1 of twin2) VP2T2(variance of pheno2 of twin2) # CVT1(covar pheno1 pheno2 of twin1) CVT2(covar pheno1 pheno2 of twin2) # WP1(covar pheno1 of twin1 twin2) WP2(covar pheno2 of twin1 twin2) # CTT1T2(pheno1 of twin1 and pheno2 of twin2) CTT2T1(pheno1 of twin2 and pheno2 of twin1) laVaMZ <- c("MZCVT1", "MZWP1", "MZCTT1T2","MZCTT2T1", "MZWP2", "MZCVT2") # labels for (co)variances for MZ twins laVaDZ <- c("DZCVT1", "DZWP1", "DZCTT1T2","DZCTT2T1", "DZWP2", "DZCVT2") # labels for (co)variances for DZ twins # Define definition variables to hold the Covariates # Matrix & Algebra for expected means (SND), Thresholds, effect of Age on Th and correlations modelMZ <- mxModel( "MZ", mxMatrix( type="Full", nrow=2, ncol=ncv, free=F, label=labdata, name="MZDefVars"), mxMatrix( type="Full", nrow=ncv, ncol=nv, free=T, values=svb,labels=labeta, name="Beta"), mxMatrix( type="Zero", nrow=1, ncol=nv, name="meanG"), # or fix to zero" mxAlgebra( expression= cbind(meanG + MZDefVars[1,] %*% Beta,meanG + MZDefVars[2,] %*% Beta), name="expMeanMZ" ), mxMatrix( type="Full", nrow=nth, ncol=ntv, free=T, values=svTh, labels=laThMZ, name="ThMZ"), mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=svCor, labels=laVaMZ, lbound=-.99, ubound=.99, name="expCovMZ"), mxData( observed=mzData, type="raw" ), mxExpectationNormal( covariance="expCovMZ", means="expMeanMZ", dimnames=selVars, thresholds="ThMZ" ), mxFitFunctionML() ) modelDZ <- mxModel( "DZ", mxMatrix( type="Full", nrow=2, ncol=ncv, free=F, label=labdata, name="DZDefVars"), mxMatrix( type="Full", nrow=ncv, ncol=nv, free=T, values=svb,labels=labeta, name="Beta"), mxMatrix( type="Zero", nrow=1, ncol=nv, name="meanG"), mxAlgebra( expression= cbind(meanG + DZDefVars[1,] %*% Beta,meanG + DZDefVars[2,] %*% Beta), name="expMeanDZ" ), mxMatrix( type="Full", nrow=nth, ncol=ntv, free=T, values=svTh,labels=laThDZ, name="ThDZ"), mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=svCor, labels=laVaDZ, lbound=-.99, ubound=.99, name="expCovDZ"), mxData( observed=dzData, type="raw" ), mxExpectationNormal( covariance="expCovDZ", means="expMeanDZ", dimnames=selVars, thresholds="ThDZ" ), mxFitFunctionML() ) # Combine Groups Conf <- mxCI (c ('MZ.expCovMZ[3,1]','DZ.expCovDZ[3,1]','MZ.expCovMZ[4,2]','DZ.expCovDZ[4,2]','MZ.expCovMZ[2,1]','DZ.expCovDZ[2,1]','MZ.expCovMZ[4,1]','DZ.expCovDZ[4,1]') )#'MZ.expCovMZ[3,1]','MZ.expCovMZ[4,2]' SatModel <- mxModel( "Sat", modelMZ, modelDZ,mxFitFunctionMultigroup(c('MZ','DZ')), Conf ) # ------------------------------------------------------------------------------------------------------------------------------- # 1) RUN Saturated Model SatFit <- mxTryHardOrdinal(SatModel, intervals=F,extraTries = 15,OKstatuscodes=c(0,1,5),exhaustive=F) #,bestInitsOutput=T,silent=F,modify start values (SatSumm <- summary(SatFit)) # ----------------------x----------------------------------------------------------------------------- # 2) Run SubModels: # Test Significance of Covariates (result:no) dropb <- omxSetParameters( SatModel, label=labeta, free=FALSE, values=0,name="dropb" ) dropbFit <- mxTryHardOrdinal(dropb, intervals=F,extraTries = 15,OKstatuscodes=c(0,1,5),exhaustive=F) mxCompare(SatFit, dropbFit ) # Constrain expected Thresholds to be equal across twin order etop1Model <- mxModel( SatModel, name="etop1" ) etop1Model <- omxSetParameters( etop1Model, label=c(laThMZ[1],laThMZ[3]), free=TRUE, values=svTh[1], newlabels=laThMZ[1] ) etop1Model <- omxSetParameters( etop1Model, label=c(laThDZ[1],laThDZ[3]), free=TRUE, values=svTh[1], newlabels=laThDZ[1] ) etop1Fit <- mxTryHardOrdinal(etop1Model, intervals=F,extraTries = 15,OKstatuscodes=c(0,1,5),exhaustive=F) mxCompare(SatFit, subs <- list(dropbFit,etop1Fit ) ) (etop1Summ <- summary(etop1Fit)) #yes# etop2Model <- mxModel( etop1Model, name="etop2" ) etop2Model <- omxSetParameters( etop2Model, label=c(laThMZ[2],laThMZ[4]), free=TRUE, values=svTh[2], newlabels=laThMZ[2] ) etop2Model <- omxSetParameters( etop2Model, label=c(laThDZ[2],laThDZ[4]), free=TRUE, values=svTh[2], newlabels=laThDZ[2] ) etop2Fit <- mxTryHardOrdinal(etop2Model, intervals=F,extraTries = 15,OKstatuscodes=c(0,1,5),exhaustive=F) mxCompare(SatFit, subs <- list(dropbFit,etop1Fit,etop2Fit ) ) #yes# (etop2Summ <- summary(etop2Fit)) # Constrain expected Thresholds to be equal across zygosity etzp1Model <- mxModel( etop2Model, name="etzp1" ) etzp1Model <- omxSetParameters( etzp1Model, label=c(laThMZ[1],laThDZ[1]), free=TRUE, values=svTh[1], newlabels=laThZ[1] ) etzp1Fit <- mxTryHardOrdinal(etzp1Model, intervals=F,extraTries = 15,OKstatuscodes=c(0,1,5),exhaustive=F) mxCompare(SatFit, subs <- list(dropbFit,etop1Fit,etop2Fit,etzp1Fit) ) # yes# (etzp1Summ <- summary(etzp1Fit)) etzp2Model <- mxModel( etzp1Model, name="etzp2" ) etzp2Model <- omxSetParameters( etzp2Model, label=c(laThMZ[2],laThDZ[2]), free=TRUE, values=svTh[2], newlabels=laThZ[2] ) etzp2Fit <- mxTryHardOrdinal(etzp2Model, intervals=F,extraTries = 15,OKstatuscodes=c(0,1,5),exhaustive=F) mxCompare(SatFit, subs <- list(dropbFit,etop1Fit,etop2Fit,etzp1Fit,etzp2Fit) ) #YES # (etzp2Summ <- summary(etzp2Fit)) # equating variances ########################################## # tea2g1 smkcur1| tea2g2 smkcur2 # tea2g1 VP1T1 | # smkcur1 CVT1* VP2T1 | #----------------------------------------- # tea2g2 WP1* CTT2T1 | VP1T2 # smkcur2 CTT1T2* WP2* | CVT2 VP2T2 ########################################## #1 MZCVT1=MZCVT2 DZCVT1=DZCVT2 ecvoModel <- mxModel( etzp2Model, name="ecvo" ) ecvoModel <- omxSetParameters( ecvoModel, label=c("MZCVT1","MZCVT2"), newlabels='MZCV',free=TRUE) ecvoModel <- omxSetParameters( ecvoModel, label=c("DZCVT1","DZCVT2"), newlabels='DZCV',free=TRUE) ecvoFit <- mxTryHardOrdinal(ecvoModel, intervals=F,extraTries = 15,OKstatuscodes=c(0,1,5),exhaustive=F) mxCompare(SatFit,subs <- list(dropbFit,etop1Fit,etop2Fit,etzp1Fit,etzp2Fit,ecvoFit)) (ecvoSumm <- summary(ecvoFit)) #2 MZCTT1T2=MZCTT2T1 DZCTT1T2=DZCTT2T1 eCTCToModel <- mxModel( ecvoModel, name="eCTCTo" ) eCTCToModel <- omxSetParameters(eCTCToModel, label=c("MZCTT2T1","MZCTT1T2"),newlabels='MZCTCT' ,free=TRUE) eCTCToModel <- omxSetParameters(eCTCToModel, label=c("DZCTT2T1","DZCTT1T2"),newlabels='DZCTCT' ,free=TRUE) eCTCToFit<- mxTryHardOrdinal(eCTCToModel, intervals=T,extraTries = 15,OKstatuscodes=c(0,1,5),exhaustive=F) mxCompare(SatFit,subs <- list(dropbFit,etop1Fit,etop2Fit,etzp1Fit,etzp2Fit,ecvoFit,eCTCToFit)) (eCTCToSumm <- summary(eCTCToFit)) # report the tetrachoric correlations from the best-fitting saturated model tecor <- round( eCTCToSumm$CI[c(2,1,3)] , 3 ) for( i in 1:length(tecor[,1])){ tecor[i,1] <- paste(tecor[i,1],"(",tecor[i,2],",",tecor[i,3],")",sep = "") } tecor_vars <- data.frame(type=c("data","-"), zy=c("MZ","DZ"), tea2g=c(tecor[1,1],tecor[2,1]), smkcur=c(tecor[3,1],tecor[4,1]), CV=c(tecor[5,1],tecor[6,1]), CTCT=c(tecor[7,1],tecor[8,1])) sink('result/tecor_biv.csv',append = T) tecor_vars sink() ################################################################################ ### ACE-model ################################################################################ ### State the Quantitative genetic model lbPa <- .0001 # start value for lower bounds lbPaD <- diag(lbPa,nv,nv) # lower bounds for diagonal of covariance matrix lbPaD[lower.tri(lbPaD)] <- -10 # lower bounds for below diagonal elements lbPaD[upper.tri(lbPaD)] <- NA # lower bounds for above diagonal elements ### Some parameters included in all "submodels": baseACE <- mxModel('Base', mxMatrix( type="Full", nrow=ncv, ncol=nv, free=T, values=svb,labels=labeta, name="Beta"), mxMatrix( type="Zero", nrow=1, ncol=nv, name="meanG" ), # or type="Zero" # Create Matrices for Variance Components mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=c(2,0.1,2), name="A" ), # label=laLower("A",nv), values=valDiag(svPa,nv),# mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=c(2,0.1,2), name="C" ), mxMatrix( type="Symm", nrow=nv, ncol=nv, free=c(F,T,F), values=c(1, 0.1, 1), name="E" ),# fix E to avoid warning,without warning, E can be freely estimated mxAlgebra( expression=A+C+E, name="V" ), # Algebra to compute standardized variance components mxAlgebra( expression=A/V, name="h2"), mxAlgebra( expression=C/V, name="c2"), mxAlgebra( expression=E/V, name="e2"), # Algebra to compute total variances and standard deviations (tea2g_ggonal only) mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I"), mxAlgebra( expression=solve(sqrt(I*V)), name="iSD"), mxAlgebra( expression=solve(sqrt(I*V))%&%V, name ="Rph"), mxAlgebra( expression=solve(sqrt(I*A))%&%A, name ="Ra" ), #cov2cor() mxAlgebra( expression=solve(sqrt(I*C))%&%C, name ="Rc" ), mxAlgebra( expression=solve(sqrt(I*E))%&%E, name ="Re" ), # the proportion of ace contributing to the overlap # IF any r is negative,I agree that it makes more sense to report it as the contribution to the phenotypic correlation rather than # the proportion. So I will not divide by the phenotypic correlation! mxAlgebra( expression=sqrt(h2[1,1])*sqrt(h2[2,2])*Ra[1,2]/Rph[1,2],name = "pa"), mxAlgebra( expression=sqrt(c2[1,1])*sqrt(c2[2,2])*Rc[1,2]/Rph[1,2],name = "pc"), mxAlgebra( expression=sqrt(e2[1,1])*sqrt(e2[2,2])*Re[1,2]/Rph[1,2],name = "pe") ) modelMZ <- mxModel( "MZ", mxMatrix( type="Full", nrow=2, ncol=ncv, free=F, label=labdata, name="MZDefVars"), mxAlgebra( expression= cbind(Base.meanG + MZDefVars[1,] %*% Base.Beta,Base.meanG + MZDefVars[2,] %*% Base.Beta), name="expMeanMZ" ), mxMatrix( type="Full", nrow=nth, ncol=ntv, free=T, values=svTh, labels=laThZ, name="ThMZ"), mxAlgebra( expression= rbind( cbind(Base.A+Base.C+Base.E , Base.A+Base.C), cbind(Base.A+Base.C , Base.A+Base.C+Base.E)), name="expCovMZ" ), mxData( observed=mzData, type="raw" ), mxExpectationNormal( covariance="expCovMZ", means="expMeanMZ", dimnames=selVars, thresholds="ThMZ" ), mxFitFunctionML() ) modelDZ <- mxModel( "DZ", mxMatrix( type="Full", nrow=2, ncol=ncv, free=F, label=labdata, name="DZDefVars"), mxAlgebra( expression= cbind(Base.meanG + DZDefVars[1,] %*% Base.Beta,Base.meanG + DZDefVars[2,] %*% Base.Beta), name="expMeanDZ" ), mxMatrix( type="Full", nrow=nth, ncol=ntv, free=T, values=svTh,labels=laThZ, name="ThDZ"), mxAlgebra( expression= rbind( cbind(Base.A+Base.C+Base.E , 0.5%x%Base.A+Base.C), cbind(0.5%x%Base.A+Base.C , Base.A+Base.C+Base.E)), name="expCovDZ" ), mxData( observed=dzData, type="raw" ), mxExpectationNormal( covariance="expCovDZ", means="expMeanDZ", dimnames=selVars, thresholds="ThDZ" ), mxFitFunctionML() ) Conf1 <- mxCI (c ('Base.h2[1,1]','Base.c2[1,1]','Base.e2[1,1]','Base.h2[2,2]','Base.c2[2,2]','Base.e2[2,2]')) Conf2 <- mxCI (c ('Base.Ra', 'Base.Rc', 'Base.Re','Base.Rph') ) Conf3 <- mxCI (c ("Base.pa","Base.pc","Base.pe") ) AceModel <- mxModel( "ACE", baseACE,modelMZ, modelDZ, Conf1, Conf2,Conf3, mxFitFunctionMultigroup(c('MZ.fitfunction','DZ.fitfunction'))) # ------------------------------------------------------------------------------ # 4) RUN AceModel AceFit <- mxTryHardOrdinal(AceModel, intervals=F,extraTries = 15,OKstatuscodes=c(0,1,5),exhaustive=F) (AceSumm <- summary(AceFit)) mxCompare(SatFit,AceFit) AceFit$Base.Ra # ------------------------------------------------------------------------------------------ # RUN SUBMODELS # Test Significance of Covariates dropbace <- omxSetParameters( AceModel, label=labeta, free=FALSE, values=0,name="dropbace" ) dropbaceFit <- mxTryHardOrdinal(dropbace, intervals=F,extraTries = 15,OKstatuscodes=c(0,1,5),exhaustive=F) mxCompare( AceFit, dropbaceFit ) (dropbaceSumm <- summary(dropbaceFit)) #no# #C for p1 dropped AE-ACE dropcp1 <- omxSetParameters(AceModel,labels=c("Base.C[1,1]","Base.C[1,2]"),free=F,values=0,name="dropcp1") dropcp1Fit <- mxTryHardOrdinal(dropcp1, intervals=F,extraTries = 15,OKstatuscodes=c(0,1,5),exhaustive=F) mxCompare(AceFit,c(dropbaceFit,dropcp1Fit)) #no# (dropcp1Summ <- summary(dropcp1Fit)) #C for p2 dropped ACE-AE dropcp2 <- omxSetParameters(AceModel,labels=c("Base.C[2,2]","Base.C[1,2]"),free=F,values=0,name="dropcp2") #not sure# dropcp2Fit <- mxTryHardOrdinal(dropcp2, intervals=F,extraTries = 15,OKstatuscodes=c(0,1,5),exhaustive=F) mxCompare(AceFit,c(dropbaceFit,dropcp1Fit,dropcp2Fit)) #no# (dropcp2Summ <- summary(dropcp2Fit)) #C for p1p2 overlap dropped dropc12 <- omxSetParameters(AceModel,labels="Base.C[1,2]",free=F,values=0,name="dropc12") dropc12Fit <- mxTryHardOrdinal(dropc12, intervals=T,extraTries = 15,OKstatuscodes=c(0,1,5),exhaustive=F) mxCompare(AceFit,c(dropbaceFit,dropcp1Fit,dropcp2Fit,dropc12Fit)) #yes# (dropc12Summ <- summary(dropc12Fit)) #A for p1 dropped dropap1 <- omxSetParameters(dropc12,labels=c("Base.A[1,1]","Base.A[1,2]"),free=F,values=0,name="dropap1") dropap1Fit <- mxTryHardOrdinal(dropap1, intervals=F,extraTries = 15,OKstatuscodes=c(0,1,5),exhaustive=F) mxCompare(AceFit,c(dropbaceFit,dropcp1Fit,dropcp2Fit,dropc12Fit,dropap1Fit)) #no# (dropap1Summ <- summary(dropap1Fit)) #A for p2 dropped dropap2 <- omxSetParameters(dropc12,labels=c("Base.A[2,2]","Base.A[1,2]"),free=F,values=0,name="dropap2") dropap2Fit <- mxTryHardOrdinal(dropap2, intervals=F,extraTries = 15,OKstatuscodes=c(0,1,5),exhaustive=F) mxCompare(AceFit,c(dropbaceFit,dropcp1Fit,dropcp2Fit,dropc12Fit,dropap1Fit,dropap2Fit)) #NO# (dropap2Summ <- summary(dropap2Fit)) #A for p1p2 dropped dropa12 <- omxSetParameters(dropc12,labels="Base.A[1,2]",free=F,values=0,name="dropa12") dropa12Fit <- mxTryHardOrdinal(dropa12, intervals=F,extraTries = 15,OKstatuscodes=c(0,1,5),exhaustive=F) mxCompare(AceFit,c(dropbaceFit,dropcp1Fit,dropcp2Fit,dropc12Fit,dropap1Fit,dropap2Fit,dropa12Fit)) #NO# (dropa12Summ <- summary(dropa12Fit)) #E for p1p2 dropped drope12 <- omxSetParameters(dropc12,labels="Base.E[1,2]",free=F,values=0,name="drope12") drope12Fit <- mxTryHardOrdinal(drope12, intervals=T,extraTries = 15,OKstatuscodes=c(0,1,5),exhaustive=F) mxCompare(AceFit,c(dropbaceFit,dropcp1Fit,dropcp2Fit,dropc12Fit,dropap1Fit,dropap2Fit,dropa12Fit,drope12Fit)) #yes# (drope12Summ <- summary(drope12Fit)) #---------------------------------------------------- #report result com_sat <- mxCompare(SatFit, c(dropbFit,etop1Fit,etop2Fit,etzp1Fit,etzp2Fit,ecvoFit,eCTCToFit,AceFit) ) com_ace <- mxCompare(AceFit, c(dropbaceFit,dropcp1Fit,dropcp2Fit,dropc12Fit,dropap1Fit,dropap2Fit,dropa12Fit,drope12Fit)) com<- rbind(com_sat,com_ace[-c(1:2),]) com[,c(4,6,7)] <- round(com[,c(4,6,7)],2) com[,9] <- ifelse(com[,9]<0.001,"<0.001",round(com[,9],3)) com[1,1:2] <- vars sink('result/biv_fit_compare.csv',append = T) com sink() ace_full <- round( AceSumm$CI[,c(2,1,3)] , 2 ) ace_best <- round( dropc12Summ$CI[,c(2,1,3)] , 2 ) ace_j <-rbind(ace_full,ace_best) for( i in c(1:length(ace_j[,1]))){ ace_j[i,1] <- paste(ace_j[i,1],"(",ace_j[i,2],",",ace_j[i,3],")",sep = "") } ace_biv <- data.frame(a=c("tea2g","smkcur","-"), b=c("A","C","E"), tea1=c(ace_j[1,1],ace_j[2,1],ace_j[3,1]), smk1=c(ace_j[4,1],ace_j[5,1],ace_j[6,1]), r1=c(ace_j[7,1],ace_j[8,1],ace_j[9,1]), rph1=c(ace_j[10,1],"",""), pro1=c(ace_j[11,1],ace_j[12,1],ace_j[13,1]), c=rep("-",3), tea2=c(ace_j[14,1],ace_j[15,1],ace_j[16,1]), smk2=c(ace_j[17,1],ace_j[18,1],ace_j[19,1]), r2=c(ace_j[20,1],ace_j[21,1],ace_j[22,1]), rph2=c(ace_j[23,1],"-","-"), pro2=c(ace_j[24,1],ace_j[25,1],ace_j[26,1])) sink('result/ace_biv_ageregion.csv',append = T) ace_biv sink()