# ----------------------------------------------------------------------------------------------------------------------
#  Script: 00_ACEvc.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)

dfBin <- df
dfBin$Twin1 <- ifelse(dfBin$Twin1 > 0, 1, 0)
dfBin$Twin2 <- ifelse(dfBin$Twin2 > 0, 1, 0)
dfBin$Sib   <- ifelse(dfBin$Sib > 0, 1, 0)

# Check the data
table(dfBin$Twin1)
table(dfBin$Twin2)
table(dfBin$Sib)

# Ensure that the data are encoded as an ordered factor
dfBin$Twin1 <- mxFactor(dfBin$Twin1, levels = 0:1)
dfBin$Twin2 <- mxFactor(dfBin$Twin2, levels = 0:1)
dfBin$Sib   <- mxFactor(dfBin$Sib  , levels = 0:1)

mzData    <- subset(dfBin, zygosity==1)
dzData    <- subset(dfBin, zygosity==2)

# Select Variables for Analysis
nv        <- 1                          # number of traits
nt        <- 2                          # number of individuals
nTH      <- 1                           # number of thresholds
ntv       <- nv*nt                      # number of total variables
selVars   <- c('Twin1', 'Twin2')        # name of traits
covVars   <- c('age1', 'age2', 'sex1', 'sex2') # list of covariates

# Set Starting Values
sVa      <- .3                       # start value for A
sVc      <- .3                       # start value for C
sVe      <- .3                       # start value for E
ThrVals  <- 0                        # start for the threshold

# ----------------------------------------------------------------------------------------------------------------------
# PREPARE MODEL

# Create Algebra for expected Mean Matrices
intercept     <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=FALSE, values=0, 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"), name="Sex" )
defAge        <- mxMatrix( type="Full", nrow=1, ncol=nt, free=FALSE, labels=c("data.age1","data.age2"), name="Age" )
expMean       <- mxAlgebra( expression = intercept + Sex%x%bS + Age%x%bA , name="expMean" )
expThr        <- mxMatrix(type="Full", nrow=nTH, ncol=ntv, free = TRUE, values = ThrVals, labels =  paste("th", 1:nTH, sep = ""), name = "expThr")

# 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" )
cons <- mxConstraint(VA+VC+VE ==1, name = "cons")      # Constraint so that the variance equals 1 (necessary for identification)

# 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),
                                           cbind(t(cMZ), V)), name="expCovMZ" )
expCovDZ  <- mxAlgebra( expression= rbind( cbind(V, cDZ),
                                           cbind(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", thresholds = "expThr", dimnames=selVars )
expDZ     <- mxExpectationNormal( covariance="expCovDZ", means="expMean", thresholds = "expThr", dimnames=selVars )
funML     <- mxFitFunctionML()

# Create Model Objects for Multiple Groups
defs      <- list( defAge, defSex )
pars      <- list( intercept, betaS, betaA, expThr, covA, covC, covE, covP )
modelMZ   <- mxModel( defs, pars, expMean, covMZ, 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( "ACEvc", modelMZ, modelDZ, multi, pars, estVC, ciACE, cons )

# ----------------------------------------------------------------------------------------------------------------------
# RUN MODEL

# Run ACE Model
fitACE    <- mxRun( modelACE, intervals=T )
sumACE    <- summary( fitACE )
sumACE

# Some different ways to extract output
fitACE$output$estimate
fitACE$VarC$result
fitACE$output$confidenceIntervals

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

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