# Simulation and Power
require(MASS)
require(OpenMx)


#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$
#$#                     Power in a univariate Twin Model                         $#$ 
#$#                           Monte Carlo Approach                               $#$
#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$ 

#############################################################
# You can change all these parameters
Nmz <- 200
Ndz <- 200


AA <- .3
CC <- .3
EE <- 1 - AA - CC
#############################################################

mzMat <- matrix(c(AA + CC + EE, AA + CC, AA + CC, AA + CC + EE),2)
dzMat <- matrix(c(AA + CC + EE, .5*AA + CC, .5*AA + CC, AA + CC + EE),2)

# Putting the OpenMx model together
nv <- 1
selVars <- paste("t", 1:2, sep = "")


A <-  mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=.6, label="A11", name="A" )
C <-  mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=.6, label="C11", name="C" )
E <-  mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=.6, label="E11", 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,Mean,expMean,expCovMZ,expCovDZ)


fun <- mxFitFunctionML()
mzExp <- mxExpectationNormal(covariance="expCovMZ", means="expMean", dimnames=selVars )
dzExp <- mxExpectationNormal(covariance="expCovDZ", means="expMean", dimnames=selVars )

# Now we are going to fit the model over and over and over ...

nreps <- 50

store <- as.data.frame(matrix(NA, nreps, 6))
colnames(store) <- c("A", "C", "E", "mean", "pNoA", "pNoC")

start <- Sys.time()

for (i in 1:nreps){
if (i %% 10 ==0)  {	print(i)}
mzData <- mvrnorm(Nmz , mu = c(0,0), mzMat, empirical = F)
dzData <- mvrnorm(Ndz , mu = c(0,0), dzMat, empirical = F)

colnames(mzData) <- colnames(dzData) <- selVars

MZdata  <-  mxData(mzData, type="raw" )
DZdata  <-  mxData(dzData, type="raw" )

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 <- suppressWarnings(mxRun(ace, silent = T, unsafe = T))

ceFit <- omxSetParameters(aceFit, labels = "A11", free = F, values = 0, name = "CE")
ceFit <- suppressWarnings(mxRun(ceFit, silent = T, unsafe = T))

aeFit <- omxSetParameters(aceFit, labels = "C11", free = F, values = 0, name = "AE")
aeFit <- suppressWarnings(mxRun(aeFit, silent = T, unsafe = T))

res <- mxCompare(aceFit, c(ceFit, aeFit))

store[i,] <- c(coef(aceFit), c(res[2:3, 9]) )

}

end <- Sys.time()
end - start

hist(store$A)
hist(store$C)
hist(store$E)

# Power is the percent of the time a parameter is significant (at .05 here but you can select the p-value threshold)
sum(store$pNoA<.05)/ (nreps) # power of A
sum(store$pNoC<.05)/ (nreps) # power of C


#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$
#$#                     Power in a univariate Twin Model                         $#$ 
#$#                     Non-Centrality Parameter Approach                        $#$
#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$ 


#############################################################
# You can change all these parameters
Nmz <- 100
Ndz <- 100

AA <- .3
CC <- .3
EE <- 1 - AA - CC
#############################################################

mzMat <- matrix(c(AA + CC + EE, AA + CC, AA + CC, AA + CC + EE),2)
dzMat <- matrix(c(AA + CC + EE, .5*AA + CC, .5*AA + CC, AA + CC + EE),2)

# Putting the OpenMx model together
nv <- 1

A <-  mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=.6, label="A11", name="A" )
C <-  mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=.6, label="C11", name="C" )
E <-  mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=.6, label="E11", 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,Mean,expMean,expCovMZ,expCovDZ)

fun <- mxFitFunctionML()
mzExp <- mxExpectationNormal(covariance="expCovMZ", means="expMean", dimnames=selVars )
dzExp <- mxExpectationNormal(covariance="expCovDZ", means="expMean", dimnames=selVars )

# Here are a big differences between the Monte Carlo and NCP approaches

mzData <- mvrnorm(Nmz , mu = c(0,0), mzMat, empirical = T)  # By using empirical = T your observed covariance matrix will be nearly identical 
dzData <- mvrnorm(Ndz , mu = c(0,0), dzMat, empirical = T)  # to what you provide to mvrnorm
 
selVars <- paste("t", 1:2, sep = "")
colnames(mzData) <- colnames(dzData) <- selVars

MZdata  <-  mxData(mzData, type="raw" )
DZdata  <-  mxData(dzData, type="raw" )

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, silent = T)
summary(aceFit)                                             # look at these parameters to make sure that are really close to your simulated parameters

# Fit the submodels to identify the parameters that you want power estimates for
ceFit <- omxSetParameters(aceFit, labels = "A11", free = F, values = 0, name = "CE")
ceFit <- mxRun(ceFit, silent = T)

aeFit <- omxSetParameters(aceFit, labels = "C11", free = F, values = 0, name = "AE")
aeFit <- mxRun(aeFit, silent = T)

# OpenMx gives you the a summary of estimates that we can use
res <- mxCompare(aceFit, c(ceFit, aeFit))
res
res <- as.data.frame(res)

# By dividing the difference in the Log-Likelihood by the total sample size you get an estimate of how much each family 
# contributes to the test statistic

# We can store these as Weighted NCP parameters

WncpA <- res$diffLL[2]/(Nmz+Ndz)
WncpC <- res$diffLL[3]/(Nmz+Ndz)

WncpA
WncpC

# The Wncp parameters scale linearly with sample size, making it possible to simply multiply the Wncp by the expected sample size you want to collect.

TargetSampSize <- 300
pcrit <- .05

ncpA <- WncpA * TargetSampSize
ncpC <- WncpC * TargetSampSize

# This is the NCP for A and C.  
ncpA
ncpC

# This Wncp parameter can then be added to a simple R statement to calculate power
1- pchisq(qchisq(1- pcrit, 1), 1, ncpA)
1- pchisq(qchisq(1- pcrit, 1), 1, ncpC)

# What happens if you want a pretty power plot

TargetSampSize <- 1500

powA <- 1- pchisq(qchisq(1- pcrit, 1), 1,  WncpA * 1:TargetSampSize)
powC <- 1- pchisq(qchisq(1- pcrit, 1), 1,  WncpC * 1:TargetSampSize)

plot(powA, type = "l", lwd = 3, col = "red", xlab = "Sample Size", ylab = "Power")
lines(powC, col = "blue", lwd = 3)

# Question: How many families do we need to have 80% power for both A and C

sum(powA < .8)+1
sum(powC < .8)+1