rm(list=ls()) #load packages library(readr) library(plyr) library(psych) library("OpenMx") ### Specify optimizer mxOption(NULL, "Default optimizer", "SLSQP") #load data mtfs <- read.csv("/Users/zelle063/Desktop/Vrieze Lab/GEDI_Cleaning/CleanData/SUDEV.csv") #create cohort variables mtfs$OC[mtfs$IDAC==17] <- 1 mtfs$OC[mtfs$IDAC==11] <- 0 mtfs$YC[mtfs$IDAC==17] <- 0 mtfs$YC[mtfs$IDAC==11 & mtfs$IDES == 9] <- 1 mtfs$YC[mtfs$IDAC==11 & mtfs$IDES != 9] <- 0 mtfs$ES[mtfs$IDES == 9] <- 0 mtfs$ES[mtfs$IDES != 9] <- 1 ##reshape to wide format mtfs <- subset(mtfs, select=c("IDYRFAM", "IDSC", "IDSEX", "ZYGOSITY", "AlcFreq14", "AlcFreq17", "AlcFreq20", "AlcFreq24", "AlcFreq29", "AlcFreq35", "AlcFreq45", "AlcQuant14", "AlcQuant17", "AlcQuant20", "AlcQuant24", "AlcQuant29", "AlcQuant35", "AlcQuant45", "CPD_TOB_age14", "CPD_TOB_age17", "CPD_TOB_age20", "CPD_TOB_age24", "CPD_TOB_age29", "CPD_TOB_age35","CPD_TOB_age45", "MJ_FREQ_14", "MJ_FREQ_17", "MJ_FREQ_20", "MJ_FREQ_24", "MJ_FREQ_29", "MJ_FREQ_35", "MJ_FREQ_45", "OC", "YC")) mtfs.w <- reshape(mtfs, v.names=c("AlcFreq14", "AlcFreq17", "AlcFreq20", "AlcFreq24", "AlcFreq29", "AlcFreq35", "AlcFreq45", "AlcQuant14", "AlcQuant17", "AlcQuant20", "AlcQuant24", "AlcQuant29", "AlcQuant35", "AlcQuant45", "CPD_TOB_age14", "CPD_TOB_age17", "CPD_TOB_age20", "CPD_TOB_age24", "CPD_TOB_age29", "CPD_TOB_age35","CPD_TOB_age45", "MJ_FREQ_14", "MJ_FREQ_17", "MJ_FREQ_20", "MJ_FREQ_24", "MJ_FREQ_29", "MJ_FREQ_35", "MJ_FREQ_45"), timevar="IDSC", idvar=c("IDYRFAM", "IDSEX", "OC", "YC"), direction="wide", sep="") ######################## ### Biometric Models ### ######################## ### ### Select variables Vars <- c("AlcFreq24", "AlcQuant24", "CPD_TOB_age24", "MJ_FREQ_24") nv <- length(Vars) ntv <- nv*2 selVars <- paste(Vars, c(rep(0,nv), rep(1,nv)), sep="") DefVars <- c("OC", "YC") ### Select Data for Analysis mtfsMzDatam <- subset(mtfs.w, ZYGOSITY==1 & IDSEX == 1, select=c(selVars, DefVars)) mtfsDzDatam <- subset(mtfs.w, ZYGOSITY==2 & IDSEX == 1, select=c(selVars, DefVars)) mtfsMzDataTm <- mtfsMzDatam[rowSums(is.na(mtfsMzDatam)) != ncol(mtfsMzDatam),] mtfsDzDataTm <- mtfsDzDatam[rowSums(is.na(mtfsDzDatam)) != ncol(mtfsDzDatam),] mtfsDataMZm <- mxData( observed=mtfsMzDataTm, type="raw" ) mtfsDataDZm <- mxData( observed=mtfsDzDataTm, type="raw" ) mtfsMzDataf <- subset(mtfs.w, ZYGOSITY==1 & IDSEX == 2, select=c(selVars, DefVars)) mtfsDzDataf <- subset(mtfs.w, ZYGOSITY==2 & IDSEX == 2, select=c(selVars, DefVars)) mtfsMzDataTf <- mtfsMzDataf[rowSums(is.na(mtfsMzDataf)) != ncol(mtfsMzDataf),] mtfsDzDataTf <- mtfsDzDataf[rowSums(is.na(mtfsDzDataf)) != ncol(mtfsDzDataf),] mtfsDataMZf <- mxData( observed=mtfsMzDataTf, type="raw" ) mtfsDataDZf <- mxData( observed=mtfsDzDataTf, type="raw" ) ### ------------------------------------------------------------------------------ ### Common Pathway ACE model ### ------------------------------------------------------------------------------ nl <- 1 # number of latent factors # Matrices ac, cc, and ec to store a, c, and e path coefficients for latent phenotype(s) Xm <- mxMatrix(type="Lower", nrow=nl, ncol=nl, free=TRUE, values=c(.6), labels=c("x11m"), name="Xm", byrow = T) Ym <- mxMatrix(type="Lower", nrow=nl, ncol=nl, free=TRUE, values=c(.6), labels=c("y11m"), name="Ym", byrow = T) Zm <- mxMatrix(type="Lower", nrow=nl, ncol=nl, free=TRUE, values=c(.6), labels=c("z11m"), name="Zm", byrow = T ) Xf <- mxMatrix(type="Lower", nrow=nl, ncol=nl, free=TRUE, values=c(.6), labels=c("x11f"), name="Xf", byrow = T) Yf <- mxMatrix(type="Lower", nrow=nl, ncol=nl, free=TRUE, values=c(.6), labels=c("y11f"), name="Yf", byrow = T) Zf <- mxMatrix(type="Lower", nrow=nl, ncol=nl, free=TRUE, values=c(.6), labels=c("z11f"), name="Zf", byrow = T ) Alm <- mxAlgebra(Xm %*% t(Xm), name="Alm") Clm <- mxAlgebra(Ym %*% t(Ym), name="Clm") Elm <- mxAlgebra(Zm %*% t(Zm), name="Elm") Alf <- mxAlgebra(Xf %*% t(Xf), name="Alf") Clf <- mxAlgebra(Yf %*% t(Yf), name="Clf") Elf <- mxAlgebra(Zf %*% t(Zf), name="Elf") A <- c(TRUE, TRUE, F, F, TRUE, TRUE, F, F, F, F, TRUE, F, F, F, F, TRUE) SpecificFree <- matrix(A, nrow=nv, ncol=nv, byrow=T) SpecificVals <- SpecificFree SpecificVals[which(SpecificVals==TRUE)] <- .4 SpecificVals[which(SpecificVals==FALSE)] <- 0 diag(SpecificVals) <- .5 labelsData <- c("AFasm", "ALCasm", NA, NA, "ALCasm", "AQasm", NA, NA, NA, NA, "CPDasm", NA, NA, NA, NA, "MJasm") labelsAsm <- matrix(labelsData, nrow=nv, ncol=nv, byrow=T) labelsCsm <- gsub("asm", "csm", labelsAsm) labelsEsm <- gsub("asm", "esm", labelsAsm) labelsAsf <- gsub("asm", "asf", labelsAsm) labelsCsf <- gsub("asm", "csf", labelsAsm) labelsEsf <- gsub("asm", "esf", labelsAsm) # Matrices as, cs, and es to store a, c, and e path coefficients for specific factors pathAsm <- mxMatrix(type="Symm", nrow=nv, ncol=nv, free=SpecificFree, values=SpecificVals, labels=labelsAsm, name="pathAsm", byrow = T ) pathCsm <- mxMatrix(type="Symm", nrow=nv, ncol=nv, free=SpecificFree, values=SpecificVals, labels=labelsCsm, name="pathCsm", byrow = T ) pathEsm <- mxMatrix(type="Symm", nrow=nv, ncol=nv, free=SpecificFree, values=SpecificVals, labels=labelsEsm, name="pathEsm", byrow = T ) pathAsf <- mxMatrix(type="Symm", nrow=nv, ncol=nv, free=SpecificFree, values=SpecificVals, labels=labelsAsf, name="pathAsf", byrow = T ) pathCsf <- mxMatrix(type="Symm", nrow=nv, ncol=nv, free=SpecificFree, values=SpecificVals, labels=labelsCsf, name="pathCsf", byrow = T ) pathEsf <- mxMatrix(type="Symm", nrow=nv, ncol=nv, free=SpecificFree, values=SpecificVals, labels=labelsEsf, name="pathEsf", byrow = T ) # Matrix and Algebra for constraint on variance of latent phenotype CovarLPm <- mxAlgebra( expression= Alm + Clm + Elm, name="CovarLPm" ) VarLPm <- mxAlgebra( expression= diag2vec(CovarLPm), name="VarLPm" ) Unitm <- mxMatrix( type="Unit", nrow=nl, ncol=1, name="Unitm") varLP1m <- mxConstraint( expression=VarLPm == Unitm, name="varLP1m") CovarLPf <- mxAlgebra( expression= Alf + Clf + Elf, name="CovarLPf" ) VarLPf <- mxAlgebra( expression= diag2vec(CovarLPf), name="VarLPf" ) Unitf <- mxMatrix( type="Unit", nrow=nl, ncol=1, name="Unitf") varLP1f <- mxConstraint( expression=VarLPf == Unitf, name="varLP1f") ### Matrix f for factor loadings on latent phenotype pathFlm <- mxMatrix(type="Full", nrow=nv, ncol=nl, free=T, values=c(.6, .6, .6, .6), labels=c("flAFm", "flAQm", "flCPDm", "flMJm"), name="pathFlm", byrow = T) pathFlf <- mxMatrix(type="Full", nrow=nv, ncol=nl, free=T, values=c(.6, .6, .6, .6), labels=c("flAFf", "flAQf", "flCPDf", "flMJf"), name="pathFlf", byrow = T) # Matrices A, C, and E compute variance components covAm <- mxAlgebra( expression=pathFlm %*% Alm %*% t(pathFlm) + pathAsm %*% t(pathAsm), name="Am" ) covCm <- mxAlgebra( expression=pathFlm %*% Clm %*% t(pathFlm) + pathCsm %*% t(pathCsm), name="Cm" ) covEm <- mxAlgebra( expression=pathFlm %*% Elm %*% t(pathFlm) + pathEsm %*% t(pathEsm), name="Em" ) covAf <- mxAlgebra( expression=pathFlf %*% Alf %*% t(pathFlf) + pathAsf %*% t(pathAsf), name="Af" ) covCf <- mxAlgebra( expression=pathFlf %*% Clf %*% t(pathFlf) + pathCsf %*% t(pathCsf), name="Cf" ) covEf <- mxAlgebra( expression=pathFlf %*% Elf %*% t(pathFlf) + pathEsf %*% t(pathEsf), name="Ef" ) covPm <- mxAlgebra( expression=Am+Cm+Em, name="Vm" ) matIm <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="matIm") invSDm <- mxAlgebra( expression=solve(sqrt(matIm*Vm)), name="iSDm") Sflm <- mxAlgebra(expression = iSDm %*% pathFlm, name = "Sflm") covPf <- mxAlgebra( expression=Af+Cf+Ef, name="Vf" ) matIf <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="matIf") invSDf <- mxAlgebra( expression=solve(sqrt(matIf*Vf)), name="iSDf") Sflf <- mxAlgebra(expression = iSDf %*% pathFlf, name = "Sflf") averageFLm <- mxAlgebra( expression = ((MZm.Sflm[1,1] + MZm.Sflm[2,1] + MZm.Sflm[3,1] + MZm.Sflm[4,1])/4), name="averageFLm") averageFLf <- mxAlgebra( expression = ((MZf.Sflf[1,1] + MZf.Sflf[2,1] + MZf.Sflf[3,1] + MZf.Sflf[4,1])/4), name="averageFLf") expMeanm <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=1, labels=c("mAFm", "mAQm", "mCPDm", "mMJm", "mAFm", "mAQm","mCPDm", "mMJm"), name="expMeanm" ) expMeanf <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=1, labels=c("mAFf", "mAQf", "mCPDf", "mMJf", "mAFf", "mAQf","mCPDf", "mMJf"), name="expMeanf" ) covMZm <- mxAlgebra( expression= rbind( cbind(Am+Cm+Em, Am+Cm), cbind(Am+Cm, Am+Cm+Em)), name="expCovMZm" ) covDZm <- mxAlgebra( expression= rbind( cbind(Am+Cm+Em, 0.5%x%Am+Cm), cbind(0.5%x%Am+Cm, Am+Cm+Em)), name="expCovDZm" ) covMZf <- mxAlgebra( expression= rbind( cbind(Af+Cf+Ef, Af+Cf), cbind(Af+Cf, Af+Cf+Ef)), name="expCovMZf" ) covDZf <- mxAlgebra( expression= rbind( cbind(Af+Cf+Ef, 0.5%x%Af+Cf), cbind(0.5%x%Af+Cf, Af+Cf+Ef)), name="expCovDZf" ) ### Combine Groups objMZm <- mxExpectationNormal( covariance="expCovMZm", means="expMeanm", dimnames=selVars ) objDZm <- mxExpectationNormal( covariance="expCovDZm", means="expMeanm", dimnames=selVars ) objMZf <- mxExpectationNormal( covariance="expCovMZf", means="expMeanf", dimnames=selVars ) objDZf <- mxExpectationNormal( covariance="expCovDZf", means="expMeanf", dimnames=selVars ) parsm <- list( Xm, Ym, Zm, Alm, Clm, Elm, CovarLPm, VarLPm, pathFlm, pathAsm, pathCsm, pathEsm, covAm, covCm, covEm, covPm, matIm, invSDm, expMeanm, Sflm, Unitm, averageFLm ) parsf <- list( Xf, Yf, Zf, Alf, Clf, Elf, CovarLPf, VarLPf, pathFlf, pathAsf, pathCsf, pathEsf, covAf, covCf, covEf, covPf, matIf, invSDf, expMeanf, Sflf, Unitf, averageFLf ) funML <- mxFitFunctionML() mtfsModelMZm <- mxModel( parsm, covMZm, mtfsDataMZm, objMZm, funML, name="MZm", mxCI(c("Alm", "Clm", "Elm"))) mtfsModelDZm <- mxModel( parsm, covDZm, mtfsDataDZm, objDZm, funML, name="DZm" ) mtfsModelMZf <- mxModel( parsf, covMZf, mtfsDataMZf, objMZf, funML, name="MZf", mxCI(c("Alf", "Clf", "Elf"))) mtfsModelDZf <- mxModel( parsf, covDZf, mtfsDataDZf, objDZf, funML, name="DZf" ) #run MN model fitML <- mxFitFunctionMultigroup(c("MZm.fitfunction","DZm.fitfunction","MZf.fitfunction","DZf.fitfunction")) Cf <- mxModel( "ComACE", parsm, parsf, varLP1m, varLP1f, mtfsDataMZm, mtfsDataDZm, mtfsDataMZf, mtfsDataDZf, mtfsModelMZm, mtfsModelDZm, mtfsModelMZf, mtfsModelDZf, fitML) Cf <- mxOption(Cf, "Standard Errors", "No") Cf <- mxOption(Cf, "Calculate Hessian", "No") #Cf <- mxOption(Cf, key="Number of Threads", value=4) CfFit <- mxRun(Cf, intervals=F) ### Run CholACE model #examine MN Output summary(CfFit) CfFit$output$status$code modelsummary <- summary(CfFit) results_mA_lower <- round(modelsummary$CIdetail[1,3], 2) results_mA_upper <- round(modelsummary$CIdetail[2,3], 2) results_mC_lower <- round(modelsummary$CIdetail[3,3], 2) results_mC_upper <- round(modelsummary$CIdetail[4,3], 2) results_mE_lower <- round(modelsummary$CIdetail[5,3], 2) results_mE_upper <- round(modelsummary$CIdetail[6,3], 2) results_fA_lower <- round(modelsummary$CIdetail[7,3], 2) results_fA_upper <- round(modelsummary$CIdetail[8,3], 2) results_fC_lower <- round(modelsummary$CIdetail[9,3], 2) results_fC_upper <- round(modelsummary$CIdetail[10,3], 2) results_fE_lower <- round(modelsummary$CIdetail[11,3], 2) results_fE_upper <- round(modelsummary$CIdetail[12,3], 2) maleA <- round(CfFit$output$algebras$MZm.Alm, 2) maleC <- round(CfFit$output$algebras$MZm.Clm, 2) maleE <- round(CfFit$output$algebras$MZm.Elm, 2) femaleA <- round(CfFit$output$algebras$MZf.Alf, 2) femaleC <- round(CfFit$output$algebras$MZf.Clf, 2) femaleE <- round(CfFit$output$algebras$MZf.Elf, 2) round(CfFit$output$matrices$MZm.pathAsm, 2) round(CfFit$output$matrices$MZm.pathCsm, 2) round(CfFit$output$matrices$MZm.pathEsm, 2) round(CfFit$output$matrices$MZf.pathAsf, 2) round(CfFit$output$matrices$MZf.pathCsf, 2) round(CfFit$output$matrices$MZf.pathEsf, 2) maleFL <- t(round(CfFit$output$algebras$MZm.Sflm, 2)) femaleFL <- t(round(CfFit$output$algebras$MZf.Sflf, 2)) avgFLmale <- round(CfFit$output$algebras$MZm.averageFLm, 2) avgFLfemale <- round(CfFit$output$algebras$MZf.averageFLf, 2) obsMeansAFM <- round(CfFit$output$matrices$MZm.expMeanm[1,1], 2) obsMeansAFF <- round(CfFit$output$matrices$MZf.expMeanf[1,1], 2) obsMeansAQM <- round(CfFit$output$matrices$MZm.expMeanm[1,2], 2) obsMeansAQF <- round(CfFit$output$matrices$MZf.expMeanf[1,2], 2) obsMeansCPDM <- round(CfFit$output$matrices$MZm.expMeanm[1,3], 2) obsMeansCPDF <- round(CfFit$output$matrices$MZf.expMeanf[1,3], 2) obsMeansMJM <- round(CfFit$output$matrices$MZm.expMeanm[1,4], 2) obsMeansMJF <- round(CfFit$output$matrices$MZf.expMeanf[1,4], 2) maleSex <- "male" femaleSex <- "female" wave <- 24 model <- "No Cohort" male <- cbind(maleSex, model, wave, results_mA_lower, maleA, results_mA_upper, results_mC_lower, maleC, results_mC_upper, results_mE_lower,maleE, results_mE_upper, maleFL, avgFLmale, obsMeansAFM, obsMeansAQM, obsMeansCPDM, obsMeansMJM) female <- cbind(femaleSex, model, wave, results_fA_lower, femaleA, results_fA_upper, results_fC_lower, femaleC, results_fC_upper, results_fE_lower,femaleE, results_fE_upper, femaleFL, avgFLfemale, obsMeansAFF, obsMeansAQF, obsMeansCPDF, obsMeansMJF) both <- rbind(male, female) colnames(both) <- c("Sex", "Model", "Wave", "Lower CI A","A", "Upper CI A", "Lower CI C","C", "Upper CI C", "Lower CI E","E", "Upper CI E", "AF FL", "AQ FL", "CPD FL", "MJ FL", "Avg FL", "AF Mean", "AQ Mean", "CPD Mean", "MJ Mean") both #take model out to run this line #write.table(both, file="/Users/zelle063/Desktop/OneAtTime_NoCohortInside.csv", sep=",", append=T, row.names = F, col.names = F) write.table(both, file="/Users/zelle063/Desktop/CohortACE.csv", sep=",", append=F, row.names = F, col.names = T)