# Simulation and Power
require(MASS)


#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$
#$#            Simulate data to test a linear Regression model                   $#$ 
#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$

# Simulating a simple linear regression model
N <- 1000

b0 <- 2
b1 <- .3
b2 <- .3

X1 <- rnorm(N)
X2 <- rnorm(N)
resid <- rnorm(N)

Y <- b0 + b1*X1 + b2*X2 + resid

# Now let's run the model

summary(lm(Y ~ X1 + X2))
summary(lm(Y ~ X1))
summary(lm(Y ~ X2))


#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$
#$#                     Simulating data using a RAM model                        $#$ 
#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$ 

# This is the algebra RAM model
# F %*% solve(I-A) %*% S %*% t(solve(I-A)) %*% t(F)

# F is a filter matrix (to remove those pesky latent variables) | Nobserved by Ntotal (Nobserved + Nlatent)
# I is an Identity matrix | Ntotal by Ntotal

# Note there are no free parameters in the F and I matrices

# A is a matrix of Asymmetrical paths (i.e. regression paths)      | Ntotal by Ntotal 
# S is a matrix of Symmetrical Paths (i.e. covariance/variance/residual paths)       | Ntotal by Ntotal

# Ok, let's simulate the data.

AA <- matrix(c(
.0, .0, .0, .0, .0, .0, .8, .0, 
.0, .0, .0, .0, .0, .0, .7, .0, 
.0, .0, .0, .0, .0, .0, .6, .0, 
.0, .0, .0, .0, .0, .0, .0, .8, 
.0, .0, .0, .0, .0, .0, .0, .7, 
.0, .0, .0, .0, .0, .0, .0, .6, 
.0, .0, .0, .0, .0, .0, .0, .0, 
.0, .0, .0, .0, .0, .0, .0, .0), 8, 8, byrow = T)

SS <- diag(diag(diag(8) - (AA %*% t(AA))))
SS[7,8] <- SS[8,7]<- .5

II <- diag(8)
FF <- cbind(diag(6), matrix(0, 6, 2))

impCov <- FF %*% solve(II-AA) %*% SS %*% t(solve(II-AA)) %*% t(FF)
impMean <- rep(0, 6)

dat <- mvrnorm(N, mu = impMean, Sigma = impCov, empirical = T)
colnames(dat) <- paste0("i", 1:6)

# Now, let's see if we can build a model to fit the data
require(OpenMx)

# Let's define a few quick parameters that we can use later
nv <- 6
ntv <- 8 
selVars <- colnames(dat) 

Afree <- matrix(F, ntv, ntv)
Afree[1:3, 7] <- T
Afree[4:6, 8] <- T

Aval <- matrix(0, ntv, ntv)
Aval[Afree==T] <- .2

Alab <- matrix(NA, ntv, ntv)
Alab[1:3, 7] <- paste0("lambda1", 1:3)
Alab[4:6, 8] <- paste0("lambda2", 1:3)

Amat <- mxMatrix( type="Full", nrow=ntv, ncol=ntv, free=Afree, values= Aval, label=Alab, name="Amat" )

Sval <- diag(ntv)
Sval[7,8] <- Sval[8,7] <- .2

Sfree <- matrix(F, ntv, ntv)
Sfree[Sval!=0] <- T
Sfree[7,7] <- Sfree[8,8] <- F

Slab <- matrix(NA, ntv, ntv)
diag(Slab) <- c(paste0("res", 1:6), "F1var", "F2var")
Slab[7,8] <- Slab[8,7] <- ("facCov")

Smat <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=Sfree, values= Sval, label=Slab, name="Smat" )

Fval <- cbind(diag(6), matrix(0, 6, 2))
Fmat <- mxMatrix( type="Full", nrow=nv, ncol=ntv, free=F, values= Fval, name="Fmat" )

Imat <- mxMatrix( type="Iden", nrow=ntv, ncol=ntv, name="Imat" )

impCov      <- mxAlgebra( expression=Fmat %*% solve(Imat - Amat) %*% Smat %*% t(solve(Imat - Amat)) %*% t(Fmat), name="impCov" )
expMean     <- mxMatrix( type="Full", nrow=1, ncol=nv, name="expMean" )

obj      <- mxExpectationNormal( covariance="impCov", means="expMean", dimnames=selVars)
fitFun   <- mxFitFunctionML()

latData  <-  mxData(dat, type="raw" )

LatModel   <- mxModel( Amat, Smat, Imat, Fmat, impCov, expMean, obj, fitFun, latData, name="2fac" )

LatFit <- mxRun(LatModel)
summary(LatFit)

LatFit$Amat
LatFit$Smat
LatFit$impCov
LatFit$expMean



#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$
#$#    Simulating gene-environment correlation in a univariate Twin Model        $#$ 
#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$ 

Nmz <- 1000
Ndz <- 1000
nv <- 1

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)

mzData <- mvrnorm(Nmz , mu = c(0,0), mzMat, empirical = T)
dzData <- mvrnorm(Ndz , mu = c(0,0), dzMat, empirical = T)

selVars <- paste("t", 1:2, sep = "")
colnames(mzData) <- colnames(dzData) <- selVars


# Putting the OpenMx model together

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

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 )

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)
summary(aceFit)



#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$
#$#                 Simulating bias univariate Twin Model                        $#$ 
#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$ 

AA <- sqrt(.3)
CC <- sqrt(.3)
EE <- sqrt(.4)

rAC <- 0
rAc <- rAC

rAE <- .2
rAe <- rAE/2



#T1   T2   A    C    E    A    C    E
Amat <- matrix(c(
0,   0,  AA,  CC,  EE,   0,   0,   0,                     # T1
0,   0,   0,   0,   0,  AA,  CC,  EE,                     # T2
0,   0,   0,   0,   0,   0,   0,   0,                     # A
0,   0,   0,   0,   0,   0,   0,   0,                     # C
0,   0,   0,   0,   0,   0,   0,   0,                     # E
0,   0,   0,   0,   0,   0,   0,   0,                     # A
0,   0,   0,   0,   0,   0,   0,   0,                     # C
0,   0,   0,   0,   0,   0,   0,   0), 8, 8, byrow = T)   # E

#T1   T2   A    C    E    A    C    E
Smz <- matrix(c(
0,   0,   0,   0,   0,   0,   0,   0,                     # T1
0,   0,   0,   0,   0,   0,   0,   0,                     # T2
0,   0,   1, rAC, rAE,   1, rAC, rAE,                     # A
0,   0, rAC,   1,   0, rAC,   1,   0,                     # C
0,   0, rAE,   0,   1, rAE,   0,   0,                     # E
0,   0,   1, rAC, rAE,   1, rAC, rAE,                     # A
0,   0, rAC,   1,   0, rAC,   1,   0,                     # C
0,   0, rAE,   0,   0, rAE,   0,   1), 8, 8, byrow = T)   # E

#T1   T2   A    C    E    A    C    E
Sdz <- matrix(c(
0,   0,   0,   0,   0,   0,   0,   0,                     # T1
0,   0,   0,   0,   0,   0,   0,   0,                     # T2
0,   0,   1, rAC, rAE,  .5, rAc, rAe,                     # A
0,   0, rAC,   1,   0, rAc,   1,   0,                     # C
0,   0, rAE,   0,   1, rAe,   0,   0,                     # E
0,   0,  .5, rAc, rAe,   1, rAC, rAE,                     # A
0,   0, rAc,   1,   0, rAC,   1,   0,                     # C
0,   0, rAe,   0,   0, rAE,   0,   1), 8, 8, byrow = T)   # E

FF <- cbind(diag(2), matrix(0, 2, 6))
II <- diag(8)

covMZ <- FF %*% solve(II - Amat) %*% Smz %*% t(solve(II - Amat)) %*% t(FF)
covDZ <- FF %*% solve(II - Amat) %*% Sdz %*% t(solve(II - Amat)) %*% t(FF)

mzData <- mvrnorm(Nmz , mu = c(0,0), covMZ, empirical = T)
dzData <- mvrnorm(Ndz , mu = c(0,0), covDZ, empirical = T)

selVars <- paste("t", 1:2, sep = "")
colnames(mzData) <- colnames(dzData) <- selVars


# Putting the OpenMx model together

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

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 )

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)
summary(aceFit)


#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$
#$#                     Power in a univariate Twin Model                         $#$ 
#$#                           Monte Carlo 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 )

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

nreps <- 10

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

for (i in 1:nreps){

mzData <- mvrnorm(Nmz , mu = c(0,0), mzMat, empirical = F)
dzData <- mvrnorm(Ndz , mu = c(0,0), dzMat, empirical = F)

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)

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)

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

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

}

# Power is the percent of the time a parameter is significant (at .05 here but you can select the p-value threshold)
table(store$pNoA<.05) # power of A
table(store$pNoC<.05) # 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 <- 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)

# 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 <- 50
pcrit <- .05

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

# 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)