# ----------------------------------------------------------------------------------------------------------------------
#  Script: oneSACEvcrH.R  
#  Author: Sarah Medland - Hermine Maes
#    Date: 29 02 2020 - 06 02 2026
#
# Twin Univariate ACE model to estimate causes of variation using twin data with genotypes
# Matrix style Model - Raw Continuous Data
# -------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|

# ----------------------------------------------------------------------------------------------------------------------
############ Section 1: LOAD LIBRARIES & FUNCTIONS
# ----------------------------------------------------------------------------------------------------------------------
############ Load Libraries & Options
library(OpenMx); library(mvtnorm) 
mxOption(NULL,"Default optimizer","SLSQP") #alternative optimizers: "NPSOL", "CSOLNP"
############ Create Output 
filename  <- "oneSACEvcrH"
sink(paste(filename,".Ro",sep=""), append=FALSE, split=TRUE)

# ----------------------------------------------------------------------------------------------------------------------
############ Load Functions used in this script
  fitGofs   <- function(fit) { summ <- summary(fit); cat(paste0("Mx:", fit$name,"  os=", summ$ob,"  ns=", summ$nu,"   ep=", summ$es,"   co=", sum(summ$cons),"  df=", summ$de, "  ll=", round(summ$Mi,4), "  cpu=", round(summ$cpu,4),"  opt=", summ$op,"  ver=", summ$mx,"  stc=", fit$output$status$code, "\n")) }
  fitEsts   <- function(fit) { print(round(fit$output$estimate,4)) }
  fitECIs   <- function(fit,fit2) { print(round(fit$output$estimate,4)); print(rbind(round(fit2$result,4),t(as.table(round(fit$output$confidenceIntervals,4)))[c(1,3),])) }

# ----------------------------------------------------------------------------------------------------------------------
############ Section 2: SIMULATE DATA for classical twin design analysis
# ----------------------------------------------------------------------------------------------------------------------
############ Set Seed, Parameters and Regression coefficients
set.seed(145)                                    # for reproducibility
Np        <- 1000                                # number of twin pairs per group
A         <- 0.5                                 # proportion for additive genetic variance
C         <- 0.2                                 # proportion for shared environmental variance
E         <- 0.3                                 # proportion for unique environmental variance
r0        <- 1                                   # total variance
bA        <- 0.001                               # regression on age
bS        <- 0.1                                 # regression on sex
############ Simulate MZ Twin Data (monozygotic - share 100% of segregating alleles)
piH       <- 1                                   # pihat - rel=1  
r1        <- piH*A + C                           # expected MZ correlation
age       <- runif( Np, min=20, max=70)          # simulate age 
sex1      <- rbinom( Np, size=1, prob=0.5)       # simulate sex1: 0=female, 1=male
sex2      <- sex1                                # sex2=sex1 for MZ twins
mzT1      <- numeric(Np); mzT2 <- numeric(Np)    # placeholders for phenotypes
for (i in 1:Np) { mzp <- rmvnorm(n=1,c(0,0),matrix(c(r0,r1,r1,r0),2,2)); mzT1[i] <- mzp[1]; mzT2[i] <- mzp[2] }
############ Create MZ Twin Dataset
mzDat     <- data.frame( fid=1:Np, piH, age, sex1, sex2, mzT1, mzT2, zyg=1, rel=1)
psych::describe(mzDat); cor(mzDat$mzT1,mzDat$mzT2); # plot(mzDat$mzT1,mzDat$mzT2)
mzDat$T1  <- mzDat$mzT1 + .001*mzDat$age + .1*mzDat$sex1
mzDat$T2  <- mzDat$mzT2 + .001*mzDat$age + .1*mzDat$sex2
cor(mzDat$T1,mzDat$T2)
############ Simulate DZ Twin Data (dizygotic - share ~50% of segregating alleles on average)
piH       <- rnorm( Np, mean=.5, sd=.03)         # pihat ~ N(0.5, 0.03^2) # hist(piH)
r2        <- piH*A + C                           # expected DZ correlation
age       <- runif( Np, min=20, max=70)          # simulate age 
sex1      <- rbinom( Np, size=1, prob=0.5)       # simulate sex1: 0=female, 1=male
sex2      <- rbinom( Np, size=1, prob=0.5)       # simulate sex2: 0=female, 1=male
dzT1      <- numeric(Np); dzT2 <- numeric(Np)    # placeholders for phenotypes
for (i in 1:Np) { dzp <- rmvnorm(n=1,c(0,0),matrix(c(r0,r2[i],r2[i],r0),2,2)); dzT1[i] <- dzp[1]; dzT2[i] <- dzp[2] }
############ Create DZ Twin Dataset
dzDat     <- data.frame(fid=(Np+1):(2*Np), piH, age, sex1, sex2, dzT1, dzT2, zyg=2, rel=.5)
psych::describe(dzDat); cor(dzDat$dzT1,dzDat$dzT2); # plot(dzDat$dzT1,dzDat$dzT2)
dzDat$T1  <- dzDat$dzT1 + .001*dzDat$age + .1*dzDat$sex1
dzDat$T2  <- dzDat$dzT2 + .001*dzDat$age + .1*dzDat$sex2
cor(dzDat$T1,dzDat$T2)
############ Combine MZ & DZ Data
head(mzDat); head(dzDat); dim(mzDat); dim(dzDat)
vars      <- c("fid","zyg","rel","piH","age","sex1","sex2","T1","T2")
tDat      <- rbind(mzDat[,vars],dzDat[,vars])
rMZ       <- cor(tDat[tDat$zyg==1,"T1"],tDat[tDat$zyg==1,"T2"]); rMZ
rDZ       <- cor(tDat[tDat$zyg==2,"T1"],tDat[tDat$zyg==2,"T2"]); rDZ
psych::describe(tDat)
############ Save Data to External File 
#write.table(tDat,"tDat.txt",eol="\n", quote=F)

# ----------------------------------------------------------------------------------------------------------------------
############ Section 3: PREPARE DATA for classical twin design analysis
# ----------------------------------------------------------------------------------------------------------------------
############ Read Data from External File
#tDat      <- read.table("tDat,txt", header=T); dim(tDat)
psych::describe(tDat, skew=F)
############ Select Variables for Analysis
nv        <- 1                                   # number of variables
ntv       <- nv*2                                # number of total variables
selVars   <- c('T1','T2')                        # list of variables names of observed variables
covVars   <- c('age','sex1','sex2',"rel","piH")  # list of variables names of covariates
############ Select Data for Analysis
mzData    <- subset(tDat, zyg==1, c(selVars,covVars))
dzData    <- subset(tDat, zyg==2, c(selVars,covVars))
cov(mzData[,selVars],use="complete"); cor(mzData[,selVars],use="complete")
cov(dzData[,selVars],use="complete"); cor(dzData[,selVars],use="complete")
psych::describeBy(tDat[,selVars],group=tDat$zyg,na.rm=TRUE)

# ----------------------------------------------------------------------------------------------------------------------
############ Section 4: PREPARE & RUN Saturated MODEL for 2 MZ/DZ groups
# ----------------------------------------------------------------------------------------------------------------------
############ Set Starting Values for Means and Variances
svm       <- mean(tDat[,selVars[1:nv]],na.rm=TRUE)
svv       <- var(tDat[,selVars[1:nv]],use="complete")
############ Create Data Objects for Multiple Groups including Observed Phenotypes and Covariates
dataMZ    <- mxData( observed=subset(tDat, zyg==1, c(selVars,covVars)), type="raw" )
dataDZ    <- mxData( observed=subset(tDat, zyg==2, c(selVars,covVars)), type="raw" )
############ Create Matrices for Covariates and linear Regression Coefficients
defSex    <- mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, labels=c("data.sex1","data.sex2"), name="defSex" )
defAge    <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age"), name="defAge" )
betaS     <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, labels="bS", name="betaS" )
betaA     <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, labels="bA", name="betaA" )
############ Create Algebra for expected Mean Matrices
meanMZ    <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svm, labels=c("mMZ1","mMZ2"), name="meanMZ" )
meanDZ    <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svm, labels=c("mDZ1","mDZ2"), name="meanDZ" )
emeanMZ   <- mxAlgebra( expression= meanMZ +betaS%*%defSex +betaA%*%cbind(defAge,defAge), name="emeanMZ" )
emeanDZ   <- mxAlgebra( expression= meanDZ +betaS%*%defSex +betaA%*%cbind(defAge,defAge), name="emeanDZ" )
############ Create Algebra for expected Variance/Covariance Matrices
covMZ     <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=c(svv,0,svv), labels=c("vMZ1","cMZ21","vMZ2"), name="covMZ" )
covDZ     <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=c(svv,0,svv), labels=c("vDZ1","cDZ21","vDZ2"), name="covDZ" )
############ Create Algebra for Maximum Likelihood Estimates of Twin Correlations
corMZ     <- mxAlgebra( cov2cor(covMZ), name="corMZ" )
corDZ     <- mxAlgebra( cov2cor(covDZ), name="corDZ" )
############ Create Expectation Objects for Multiple Groups linking model-implied Means and Covariances to Data
expMZ     <- mxExpectationNormal( covariance="covMZ", means="emeanMZ", dimnames=selVars )
expDZ     <- mxExpectationNormal( covariance="covDZ", means="emeanDZ", dimnames=selVars )
funML     <- mxFitFunctionML()
############ Create Model Objects for Multiple Groups
# -------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|
pars      <- list( betaS, betaA )                # parameters common to multiple groups
defs      <- list( defSex, defAge )              # definition variables (and algebras incorporating them)
modMZ     <- mxModel( pars, defs, meanMZ, emeanMZ, covMZ, corMZ, dataMZ, expMZ, funML, name="MZ" )
modDZ     <- mxModel( pars, defs, meanDZ, emeanDZ, covDZ, corDZ, dataDZ, expDZ, funML, name="DZ" )
multi     <- mxFitFunctionMultigroup( c("MZ","DZ") )
############ Create Algebra for Correlations & Confidence Intervals
stR       <- mxAlgebra( expression=cbind(corMZ[2,1],corDZ[2,1]), name="stR", dimnames=list("stR",c("rMZ","rDZ")) )  
ciR       <- mxCI( reference="stR" )                    
############ Build & Run Saturated Model [with Confidence Intervals if intervals=T] - 2 groups
modSAT    <- mxModel( "SAT", pars, covMZ, covDZ, corMZ, corDZ, modMZ, modDZ, multi, stR, ciR )
fitSAT    <- mxRun( modSAT, intervals=F )
############ Print Output & Maximum Likelihood Estimates of Twin Correlations
fitGofs(fitSAT); fitEsts(fitSAT)
print(c("erMZ"=round(fitSAT$MZ$corMZ$result[2,1],2), "erDZ"=round(fitSAT$DZ$corDZ$result[2,1],2)),quote=F)
#mxGetExpected(fitSAT, "covariance"); mxGetExpected(fitSAT, "means"); fitSAT$MZ$meanMZ; fitSAT$DZ$meanDZ

############ Section 4b: TEST data assumptions
############ Run Submodels: Constrain expected Means/Variances to be equal across Twin Order [EMO & EMVO] & Zygosity [EMVZ]
fitEMO    <- mxRun( omxSetParameters( fitSAT, label=c("mMZ1","mDZ1","mMZ2","mDZ2"), free=TRUE, values=svm, newlabels=c("mMZ","mDZ"), name="EMO" ))
fitEMVO   <- mxRun( omxSetParameters( fitEMO, label=c("vMZ1","vDZ1","vMZ2","vDZ2"), free=TRUE, values=svv, newlabels=c("vMZ","vDZ"), name="EMVO" ))
fitEMVZ   <- mxRun( omxSetParameters ( omxSetParameters( fitEMVO, label=c("mMZ","mDZ"), free=TRUE, values=svm, newlabels="mZ"), label=c("vMZ","vDZ"), free=TRUE, values=svv, newlabels='vZ', name="EMVZ" ), intervals=T )
fitGofs(fitEMO); fitEsts(fitEMO); fitGofs(fitEMVO); fitEsts(fitEMVO);  fitGofs(fitEMVZ); fitECIs(fitEMVZ, fitEMVZ$stR)
############ Print Comparative Fit Statistics
print(mxCompare( fitSAT, subs <- list(fitEMO, fitEMVO, fitEMVZ) ) )
#mxGetExpected(fitEMVZ, "covariance"); mxGetExpected(fitEMVZ, "means"); fitEMVZ$MZ$meanMZ; fitEMVZ$DZ$meanDZ

# ----------------------------------------------------------------------------------------------------------------------
############ Section 5: PREPARE & RUN ACE MODEL 1: univariate twin multigroup model + covariates 
# ----------------------------------------------------------------------------------------------------------------------
############ Set Starting Values
svp       <- svv*.1                              # start value for variance components
sve       <- svv*.8                              # start value for e variance components
############ Create Algebra for expected Mean Matrices
meanG     <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svm, labels="mean", name="meanG" )
emeanG    <- mxAlgebra( expression= meanG +betaS%*%defSex +betaA%*%cbind(defAge,defAge), name="emeanG" )
############ Create Matrices for Variance Components
covA      <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svp, label="VA11", name="VA" ) 
covC      <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svp, label="VC11", name="VC" )
covE      <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=sve, label="VE11", name="VE" )
############ 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" )
ecovMZ    <- mxAlgebra( expression= rbind( cbind(V, cMZ), cbind(t(cMZ), V)), name="ecovMZ" )
ecovDZ    <- mxAlgebra( expression= rbind( cbind(V, cDZ), cbind(t(cDZ), V)), name="ecovDZ" )
############ Create Expectation Objects for Multiple Groups
expMZ     <- mxExpectationNormal( covariance="ecovMZ", means="emeanG", dimnames=selVars )
expDZ     <- mxExpectationNormal( covariance="ecovDZ", means="emeanG", dimnames=selVars )
############ Create Model Objects for Multiple Groups
pars      <- list( meanG, betaS, betaA, covA, covC, covE, covP )
defs      <- list( defSex, defAge, emeanG )
modMZ     <- mxModel( pars, defs, covMZ, ecovMZ, dataMZ, expMZ, funML, name="MZ" )
modDZ     <- mxModel( pars, defs, covDZ, ecovDZ, dataDZ, expDZ, funML, name="DZ" )
multi     <- mxFitFunctionMultigroup( c("MZ","DZ") )
############ Create Algebra for Unstandardized and Standardized Variance Components & Confidence Intervals
stVC      <- mxAlgebra( expression=cbind(VA,VC,VE,VA/V,VC/V,VE/V,V), name="stVC", dimnames=list("stVC",c('VA','VC','VE','SA','SC','SE','V')) )  
ciVC      <- mxCI( reference= "stVC" )
############ Build & Run Model [with Confidence Intervals if intervals=T] - 2 groups
modACE    <- mxModel( "ACE1", pars, modMZ, modDZ, multi, stVC, ciVC )
fitACE    <- mxRun( modACE, intervals=TRUE )
fitGofs(fitACE); fitECIs(fitACE,fitACE$stVC)
mxCompare( fitSAT, fitACE)

############ Section 5b: TEST Significance of C, A, & A&C combined
############ Run Submodels: Fit AE Model, CE Model, & E Model
fitAE     <- mxRun( omxSetParameters( fitACE, labels="VC11", free=FALSE, values=0, name="AE" ) )
fitCE     <- mxRun( omxSetParameters( fitACE, labels="VA11", free=FALSE, values=0, name="CE" ) )
fitE      <- mxRun( omxSetParameters( fitAE, labels="VA11", free=FALSE, values=0, name="E" ) )
fitGofs(fitAE); fitEsts(fitAE);  fitGofs(fitCE); fitEsts(fitCE);  fitGofs(fitE); fitEsts(fitE)
############ Print Comparative Fit Statistics & ML Estimates of Model Parameters
print(mxCompare( fitACE, nested <- list(fitAE, fitCE, fitE) ) )
print( cbind( rbind(fitACE$name, fitAE$name, fitCE$name, fitE$name),
 round( rbind(fitACE$stVC$result,fitAE$stVC$result,fitCE$stVC$result,fitE$stVC$result),4)), quote=F)

# ----------------------------------------------------------------------------------------------------------------------
############ Section 6: PREPARE & RUN ACE MODEL 2: univariate twin model - alternate parameterization 
# ----------------------------------------------------------------------------------------------------------------------
############ Create Algebra for expected Variance/Covariance Matrices in MZ & DZ twins
relAmz    <- mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=FALSE, values=1, name="relAmz" ) 
relAdz    <- mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=FALSE, values=.5, name="relAdz" ) 
relC      <- mxMatrix( type="Unit", nrow=ntv, ncol=ntv, free=FALSE, name="relC" ) 
relE      <- mxMatrix( type="Iden", nrow=ntv, ncol=ntv, free=FALSE, name="relE" )
ecovMZ    <- mxAlgebra( expression= VA%x%relAmz + VC%x%relC + VE%x%relE, name="ecovMZ" )
ecovDZ    <- mxAlgebra( expression= VA%x%relAdz + VC%x%relC + VE%x%relE, name="ecovDZ" )
############ Create Model Objects for Multiple Groups
pars      <- list( meanG, betaS, betaA, covA, covC, covE, covP, relAmz, relAdz, relC, relE )
modMZ     <- mxModel( pars, defs, covMZ, ecovMZ, dataMZ, expMZ, funML, name="MZ" )
modDZ     <- mxModel( pars, defs, covDZ, ecovDZ, dataDZ, expDZ, funML, name="DZ" )
multi     <- mxFitFunctionMultigroup( c("MZ","DZ") )
############ Build Model [with Confidence Intervals if intervals=T]
modACE2   <- mxModel( "ACE2", pars,  modMZ, modDZ, multi, stVC, ciVC )
fitACE2   <- mxRun( modACE2, intervals=FALSE )
fitGofs(fitACE2); fitEsts(fitACE2)

# ----------------------------------------------------------------------------------------------------------------------
############ Section 7: PREPARE & RUN ACE MODEL 3: univariate twin model - expected relatedness as definition variable
# ----------------------------------------------------------------------------------------------------------------------
############ Create Data Objects for Single Group
dataTW    <- mxData( observed=subset(tDat,,c(selVars,covVars)), type="raw" )
############ Create Algebra for expected Variance/Covariance Matrices in twins
relA      <- mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=FALSE, labels=c("data.rel"), name="relA" ) 
ecovTW    <- mxAlgebra( expression= VA%x%relA + VC%x%relC + VE%x%relE, name="ecovTW" )
############ Create Expectation Objects for Single Group
expTW     <- mxExpectationNormal( covariance="ecovTW", means="emeanG", dimnames=selVars )
funML2    <- mxFitFunctionML()
############ Create Model Objects for Single Group
defs      <- list( defAge, defSex, emeanG, relA)
pars      <- list( meanG, betaS, betaA, covA, covC, covE, covP, relC, relE )
############ Build Model with Confidence Intervals
modACE3   <- mxModel( "ACE3", pars, defs, ecovTW, dataTW, expTW, funML2, stVC, ciVC )
fitACE3   <- mxRun( modACE3, intervals=TRUE )
fitGofs(fitACE3); fitEsts(fitACE3)

# ----------------------------------------------------------------------------------------------------------------------
############ Section 8: PREPARE & RUN ACE MODEL 4: univariate twin model - actual relatedness as definiation variable
# ----------------------------------------------------------------------------------------------------------------------
############ Create Algebra for expected Variance/Covariance Matrices in twins
relA      <- mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=FALSE, labels=c("data.piH"), name="relA" ) 
defs      <- list( defAge, defSex, emeanG, relA)
############ Build Model with Confidence Intervals
modACE4   <- mxModel( "ACE4", pars, defs, ecovTW, dataTW, expTW, funML2, stVC, ciVC )
fitACE4   <- mxRun( modACE4, intervals=TRUE )
fitGofs(fitACE4); fitECIs(fitACE4,fitACE4$stVC)

# ----------------------------------------------------------------------------------------------------------------------
############ Section 9: PREPARE & RUN ACE MODEL 5: univariate twin model - actual relatedness for DZs only
# ----------------------------------------------------------------------------------------------------------------------
############ Build Model with Confidence Intervals
modACE5   <- mxModel( "ACE5", pars, defs, ecovTW, dataDZ, expTW, funML, stVC, ciVC )
fitACE5   <- mxRun( modACE5, intervals=TRUE )
fitGofs(fitACE5); fitECIs(fitACE5,fitACE5$stVC)
############ Run AE Submodel
fitAE5    <- mxRun( omxSetParameters( fitACE5, labels="VC11", free=FALSE, values=0, name="AE" ) )
fitGofs(fitAE5); fitEsts(fitAE5)

# ----------------------------------------------------------------------------------------------------------------------
############ Section 10: SUMMARIZE RESULTS: univariate twin models 
# ----------------------------------------------------------------------------------------------------------------------
############ Print Goodness-of-Fit Statistics
mxCompare( fitSAT, subs <- list(fitEMO, fitEMVO, fitEMVZ, fitACE, fitAE, fitCE, fitE, fitACE2, fitACE3, fitACE4))
############ Print Parameter Estimates
print( cbind( rbind(fitACE$name, fitAE$name, fitCE$name, fitE$name, fitACE2$name, fitACE3$name, fitACE4$name, fitACE5$name),
 round( rbind(fitACE$stVC$result,fitAE$stVC$result,fitCE$stVC$result,fitE$stVC$result,
  fitACE2$stVC$result,fitACE3$stVC$result,fitACE4$stVC$result,fitACE5$stVC$result),4)), quote=F)

# ----------------------------------------------------------------------------------------------------------------------
sink()