# ---------------------------------------------------------------------------------------------------------------------- # Script: 01_extrasib.R # Author: Sarah Medland # Date: 29 02 2020 # # Twin Univariate ACE model to estimate causes of variation across multiple groups # Matrix style model - Raw data - Continuous data # -------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------| # Set Working Directory # Clear Workspace rm(list=ls()) # Load Libraries & Options library(OpenMx) source("powerFun.R") # ---------------------------------------------------------------------------------------------------------------------- # PREPARE DATA # Load Data df <- read.table("TwinsSibsDataSmallN.txt", header=T) # Select Variables for Analysis nv <- 1 # number of traits nt <- 3 # number of individuals ntv <- nv*nt # number of total variables selVars <- c('Twin1', 'Twin2', 'Sib') # name of traits covVars <- c('age1', 'age2', 'age3', 'sex1', 'sex2', 'sex3') # list of covariates # Select Data for Analysis mzData <- subset(df, zygosity==1, c(selVars,covVars)) dzData <- subset(df, zygosity==2, c(selVars,covVars)) cov(mzData[,selVars],use="pairwise.complete.obs") cov(dzData[,selVars],use="pairwise.complete.obs") # Set Starting Values sMu <- 0 # start value for means sVa <- .3 # start value for A sVc <- .3 # start value for C sVe <- .3 # start value for E # ---------------------------------------------------------------------------------------------------------------------- # PREPARE MODEL # Create Algebra for expected Mean Matrices intercept <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=sMu, labels="interC", name="intercept" ) betaS <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=0, labels="betaS", name="bS" ) betaA <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=0, labels="betaA", name="bA" ) defSex <- mxMatrix( type="Full", nrow=1, ncol=nt, free=FALSE, labels=c("data.sex1","data.sex2","data.sex3"), name="Sex" ) defAge <- mxMatrix( type="Full", nrow=1, ncol=nt, free=FALSE, labels=c("data.age1","data.age2","data.age3"), name="Age" ) expMean <- mxAlgebra( expression = intercept + Sex%x%bS + Age%x%bA , name="expMean" ) # Create Matrices for Variance Components covA <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=sVa, label="VA11", name="VA" ) covC <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=sVc, label="VC11", name="VC" ) covE <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=sVe, label="VE11", name="VE" ) # Create Algebra for expected Variance/Covariance Matrices in MZ & DZ twins covP <- mxAlgebra( expression= VA+VC+VE, name="V" ) covMZ <- mxAlgebra( expression= VA+VC, name="cMZ" ) covDZ <- mxAlgebra( expression= 0.5%x%VA+ VC, name="cDZ" ) expCovMZ <- mxAlgebra( expression= rbind( cbind(V, cMZ, cDZ), cbind(t(cMZ), V, cDZ), cbind(t(cDZ), t(cDZ), V)), name="expCovMZ" ) expCovDZ <- mxAlgebra( expression= rbind( cbind(V, cDZ, cDZ), cbind(t(cDZ), V, cDZ), cbind(t(cDZ), t(cDZ), V)), name="expCovDZ" ) # 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=selVars ) expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars ) funML <- mxFitFunctionML() # Create Model Objects for Multiple Groups defs <- list(defAge, defSex) pars <- list( intercept, betaS, betaA, covA, covC, covE, covP ) modelMZ <- mxModel( defs, pars, expMean, covMZ, covDZ, expCovMZ, dataMZ, expMZ, funML, name="MZ" ) modelDZ <- mxModel( defs, pars, expMean, covDZ, expCovDZ, dataDZ, expDZ, funML, name="DZ" ) multi <- mxFitFunctionMultigroup( c("MZ","DZ") ) # Create Algebra for Variance Components rowVC <- rep('VC',nv) colVC <- rep(c('VA','VC','VE','SA','SC','SE'),each=nv) estVC <- mxAlgebra( expression=cbind(VA,VC,VE,VA/V,VC/V,VE/V), name="VarC", dimnames=list(rowVC,colVC) ) # Create Confidence Interval Objects ciACE <- mxCI( "VarC[1,4:6]" ) # Build Model with Confidence Intervals modelACE <- mxModel( "ACEsib", modelMZ, modelDZ, multi, pars, estVC, ciACE ) # ---------------------------------------------------------------------------------------------------------------------- # RUN MODEL # Run ACE Model fitACE <- mxRun( modelACE, intervals=T ) sumACE <- summary( fitACE ) sumACE fitAge <- omxSetParameters(fitACE, labels = "betaA", values = 0, free = F, name = "Age") fitAge <- mxRun( fitAge, intervals=T ) sumAge <- summary( fitAge ) sumAge mxCompare(fitACE,fitAge) # Run CE Model to obtain the -2ll test for A fitCE <- omxSetParameters(fitACE, labels = "VA11", values = 0, free = F, name = "CE") fitCE <- mxRun( fitCE, intervals=F ) # Run AE Model to obtain the -2ll test for C fitAE <- omxSetParameters(fitACE, labels = "VC11", values = 0, free = F, name = "AE") fitAE <- mxRun( fitAE, intervals=F ) # Compare the ACE model with the CE and AE models to obtain the -2LL values for A and C ModelFitTable <- as.data.frame(mxCompare(fitACE, list(fitCE, fitAE))) # Use the -2LL to calculate Wncp for A and C # First calculate the observed sample size by counting the number of MZ and DZ Families sample.size <- dim(mzData)[1] + dim(dzData)[1] # Second calculate the Wncp for A and C by dividing the -2ll by the observed sample size Wncp.A <- ModelFitTable$diffLL[!is.na(ModelFitTable$comparison)&ModelFitTable$comparison =="CE" ] / sample.size Wncp.C <- ModelFitTable$diffLL[!is.na(ModelFitTable$comparison)&ModelFitTable$comparison =="AE" ] / sample.size # Conduct the post hoc power analysis powerValue(N = sample.size, Wncp = Wncp.A) powerValue(N = sample.size, Wncp = Wncp.C)