require(OpenMx)
require(psych)
source( "http://www.vipbg.vcu.edu/~vipbg/Tc24/GenEpiHelperFunctions.R")
setwd("C:/Users/user/Desktop/R")
mxOption(NULL,"Default optimizer","NPSOL")
mydata <- read.delim("FinalCholesky.dat", header=T, sep="\t", na.strings=".")
library(reshape)

mydata$caster1  <- mydata$VAR1
mydata$caster2  <- mydata$VAR1b
mydata$ogo1  <- mydata$VAR2
mydata$ogo2  <- mydata$VAR2b
mydata <- rename(mydata, c(AGEADOL="ageT1"))
mydata <- rename(mydata, c(AGEADOLb="ageT2"))
mydata <- rename(mydata, c(AGEADUL="age2T1"))
mydata <- rename(mydata, c(AGEADULb="age2T2"))
mydata <- rename(mydata, c(SEX="sexT1"))
mydata <- rename(mydata, c(SEXb="sexT2"))

mydata$ageT1[ is.na(mydata$ageT1) ] <- -999
mydata$ageT2[ is.na(mydata$ageT2) ] <- -999

mydata$age2T1[ is.na(mydata$age2T1) ] <- -999
mydata$age2T2[ is.na(mydata$age2T2) ] <- -999

mydata$sexT1[ is.na(mydata$sexT1) ] <- -999
mydata$sexT2[ is.na(mydata$sexT2) ] <- -999

mydata[mydata$ageT1==-999,]$caster1 <- NA
mydata[mydata$ageT2==-999,]$caster2 <- NA

mydata[mydata$age2T1==-999,]$ogo1 <- NA
mydata[mydata$age2T2==-999,]$ogo2 <- NA


mzData    <- subset(mydata, ZYGOSITY=="0")
dzData    <- subset(mydata, ZYGOSITY=="1")

# Select Variables for Analysis
Vars      <- c('caster', 'ogo') 
nv        <- 2       # number of variables
ntv       <- nv*2    # number of total variables
selVars   <- c("caster1", "ogo1",  "caster2", "ogo2")

# Set Starting Values
svMe      <- 5                      # start value for means
svVa      <- .4                        # start value for variance
lbVa      <- .0001                     # start value for lower bounds
svVas     <- diag(svVa,ntv,ntv)        # start values on diagonal of covariance matrix
lbVas     <- diag(lbVa,ntv,ntv)        # lower bound for variances

laMeMZ    <- c("MZ1_1","MZ1_2","MZ2_1","MZ2_2")          # labels for means for MZ twins
laMeDZ    <- c("DZ1_1","DZ1_2","DZ2_1","DZ2_2")          # labels for means for DZ twins

laVaMZ    <- c("MZVP1T1", "MZCVT1", "MZWP1", "MZCTT1T2", "MZVP2T1", "MZCTT2T1", "MZWP2", "MZVP1T2", "MZCVT2", "MZVP2T2")  # labels for (co)variances for MZ twins
laVaDZ    <- c("DZVP1T1", "DZCVT1", "DZWP1", "DZCTT1T2", "DZVP2T1", "DZCTT2T1", "DZWP2", "DZVP1T2", "DZCVT2", "DZVP2T2") # labels for (co)variances for DZ twins

# PREPARE MODEL
 
BAGE <- mxMatrix( type="Full", nrow=1, ncol=2, free=TRUE, values= .6, label=c("bage_caster", "bage_ogo"), name="betaage")
BSEX <- mxMatrix( type="Full", nrow=1, ncol=2, free=TRUE, values= .6, label=c("bsex_caster", "bsex_ogo"), name="betasex")
BAGE2 <- mxMatrix( type="Full", nrow=1, ncol=2, free=TRUE, values= .6, label=c("bage2_caster", "bage2_ogo"), name="betaage2")

OBSAGE <- mxMatrix(  type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.ageT1", "data.ageT2"), name="Age"  )
OBSAGE2 <- mxMatrix(  type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.age2T1", "data.age2T2"), name="Age2"  )
OBSEX <- mxMatrix(  type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.sexT1", "data.sexT2"), name="Sex"  )

meanMZ    <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=.6, labels=laMeMZ, name="meanMZ" )
meanDZ    <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=.6, labels=laMeDZ, name="meanDZ" )

expMeanMZ <- mxAlgebra( expression= meanMZ + betasex %x%Sex + betaage %x%Age + betaage2 %x%Age2, name="expMeanMZ" )
expMeanDZ <- mxAlgebra( expression= meanDZ + betasex %x%Sex + betaage %x%Age + betaage2 %x%Age2, name="expMeanDZ" )

# Algebra for expected Variance/Covariance Matrices in MZ & DZ twins
expCovMZ  <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=svVas, labels=laVaMZ, name="expCovMZ" )
expCovDZ  <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=svVas, labels=laVaDZ, name="expCovDZ" )

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

inverseMZ <- mxAlgebra(solve(sqrt(I * expCovMZ)), name = "inverseMZ")
expCorMZ <- mxAlgebra(inverseMZ %*% expCovMZ %*% inverseMZ, name = "expCorMZ")

inverseDZ <- mxAlgebra(solve(sqrt(I * expCovDZ)), name = "inverseDZ")
expCorDZ <- mxAlgebra(inverseDZ %*% expCovDZ %*% inverseDZ, name = "expCorDZ")

# Data objects for Multiple Groups
dataMZ    <- mxData( observed=mzData, type="raw" )
dataDZ    <- mxData( observed=dzData, type="raw" )

# Objective objects for Multiple Groups
expMZ     <- mxExpectationNormal( covariance="expCovMZ", means="expMeanMZ", dimnames=selVars )
expDZ     <- mxExpectationNormal( covariance="expCovDZ", means="expMeanDZ", dimnames=selVars )
funML     <- mxFitFunctionML()

# Combine Groups
modelMZ   <- mxModel( "MZ", BAGE, BAGE2, BSEX, OBSAGE, OBSAGE2, OBSEX, meanMZ, expMeanMZ, expCovMZ, expCorMZ, dataMZ, inverseMZ, I, expMZ, funML )
modelDZ   <- mxModel( "DZ", BAGE, BAGE2, BSEX, OBSAGE, OBSAGE2, OBSEX, meanDZ, expMeanDZ, expCovDZ, expCorDZ, dataDZ, inverseDZ, I, expDZ, funML )

multi     <- mxFitFunctionMultigroup( c("MZ","DZ") )
ciCov     <- mxCI( c('MZ.expCovMZ', 'DZ.expCovDZ'))
ciMean    <- mxCI( c('MZ.expMeanMZ','DZ.expMeanDZ')) 
ciCor <- mxCI(c('MZ.expCorMZ','DZ.expCorDZ'))

twinSatModel   <- mxModel( "twinSat", modelMZ, modelDZ, multi, ciCov, ciMean, ciCor)

# Run Saturated Model
twinSatFit     <- mxRun( twinSatModel, intervals=F )

twinSat2Model   <- omxSetParameters( twinSatModel, labels=c("bage_ogo", "bage2_caster"), free=FALSE, values=0 )

# RUN SUBMODELS

# Constrain expected means to be equal across twin order

eqMeansTwinModel    <- mxModel(twinSat2Model, name="eqMeansTwin" )
eqMeansTwinModel    <- omxSetParameters( eqMeansTwinModel, label=c("MZ1_1","MZ2_1"), free=TRUE, values=svMe, newlabels='MZ1' )
eqMeansTwinModel    <- omxSetParameters( eqMeansTwinModel, label=c("MZ1_2","MZ2_2"), free=TRUE, values=svMe, newlabels='MZ2' )
eqMeansTwinModel    <- omxSetParameters( eqMeansTwinModel, label=c("DZ1_1","DZ2_1"), free=TRUE, values=svMe, newlabels='DZ1' )
eqMeansTwinModel    <- omxSetParameters( eqMeansTwinModel, label=c("DZ1_2","DZ2_2"), free=TRUE, values=svMe, newlabels='DZ2' )

eqMeansTwinFit      <- mxRun( eqMeansTwinModel, intervals=F )

mxCompare(twinSat3Fit, eqMeansTwinFit2)

# Constrain expected Covariances, and Variances to be equal across twin order 
eqMVarsTwinModel    <- mxModel(eqMeansTwinModel, name="eqMVarsTwin" )
eqMVarsTwinModel    <- omxSetParameters( eqMVarsTwinModel, label=c("MZVP1T1","MZVP1T2"), free=TRUE, values=2, newlabels='vMZP1' )
eqMVarsTwinModel    <- omxSetParameters( eqMVarsTwinModel, label=c("MZVP2T1","MZVP2T2"), free=TRUE, values=2, newlabels='vMZP2' )
eqMVarsTwinModel    <- omxSetParameters( eqMVarsTwinModel, label=c("MZCVT1","MZCVT2"), free=TRUE, values=.5, newlabels='CVMZ' )

eqMVarsTwinModel    <- omxSetParameters( eqMVarsTwinModel, label=c("DZVP1T1","DZVP1T2"), free=TRUE, values=2, newlabels='vDZP1' )
eqMVarsTwinModel    <- omxSetParameters( eqMVarsTwinModel, label=c("DZVP2T1","DZVP2T2"), free=TRUE, values=2, newlabels='vDZP2' )
eqMVarsTwinModel    <- omxSetParameters( eqMVarsTwinModel, label=c("DZCVT1","DZCVT2"), free=TRUE, values=.5, newlabels='CVDZ' )

eqMVarsTwinFit      <- mxRun( eqMVarsTwinModel, intervals=F )

subs <- list(eqMeansTwinFit2, eqMVarsTwinFit)
mxCompare(twinSat3Fit, subs)

# Constrain expected Means, Covariances, and Variances to be equal across twin order and zygosity
eqMVarsZygModel     <- mxModel(eqMVarsTwinModel, name="eqMVarsZyg" )
eqMVarsZygModel     <- omxSetParameters( eqMVarsZygModel, label=c("MZ1","DZ1"), free=TRUE, values=svMe, newlabels='mZ1' )
eqMVarsZygModel     <- omxSetParameters( eqMVarsZygModel, label=c("MZ2","DZ2"), free=TRUE, values=svMe, newlabels='mZ2' )
eqMVarsZygModel     <- omxSetParameters( eqMVarsZygModel, label=c("vMZP1","vDZP1"), free=TRUE, values=svMe, newlabels='vZ1' )
eqMVarsZygModel     <- omxSetParameters( eqMVarsZygModel, label=c("vMZP2","vDZP2"), free=TRUE, values=svMe, newlabels='vZ2' )
eqMVarsZygModel     <- omxSetParameters( eqMVarsZygModel, label=c("CVMZ","CVDZ"), free=TRUE, values=svMe, newlabels='CV' )

eqMVarsZygFit       <- mxRun( eqMVarsZygModel, intervals=F )

# Constrain CTCT to be equal 
eqMVarsCTCTZygModel     <- mxModel(eqMVarsCTCTTwinModel, name="eqMCTCTVarsZyg" )
eqMVarsCTCTZygModel     <- omxSetParameters( eqMVarsCTCTZygModel, label=c("MZCTT2T1"," MZCTT1T2"), free=TRUE, values=svMe, newlabels='MZCT' )
eqMVarsCTCTZygModel     <- omxSetParameters( eqMVarsCTCTZygModel, label=c("DZCTT2T1"," DZCTT1T2"), free=TRUE, values=svMe, newlabels='DZCT' )

eqMVarsZygCTCTFit       <- mxRun( eqMVarsCTCTZygModel, intervals=F )
summary(eqMVarsZygCTCTFit)

subs <- list(eqMeansTwinFit, eqMVarsTwinFit, eqMVarsZygFit, eqMVarsZygCTCTFit)
mxCompare(twinSatFit, subs)

eqMVarsZygFit       <- mxRun( eqMVarsZygModel, intervals=T )
summary(eqMVarsZygFit)