require(OpenMx) require(MASS) mxOption(NULL, 'Default optimizer', 'NPSOL') source("powerFun.R") # Simulate Data set.seed(12345) # Setting the seed will give us all the same numbers N <- 1000 # Setting a sample size aa <- .6 # Specifying a value for the "a" path coefficient cc <- .3 # Specifying a value for the "c" path coefficient ee <- sqrt(1 - aa^2 - cc^2) # Calculating the value for "e" such that the variance is 1 AA <- aa^2 # Calculating the Genetic Variance Component CC <- cc^2 # Calculating the Common Environment Variance Component EE <- ee^2 # Calculating the Unique Environment Variance Component mzMat <- matrix(c(AA + CC + EE, AA + CC, AA + CC, AA + CC + EE),2) # Calculating the mz covariance matrix (to be simulated) dzMat <- matrix(c(AA + CC + EE, .5*AA + CC, .5*AA + CC, AA + CC + EE),2) # Calculating the dz covariance matrix (to be simulated) mzData <- mvrnorm(N , mu = c(0,0), mzMat, empirical = T) # Simulate data for MZ Twins dzData <- mvrnorm(N , mu = c(0,0), dzMat, empirical = T) # Simulate data for DZ Twins selVars <- paste("t", 1:2, sep = "") # Generate labels colnames(mzData) <- colnames(dzData) <- selVars # Label the MZ and DZ data # Fit the model using an ACE model. # Set up a general twin model nv <- 1 ntv <- nv*2 a <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6, label="a11", name="a" ) c <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6, label="c11", name="c" ) e <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6, label="e11", name="e" ) A <- mxAlgebra( expression=a %*% t(a), name="A" ) C <- mxAlgebra( expression=c %*% t(c), name="C" ) E <- mxAlgebra( expression=e %*% t(e), name="E" ) Mean <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values= 0, label="mean", name="Mean" ) expMean <- mxAlgebra( expression= cbind(Mean,Mean), name="expMean") expCovMZ <- mxAlgebra( expression= rbind ( cbind(A+C+E , A+C), cbind(A+C , A+C+E)), name="expCovMZ" ) expCovDZ <- mxAlgebra( expression= rbind ( cbind(A+C+E , 0.5%x%A+C), cbind(0.5%x%A+C , A+C+E)), name="expCovDZ" ) obs <- list(a,c,e,A,C,E,Mean,expMean,expCovMZ,expCovDZ) MZdata <- mxData(mzData, type="raw" ) DZdata <- mxData(dzData, type="raw" ) fun <- mxFitFunctionML() mzExp <- mxExpectationNormal(covariance="expCovMZ", means="expMean", dimnames=selVars ) dzExp <- mxExpectationNormal(covariance="expCovDZ", means="expMean", dimnames=selVars ) MZ <- mxModel("MZ", obs, MZdata, fun, mzExp) DZ <- mxModel("DZ", obs, DZdata, fun, dzExp) aceFun <- mxFitFunctionMultigroup(c("MZ","DZ")) ace <- mxModel("ACE", MZ, DZ, fun, aceFun) aceFit <- mxRun(ace) # Fit the model summary(aceFit) ceFit <- omxSetParameters(aceFit, labels = "a11", free = F, values = 0) # Fix the a11 path to 0 ceFit <- mxRun(ceFit) # Fit the CE model mxCompare(aceFit, ceFit) # compare the models ncpBYn <- (ceFit$output$Minus2LogLikelihood - aceFit$output$Minus2LogLikelihood)/(dim(mzData)[1] + dim(dzData)[1]) # Calculated a weighted chi squared value n<- 1:2000 # generate a sequence of sample sizes ncp <- ncpBYn * n # calculated the expected chi sq value at each n # Plot the power par( mar = c(4, 5, 2, 1) + 0.1) plot(1- pchisq(qchisq(1- .05, 1), 1, ncp), type = "l", lwd = 3, col = "red", main = " ", ylab = "Power to Detect a Significant \n Additive Genetic Variance", xlab = "Sample Size", ylim = c(0, 1), xlim = c(0, 1000)) # Calculate the sample size necessary to achieve 80% power sum(1- pchisq(qchisq(1- .05, 1), 1, ncp)<.8) + 1 # What is the Power of A the Univariate Model as C increases? modA1 <- acePow(add = .6, com = .3, Nmz = 1000, Ndz = 1000) modA2 <- acePow(add = .6, com = .5, Nmz = 1000, Ndz = 1000) modA3 <- acePow(add = .6, com = .7, Nmz = 1000, Ndz = 1000) par( mar = c(4, 5, 2, 1) + 0.1) powerPlot(maxN = 2000, Wncp = c(modA1$WncpA, modA2$WncpA, modA3$WncpA)) legText <- c("a = .6 & c = .3", "a = .6 & c = .5", "a = .6 & c = .7") legend("right", legend = legText, col = 1:3, ncol = 1, cex = 1, lwd = 3, box.lwd = 0) # What is the Power of C the Univariate Model as A increases? modC1 <- acePow(add = .3, com = .6, Nmz = 1000, Ndz = 1000) modC2 <- acePow(add = .5, com = .6, Nmz = 1000, Ndz = 1000) modC3 <- acePow(add = .7, com = .6, Nmz = 1000, Ndz = 1000) par( mar = c(4, 5, 2, 1) + 0.1) powerPlot(maxN = 2000, Wncp = c(modC1$WncpC, modC2$WncpC, modC3$WncpC)) legText <- c("a = .3 & c = .6", "a = .5 & c = .6", "a = .7 & c = .6") legend("right", legend = legText, col = 1:3, ncol = 1, cex = 1, lwd = 3, box.lwd = 0) # What is the Power of A and C in the Univariate Model if the MZ and DZ groups are of different sizes? modN1 <- acePow(add = .6, com = .6, Nmz = 5000, Ndz = 1000) modN2 <- acePow(add = .6, com = .6, Nmz = 3000, Ndz = 3000) modN3 <- acePow(add = .6, com = .6, Nmz = 1000, Ndz = 5000) par(mfrow = c(1,2), mar = c(4, 5, 2, 1) + 0.1) powerPlot(maxN = 750, Wncp = c(modN1$WncpA, modN2$WncpA, modN3$WncpA)) legText <- c("MZ:DZ = 5:1", "MZ:DZ = 1:1", "MZ:DZ = 1:5") legend("right", legend = legText, col = 1:3, ncol = 1, cex = 1, lwd = 3, box.lwd = 0) mtext("Additive Genetic Variance", cex = 1.5) powerPlot(maxN = 750, Wncp = c(modN1$WncpC, modN2$WncpC, modN3$WncpC)) mtext("Common Environmental Variance", cex = 1.5) # What is the Power of A the Univariate Model for Binomial Data (relative to a continuous variable) con <- acePow( add = .6, com = .6, Nmz = 1000, Ndz = 1000) bin5 <- acePowBiv( add = .6, com = .6, Nmz = 1000, Ndz = 1000, prev = .5) bin4 <- acePowBiv( add = .6, com = .6, Nmz = 1000, Ndz = 1000, prev = .4) bin3 <- acePowBiv( add = .6, com = .6, Nmz = 1000, Ndz = 1000, prev = .3) bin2 <- acePowBiv( add = .6, com = .6, Nmz = 1000, Ndz = 1000, prev = .2) bin1 <- acePowBiv( add = .6, com = .6, Nmz = 1000, Ndz = 1000, prev = .1) bin05 <- acePowBiv( add = .6, com = .6, Nmz = 1000, Ndz = 1000, prev = .05) powerPlot(maxN = 5000, Wncp = c(con$WncpA, bin5$WncpA, bin4$WncpA, bin3$WncpA, bin2$WncpA, bin1$WncpA, bin05$WncpA)) legText <- c("Continuous", "Binary (Prevalence = .5)","Binary (Prevalence = .4)","Binary (Prevalence = .3)", "Binary (Prevalence = .2)", "Binary (Prevalence = .1)", "Binary (Prevalence = .05)") legend("bottomright", legend = legText, col = 1:7, ncol = 1, cex = 1, lwd = 3, box.lwd = 0) # What is the Power to detect Rg in the Bivariate Model? rg1 <- bivPow(A1 = .3, A2 = .3, rg = .10, C1 = .3, C2 = .3, rc = .15, Nmz = 1000, Ndz = 1000) rg2 <- bivPow(A1 = .3, A2 = .3, rg = .30, C1 = .3, C2 = .3, rc = .30, Nmz = 1000, Ndz = 1000) rg3 <- bivPow(A1 = .3, A2 = .3, rg = .50, C1 = .3, C2 = .3, rc = .45, Nmz = 1000, Ndz = 1000) rg4 <- bivPow(A1 = .4, A2 = .4, rg = .10, C1 = .3, C2 = .3, rc = .15, Nmz = 1000, Ndz = 1000) rg5 <- bivPow(A1 = .4, A2 = .4, rg = .30, C1 = .3, C2 = .3, rc = .30, Nmz = 1000, Ndz = 1000) rg6 <- bivPow(A1 = .4, A2 = .4, rg = .50, C1 = .3, C2 = .3, rc = .45, Nmz = 1000, Ndz = 1000) rg7 <- bivPow(A1 = .5, A2 = .5, rg = .10, C1 = .3, C2 = .3, rc = .15, Nmz = 1000, Ndz = 1000) rg8 <- bivPow(A1 = .5, A2 = .5, rg = .30, C1 = .3, C2 = .3, rc = .30, Nmz = 1000, Ndz = 1000) rg9 <- bivPow(A1 = .5, A2 = .5, rg = .50, C1 = .3, C2 = .3, rc = .45, Nmz = 1000, Ndz = 1000) par(mfrow = c(3,2), mar = c(4, 5, 2, 1) + 0.1) powerPlot(maxN = 750, Wncp = c(rg1$WncpA1, rg2$WncpA1, rg3$WncpA1)) legText <- c("a = .3 & rg = .1", "a = .3 & rg = .3", "a = .3 & rg = .5") legend("right", legend = legText, col = 1:3, ncol = 1, cex = 1, lwd = 3, box.lwd = 0) mtext("Additive Genetic Variance", cex = 1) powerPlot(maxN = 2000, Wncp = c(rg1$WncpRg, rg2$WncpRg, rg3$WncpRg)) mtext("Genetic Correlation", cex = 1) powerPlot(maxN = 750, Wncp = c(rg4$WncpA1, rg5$WncpA1, rg6$WncpA1)) legText <- c("a = .4 & rg = .1", "a = .4 & rg = .3", "a = .4 & rg = .5") legend("right", legend = legText, col = 1:3, ncol = 1, cex = 1, lwd = 3, box.lwd = 0) mtext("Additive Genetic Variance", cex = 1) powerPlot(maxN = 2000, Wncp = c(rg4$WncpRg, rg5$WncpRg, rg6$WncpRg)) mtext("Genetic Correlation", cex = 1) powerPlot(maxN = 750, Wncp = c(rg7$WncpA1, rg8$WncpA1, rg9$WncpA1)) legText <- c("a = .5 & rg = .1", "a = .5 & rg = .3", "a = .5 & rg = .5") legend("right", legend = legText, col = 1:3, ncol = 1, cex = 1, lwd = 3, box.lwd = 0) mtext("Additive Genetic Variance", cex = 1) powerPlot(maxN = 2000, Wncp = c(rg7$WncpRg, rg8$WncpRg, rg9$WncpRg)) mtext("Genetic Correlation", cex = 1) # What is the Power to detect Rc in the Bivariate Model? Figure it out on your own!!! ### This takes a little while to Run Apow <- NULL Cpow <- NULL for(i in 1:99){ mod <- acePow(add = .6, com = .6, Nmz = i*100, Ndz = (100-i)*100) n <- 1:50000 Ancp <- mod$WncpA*n Ccnp <- mod$WncpC*n Apow[i] <- sum(1- pchisq(qchisq(1- .05, 1), 1, Ancp)<.8) + 1 Cpow[i] <- sum(1- pchisq(qchisq(1- .05, 1), 1, Ccnp)<.8) + 1 } par( mar = c(4, 5, 2, 1) + 0.1) plot(1:99, Apow, type = "l", lwd = 3, col = "red", main = " ", ylab = "Required Sample Size for 80% Power to Detect a \n Significant Variance Component", xlab = "MZ:DZ Ratio", ylim = c(0, 3000), xlim = c(0, 100), xaxt = "n") lines(1:99, Cpow, lwd = 3, col = "blue") axis(side = 1, at=c(1, 10, 20, 30, 40, 50, 60, 70, 80, 90,99), labels=c("1:99","10:90","20:80","30:70","40:60","50:50","60:40","70:30","80:20","90:10","99:1"), tck= -.02) legend("top", legend = c("Additive Genetic","Common Environment"), col = c("red","blue"), ncol = 1, cex = 1, lwd = 3, box.lwd = 0)