# ############# ----------------------------------------------
# #############  PART 1  NB: PAGE NUMBERS are approximate
# ############# ----------------------------------------------
#
# -> page 4 
# ->
rm(list=ls(all=TRUE))   # clear the memory
library(psych)   
library(OpenMx)
#mxOption(NULL, "Default optimizer","NPSOL") 
mxOption(NULL, "Default optimizer","CSOLNP")
skinfold=read.table(file='skintwindata')
colnames(skinfold)
# <-
# <- page 4
# ############# ----------------------------------------------
# -> page 5 
# mz dz data frames ->
selvars=c(
"bic_T1","ssc_T1","sil_T1","tri_T1",
"bic_T2","ssc_T2","sil_T2","tri_T2", 
"sex_T1","sex_T2","zyg2")
# create the MZ and DZ data frames
mzskinfold=skinfold[skinfold$zyg2==1,selvars]
dzskinfold=skinfold[skinfold$zyg2==2,selvars]
# use the psych function describe () to get descriptives
describe(mzskinfold, skew=F, range=F)
describe(dzskinfold, skew=F, range=F)
# <-
# <- page 5
# ############# ----------------------------------------------
# -> page 7
# ->
selvars = c("bic_T1","tri_T1","ssc_T1","sil_T1","bic_T2","tri_T2","ssc_T2","sil_T2")  # the array with the variable names
Smz4=cov(mzskinfold[,selvars], use='complete')
Rmz4=cor(mzskinfold[,selvars], use='complete')
Sdz4=cov(dzskinfold[,selvars], use='complete')
Rdz4=cor(dzskinfold[,selvars], use='complete')
round(Smz4,3)
round(Rmz4,3)
round(Sdz4,3)
round(Rdz4,3)
# <-
# <- page 7
# ############# ----------------------------------------------
# -> page 7
# express covariance matrix using stdevs and correlation matrix
Vmz=(diag(Smz4))  # variances from diagonal to vector Vmz
Sdmz=diag(sqrt(Vmz))  # standard deviation into diagonal
Smz = Sdmz%*%Rmz4%*%t(Sdmz)  #### IMPORTANT TO KNOW
round(Smz,3)
round(Smz4-Smz,3) # compare with previous calculation
# <-
# <- page 7
# ############# ----------------------------------------------
# -> page 8
# ->
Rmz=cov2cor(Smz4)
round(Rmz,3)
#
Sdmzi=diag(sqrt(1/Vmz))   # Sdmzi diagonal with 1/std
Rmz=Sdmzi%*%Smz4%*%t(Sdmzi)
round(Rmz,3)
round(Rmz4-Rmz,3) # compare with previous calculation
# <-
# <- page 8
# ############# ----------------------------------------------
# -> page 8
# ->
# regression analyses in twin 1 members (MZ and DZ)
r11=summary(lm(bic_T1~sex_T1, data=skinfold))
r12=summary(lm(tri_T1~sex_T1, data=skinfold))
r13=summary(lm(ssc_T1~sex_T1, data=skinfold))
r14=summary(lm(sil_T1~sex_T1, data=skinfold))
# regression analyses in twin 2 members (MZ and DZ)
r21=summary(lm(bic_T2~sex_T2, data=skinfold))
r22=summary(lm(tri_T2~sex_T2, data=skinfold))
r23=summary(lm(ssc_T2~sex_T2, data=skinfold))
r24=summary(lm(sil_T2~sex_T2, data=skinfold))
# see the R2 – proportion of explained variance
print(c(r11$r.squared, r12$r.squared, r13$r.squared, r14$r.squared))
print(c(r21$r.squared, r22$r.squared, r23$r.squared, r24$r.squared))
# <-
# <- page 8
# ############# ----------------------------------------------
# page 9 ->
# ->
# --------------------------------- SAT model with sex as covariate
#
selVars = c("bic_T1","tri_T1","ssc_T1","sil_T1","bic_T2","tri_T2","ssc_T2","sil_T2")  # the array with the variable names
#
nv=4  # number of phenotypes
ntv=2*nv  # nv X 2 ... 4 phenotypes observed in twin pairs
svRmz=Rmz4  # starting values correlation matrix MZ
svSdmz=diag(sqrt(diag(Smz4))) ## starting values stdevs MZ
#
svRdz=Rdz4                    # starting values correlation matrix DZ
svSddz=diag(sqrt(diag(Sdz4))) ## starting values stdevs DZ
#
svb0=c(18,21,20,20)    #  intercepts b0 in skinfold=b0+b1*sex + residual
svb1=c(.01,.01,.01,.01)    # slope   b1 in skinfold=b0+b1*sex + residual
# Create regression model for expected Mean Matrices
B0_     <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=svb0, labels=c("b0b","b0t","b0s","b0si"), name="b0sex" )
B1_     <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=svb1, labels=c("b1b","b1t","b1s","b1si"), name="b1sex" )
Sex1 <- mxMatrix(type="Full", nrow=1, ncol=1, free=FALSE, labels=c('data.sex_T1'),name="sex1")  # sex twin 1
Sex2 <- mxMatrix(type="Full", nrow=1, ncol=1, free=FALSE, labels=c('data.sex_T2'),name="sex2")  # sex twin 2
#
ExpMean <- mxAlgebra(expression=cbind(b0sex+b1sex*sex1,b0sex+b1sex*sex2), name='expMean')
# 
# Model covariance matrix as function of correlation matrix (RMZ RDZ) and stdevs 
corMZ      <- mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=TRUE, values=svRmz, 
                        name="RMZ" ) # 
corDZ      <- mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=TRUE, values=svRdz, 
                        name="RDZ" ) #
sdMZ      <- mxMatrix( type="Diag", nrow=ntv, ncol=ntv, free=TRUE, values=svSdmz, 
                       name="sdMZ" )
sdDZ      <- mxMatrix( type="Diag", nrow=ntv, ncol=ntv, free=TRUE, values=svSddz, 
                       name="sdDZ" )
# 
ExpCovMZ     <- mxAlgebra( expression=sdMZ%*%RMZ%*%t(sdMZ), name="expCovMZ" )
ExpCovDZ     <- mxAlgebra( expression=sdDZ%*%RDZ%*%t(sdDZ), name="expCovDZ" )
#
dataMZ    <- mxData( observed=mzskinfold, type="raw" )
dataDZ    <- mxData( observed=dzskinfold, type="raw" )
#
expMZ     <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars )
expDZ     <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars )
funML     <- mxFitFunctionML()
# model
modelMZ   <- mxModel(B0_, B1_, corMZ, sdMZ, ExpMean,  Sex1, Sex2, ExpCovMZ, ExpMean, dataMZ, expMZ, funML, name="MZ" )
modelDZ   <- mxModel(B0_, B1_, corDZ, sdDZ, ExpMean,  Sex1, Sex2, ExpCovDZ, ExpMean, dataDZ, expDZ, funML, name="DZ" )	
multi     <- mxFitFunctionMultigroup( c("MZ","DZ") )
#
modelSAT4  <- mxModel("SAT", modelMZ, modelDZ, multi)
# Run SAT Model
#fitSAT4    <- mxRun( modelSAT4)
fitSAT4    <- mxTryHard( modelSAT4,20)
sumSAT4    <- summary( fitSAT4 )
#
round(fitSAT4$MZ$RMZ$values,3)
round(fitSAT4$DZ$RDZ$values,3)
round(fitSAT4$MZ$b1sex$values,3)
round(fitSAT4$MZ$b0sex$values,3)
# ------------------------------end  SAT model with sex as covariate
# <-
# <- page 9
# ############# ----------------------------------------------
#-> page 12
# ->
F_ADE = function(rmz, rdz) { c(4*rdz-rmz, 2*rmz-4*rdz) }
rmz=c(.770, .796, .769, .808)
rdz=c(.309, .267, .293, .297)
print(F_ADE(rmz,rdz))
# <-
#<- page 12
# ############# ----------------------------------------------
# -> page 14/15
# ->
## ADE model with sex as covariate
# ----------------------------------------------------------------------
selVars = c("bic_T1","tri_T1","ssc_T1","sil_T1","bic_T2","tri_T2","ssc_T2","sil_T2")  # the array with the variable names
# PREPARE MODEL
#
nv=4  # number of variable
ntv=2*nv  # number of variables times 2 (twins)
# starting values based on the phenotypic cov matrix
svS=matrix(c(  
  14,     12,     15,     18,
  12,     12,     14,     17,
  15,     14,     23,     25,
  18,     17,     25,     31),4,4)
#
svVA=svS*.4  # rough guess of A cov matrix
svVE=svS*.2  # rough guess of E cov matrix
svVD=svS*.4  # rough guess of D cov matrix  
rvA=round(cov2cor(svVA),2) # rough guess of A correlation matrix
rvD=round(cov2cor(svVE),2) # rough guess of E correlation matrix
rvE=round(cov2cor(svVD),2) # rough guess of D correlation matrix
svsA=sqrt(diag(svVA)) # A standard deviations 
svsD=sqrt(diag(svVD)) # D standard deviations
svsE=sqrt(diag(svVE)) # E standard deviations
#
svb0=c(18,21,20,20)   
svb1=c(.0,.0,.0,.0)
#
# Create regression model for expected Mean Matrices
B0_     <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=svb0, labels=c("b0b","b0t","b0s","b0si"), name="b0sex" )
B1_     <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=svb1, labels=c("b1b","b1t","b1s","b1si"), name="b1sex" )
Sex1 <- mxMatrix(type="Full", nrow=1, ncol=1, free=FALSE, labels=c('data.sex_T1'),name="sex1")  # sex twin 1
Sex2 <- mxMatrix(type="Full", nrow=1, ncol=1, free=FALSE, labels=c('data.sex_T2'),name="sex2")  # sex twin 2
#
ExpMean <- mxAlgebra(expression=cbind(b0sex+b1sex*sex1,b0sex+b1sex*sex2), name='expMean')
# 
# Create Matrices for Variance Components
corA      <- mxMatrix( type="Stand", nrow=nv, ncol=nv, free=TRUE, values=rvA, label=c("rA21","rA31","rA41","rA32","rA42","rA43"), name="RA" )
corD      <- mxMatrix( type="Stand", nrow=nv, ncol=nv, free=TRUE, values=rvD, label=c("rD21","rD31","rD41","rD32","rD42","rD43"), name="RD" )
corE      <- mxMatrix( type="Stand", nrow=nv, ncol=nv, free=TRUE, values=rvE, label=c("rE21","rE31","rE41","rE32","rE42","rE43"), name="RE" )
#
sdA      <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=svsA, label=c("sA1","sA2","sA3","sA4"), name="sA" )
sdD      <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=svsD, label=c("sD1","sD2","sD3","sD4"), name="sD" )
sdE      <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=svsE, label=c("sE1","sE2","sE3","sE4"), name="sE" )
#
# Create Algebra for expected Variance/Covariance Matrices in MZ & DZ twins
#
VA      <- mxAlgebra( expression= sA%*%RA%*%t(sA), name="VA" )
VD      <- mxAlgebra( expression= sD%*%RD%*%t(sD), name="VD" )
VE      <- mxAlgebra( expression= sE%*%RE%*%t(sE), name="VE" )
covP      <- mxAlgebra( expression= VA+VD+VE, name="V" )
covMZ     <- mxAlgebra( expression= VA+VD, name="cMZ" )
covDZ     <- mxAlgebra( expression= 0.5%x%VA+.25%x%VD, 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=mzskinfold, type="raw" )
dataDZ    <- mxData( observed=dzskinfold, type="raw" )
#
# Create Expectation Objects for Multiple Groups
expMZ     <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars )
expDZ     <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars )
funML     <- mxFitFunctionML()
#
pars      <- list(B0_, B1_, corA, corD, corE, sdA, sdD, sdE, VA, VD, VE, covP )
modelMZ   <- mxModel( pars, ExpMean,  Sex1, Sex2, covMZ, expCovMZ, ExpMean, dataMZ, expMZ, funML, name="MZ" )
modelDZ   <- mxModel( pars, ExpMean,  Sex1, Sex2, covDZ, expCovDZ, ExpMean, dataDZ, expDZ, funML, name="DZ" )	
multi     <- mxFitFunctionMultigroup( c("MZ","DZ") )
#
modelADE4  <- mxModel("ACE4", pars, modelMZ, modelDZ, multi)
# Run Model
fitADE4    <- mxRun( modelADE4 )
sumADE4    <- summary( fitADE4 )
mxCompare(fitSAT4, fitADE4)
# <- page 14 / page 15
# ############# ----------------------------------------------
# page 15 ->
VA_ = fitADE4$VA$result
RA_ = fitADE4$RA$values
sA_ = fitADE4$sA$values
# <- page 15
#
# -> page 16
print(VA_)
sA_%*%RA_%*%sA_
# <- page 16
# ############# ----------------------------------------------
# -> page 16
VA_ = fitADE4$VA$result[1:2,1:2]
VD_ = fitADE4$VD$result[1:2,1:2]
VE_ = fitADE4$VE$result[1:2,1:2]
# <- page 16
# ############# ----------------------------------------------
# -> page 17
VA_ = fitADE4$VA$result[1:2,1:2]
RA_ = fitADE4$RA$values[1:2,1:2]
sA_ = fitADE4$sA$values[1:2,1:2]
# <- page 17
# -> page 18
VPh_ = VA_ + VD_ + VE_
diag(VA_ / VPh_) # the standardized A variance component
diag(VD_ / VPh_) # the standardized D variance component
diag(VE_ / VPh_) # the standardized E variance component
# <- page 18
# ############# ----------------------------------------------
# -> page 20
cov2cor(VA_)  # correlations
sqrt(diag(VA_)) # st devs
cov2cor(VD_)  # correlations
sqrt(diag(VD_)) # st devs
cov2cor(VE_)  # correlations
sqrt(diag(VE_)) # st devs
# <- page 20
# ############# ----------------------------------------------
# -> page 18
VPh_ = VA_ + VD_ + VE_
diag(VA_ / VPh_) # the standardized A variance component
diag(VD_ / VPh_) # the standardized D variance component
diag(VE_ / VPh_) # the standardized E variance component
# <- page 18
# ############# ----------------------------------------------
# -> page 19
RD_ = fitADE4$RD$values
# <- oage 19
# ############# ---------------------------------------------- 
# -> page 19
RD_ = fitADE4$RD$values
sD_ = fitADE4$sD$values
sD_%*%RD_%*%sD_
# <- page 19
# ############# ----------------------------------------------
# ->
# -> page 20 
tmp=omxSetParameters(fitADE4, labels=c("rD21","rD31","rD41","rD32","rD42","rD43"), values=1.0, free=FALSE)
fittmp=mxRun(tmp)
mxCompare(fitADE4, fittmp)
# <- page 20 
#
# ############# ----------------------------------------------
# ############# END PART 1
# ############# ----------------------------------------------
#
# -> page 21
# ->
rm(list=ls(all=TRUE))   # clear the memory
library(psych)   
library(OpenMx)
#mxOption(NULL, "Default optimizer","NPSOL")
#mxOption(NULL, "Default optimizer","CSOLNP")
N4=read.table(file='ntwindata')
colnames(N4)
#
N4mz=subset(N4,zyg==1, select =c("N1_4","N1_6","N1_7","N1_10","N2_4","N2_6","N2_7","N2_10"))
N4dz=subset(N4,zyg==2, select =c("N1_4","N1_6","N1_7","N1_10","N2_4","N2_6","N2_7","N2_10"))
#
describe(N4mz, skew=F)
describe(N4dz, skew=F)
# <-
# <- page 21
## ############# ----------------------------------------------
# -> page 22
rmzs=round(cor(N4mz),3)
rdzs=round(cor(N4dz),3)
print(rmzs)
print(rdzs)
# <- page 22
# ############# ----------------------------------------------
# -> page 23
# ->
selVars = c("N1_4","N1_6","N1_7","N1_10","N2_4","N2_6","N2_7","N2_10")  # the array with the variable names
#
# unconstrained .... get good values
rmzs=round( cor(N4mz),2)   # correlations
rdzs=round( cor(N4dz),2)
sdmz=round( sqrt(diag(cov(N4mz))),2)
sddz=round( sqrt(diag(cov(N4dz))),2)
memz=apply(N4mz[,selVars[1:4]],2,mean)+apply(N4mz[,selVars[5:8]],2,mean)
memz=round( memz/2,2)
medz=apply(N4dz[,selVars[1:4]],2,mean)+apply(N4dz[,selVars[5:8]],2,mean)
medz=round( medz/2,2)
# 
nv=4
ntv=2*nv
#
# Create regression model for expected Mean Matrices
#
B0_mz <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=memz, labels=c("b01mz","b02mz","b03mz","b04mz"), name="b0MZ" )
ExpMean_mz <- mxAlgebra(expression=cbind(b0MZ,b0MZ), name='expMeanMZ')
corPh_mz      <- mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=TRUE, values=rmzs, name="RPhMZ" )
sdPh_mz      <- mxMatrix( type="Diag", nrow=ntv, ncol=ntv, free=TRUE, values=sdmz, name="sPhMZ" )
evalPh_mz <- mxAlgebra(eigenval(RPhMZ), name='eval_PhMZ')
ExpCov_mz      <- mxAlgebra( expression= sPhMZ%*%RPhMZ%*%t(sPhMZ), name="expCovMZ" )
#
#
B0_dz <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=memz, labels=c("b01dz","b02dz","b03dz","b04dz"), name="b0DZ" )
ExpMean_dz <- mxAlgebra(expression=cbind(b0DZ,b0DZ), name='expMeanDZ')
corPh_dz      <- mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=TRUE, values=rdzs, name="RPhDZ" )
sdPh_dz      <- mxMatrix( type="Diag", nrow=ntv, ncol=ntv, free=TRUE, values=sddz, name="sPhDZ" )
evalPh_dz <- mxAlgebra(eigenval(RPhDZ), name='eval_PhDZ')
ExpCov_dz      <- mxAlgebra( expression= sPhDZ%*%RPhDZ%*%t(sPhDZ), name="expCovDZ" )
#
#
dataMZ    <- mxData( observed=N4mz, type="raw" )
dataDZ    <- mxData( observed=N4dz, type="raw" )
#
# Create Expectation Objects for Multiple Groups
expMZ     <- mxExpectationNormal( covariance="expCovMZ", means="expMeanMZ", dimnames=selVars )
expDZ     <- mxExpectationNormal( covariance="expCovDZ", means="expMeanDZ", dimnames=selVars )
funML     <- mxFitFunctionML()
#
modelMZ   <- mxModel(dataMZ, B0_mz, ExpMean_mz, corPh_mz, sdPh_mz, ExpCov_mz, evalPh_mz, expMZ, funML, name="MZ" )
modelDZ   <- mxModel(dataDZ, B0_dz, ExpMean_dz, corPh_dz, sdPh_dz, ExpCov_dz, evalPh_dz, expDZ, funML, name="DZ" )	
multi     <- mxFitFunctionMultigroup( c("MZ","DZ") )
#
modelSAT4A  <- mxModel("Sat", modelMZ, modelDZ, multi)
# ---------------------------------------------------------------------------
# RUN sat MODEL
#
fitSAT4A    <- mxRun( modelSAT4A )
sfitSAT4A = summary( fitSAT4A) 
#fitSAT4A    <- mxTryHard( modelSAT4A, 20)
#modelSAT4B <- omxSetParameters(fitSAT4A,labels=c("b01mz","b02mz","b03mz","b04mz","b01dz","b02dz","b03dz","b04dz"), 
#newlabels=c("b01","b02","b03","b04","b01","b02","b03","b04"), free=T, values=c(2))
#fitSAT4B    <- mxTryHard( modelSAT4B, 20)
#fitSAT4B    <- mxTryHard( modelSAT4B, 20)
#mxCompare(fitSAT4A, fitSAT4B)
# <-
# <- page 23
# ############# ----------------------------------------------
# -> page 24
round(fitSAT4A$MZ$RPhMZ$values,3)
round(fitSAT4A$DZ$RPhDZ$values,3)
# <- page 24
# ############# ----------------------------------------------
# -> page 24
# impose equality contraints on the means
modelSAT4B <- omxSetParameters(fitSAT4A,
     labels=c("b01mz","b02mz","b03mz","b04mz","b01dz","b02dz","b03dz","b04dz"),                                
  newlabels=c("b01","b02","b03","b04","b01","b02","b03","b04"), 
       free=T, values=c(2))
#
fitSAT4B    <- mxRun( modelSAT4B)
#fitSAT4B    <- mxTryHard( modelSAT4B, 20)
mxCompare(fitSAT4A, fitSAT4B)
# <- page 24
# ############# ----------------------------------------------
# -> page 25 / 26
# ->
# ADE model  version 2 ---------------------------------------
nv=4
ntv=2*nv
svS=matrix(.5,4,4); diag(svS)=1  # approx S
svVA=svS*.4  # starting value  SA
svVE=svS*.5  # starting values SE
svVD=svS*.1   #  starting values SD
#
svb0=c(rep(2.5,4))   # means
#
# Create regression model for expected Mean Matrices
#
B0_     <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=svb0, labels=c("b01","b02","b03","b04"), name="b0" )
#
ExpMean <- mxAlgebra(expression=cbind(b0,b0), name='expMean')
#
# ADE model
# Create Matrices for Variance Components
#
SA      <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svVA, label=c("SA11","SA21","SA31","SA41","SA22","SA32","SA42","SA33","SA43","SA44"), 
        name="VA" )
SE      <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svVE, label=c("SE11","SE21","SE31","SE41","SE22","SE32","SE42","SE33","SE43","SE44"), 
         name="VE" )
# 
SD      <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svVD, label=c("SD11","SD21","SD31","SD41","SD22","SD32","SD42","SD33","SD43","SD44"), 
name="VD" )
#
RA <- mxAlgebra( expression = cov2cor(VA), name="RA")
RE <- mxAlgebra( expression = cov2cor(VE), name="RE")
RD <- mxAlgebra( expression = cov2cor(VD), name="RD")
#
covP      <- mxAlgebra( expression= VA+VD+VE, name="V" )
covMZ     <- mxAlgebra( expression= VA+VD, name="cMZ" )
covDZ     <- mxAlgebra( expression= 0.5%x%VA+.25%x%VD, 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" )
#ADE model
# Create Data Objects for Multiple Groups
dataMZ    <- mxData( observed=N4mz, type="raw" )
dataDZ    <- mxData( observed=N4dz, type="raw" )
# Create Expectation Objects for Multiple Groups
expMZ     <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars )
expDZ     <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars )
funML     <- mxFitFunctionML()
pars      <- list(B0_, SA, SD, SE, RA, RE, RD, covP )
modelMZ   <- mxModel( pars, ExpMean,covMZ, expCovMZ, ExpMean, dataMZ, expMZ, funML, name="MZ" )
modelDZ   <- mxModel( pars, ExpMean, covDZ, expCovDZ, ExpMean, dataDZ, expDZ, funML, name="DZ" )	
multi     <- mxFitFunctionMultigroup( c("MZ","DZ") )
#
modelADE4V  <- mxModel("ADE4V", pars, modelMZ, modelDZ, multi)
# ------------------------
# Run ADE Model
fitADE4V    <- mxTryHard( modelADE4V, 20)
sumADE4V    <- summary( fitADE4V )
# ----------------------------------------------------------- 
# end ADE Model version 1
# <-
# <- page 25/26
# ############# ----------------------------------------------
# -> page 26
# ->
mxCompare(fitSAT4B, fitADE4V)
# <-
# <- page 26 
# ############# ----------------------------------------------
# -> page 26 / q 3.3
modelAE4V <- omxSetParameters(fitADE4V,
label=c("SD11","SD21","SD31","SD41","SD22","SD32","SD42","SD33","SD43","SD44"),                                 
       free=FALSE, values=c(0))
fitAE4V    <- mxTryHard( modelAE4V, 20)
# <-  page 26
# ############# ----------------------------------------------
# -> page  26 / q 3.5
mxCompare(fitADE4V, fitAE4V) # -------------------------- version 2
# <- page 26
# ############# ----------------------------------------------
#
# -> page 27
# ->
rA=round(fitAE4V$RA$result,3)   # A correlations 
sA=diag(sqrt(round(fitAE4V$VA$values,3)))  # A standard deviations
SA=round(fitAE4V$VA$values,3) # A cov matrix
#
rE=round(fitAE4V$RE$result,3)  # E correlations
sE=diag(sqrt(round(fitAE4V$VE$values,3))) # E standard deviations
SE=round(fitAE4V$VE$values,3) # E cov matrix 
rA # correlation matrix A
sA
SA
rE
sE
SE
# <-
# <- page 27
# ->
SA_=diag(sA)%*%rA%*% diag(sA)  # A covariance matrix
SE_= diag(sE)%*%rE%*% diag(sE)     # E covariance matrix
#SA_
#SE_
SPh=SA_+SE_  # expected phenotypic covariance matrix
round(SA_/SPh,3)       # proportions
round(SE_/SPh,3)      # proportions 
# <-
# <- page 27
# ############# ----------------------------------------------
# -> page 28
LA=matrix(c(.6,.5,.7,.6),4,1)
TA=diag(c(.65,.75,.55,.65))
S_ = LA%*%t(LA) + TA
#
LA   # Factor loadings
TA    #  Residuals
S_      # expected covariance matrix
cov2cor(S_)  # correlation matrix
# <- page 28
### ---------------------------------------
# -> page 29
# ->
# AE IPM
nv=4
#
svb0=c(2.5,2.5,2.5,2.5)   
#
# Create regression model for expected Mean Matrices
B0_     <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=svb0, labels=c("b01","b02","b03","b04"), name="b0" )
#
ExpMean <- mxAlgebra(expression=cbind(b0,b0), name='expMean')
#AE model
# Create Matrices for Variance Components
LA      <- mxMatrix( type="Full", nrow=nv, ncol=1, free=TRUE, values=c(.5,.5,.5,.5), label=c("LA1","LA2","LA3","LA4"), name="LA" )
LD      <- mxMatrix( type="Full", nrow=nv, ncol=1, free=FALSE, values=0, label=c("LD1","LD2","LD3","LD4"), name="LD" )
LE      <- mxMatrix( type="Full", nrow=nv, ncol=1, free=TRUE, values=c(.7,.7,.7,.7), label=c("LE1","LE2","LE3","LE4"), name="LE"  )
#
TA      <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=c(.5,.5,.5,.5), label=c("tA1","tA2","tA3","tA4"), name="TA" )
TD      <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=FALSE, values=0, label=c("tD1","tD2","tD3","tD4"), name="TD" )
TE      <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=c(.7,.7,.7,.7), label=c("tE1","tE2","tE3","tE4"), name="TE" )
#
# Create Algebra for expected Variance/Covariance Matrices in MZ & DZ twins
#
VA      <- mxAlgebra( expression= LA%*%t(LA) + TA, name="VA" )
VD      <- mxAlgebra( expression= LD%*%t(LD) + TD, name="VD" )
VE      <- mxAlgebra( expression= LE%*%t(LE) + TE, name="VE" )
#
covP      <- mxAlgebra( expression= VA+VD+VE, name="V" )
covMZ     <- mxAlgebra( expression= VA+VD, name="cMZ" )
covDZ     <- mxAlgebra( expression= 0.5%x%VA+.25%x%VD, 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" )
#AE IPM model
dataMZ    <- mxData( observed=N4mz, type="raw" )
dataDZ    <- mxData( observed=N4dz, type="raw" )
#AE IPM model
# Create Expectation Objects for Multiple Groups
expMZ     <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars )
expDZ     <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars )
funML     <- mxFitFunctionML()
#AE IPM model
pars      <- list(B0_, LA, TA, LD, TD, LE, TE, VA, VD, VE, covP)
modelMZ   <- mxModel( pars, ExpMean,  covMZ, expCovMZ, ExpMean, dataMZ, expMZ, funML, name="MZ" )
modelDZ   <- mxModel( pars, ExpMean,   covDZ, expCovDZ, ExpMean, dataDZ, expDZ, funML, name="DZ" )	
multi     <- mxFitFunctionMultigroup( c("MZ","DZ") )
#AE IPM model
modelAE4IPM  <- mxModel('AE_IPM', pars, modelMZ, modelDZ, multi)
# ----------------------------------------------------------------------------------------------------------------------
# Run AE IPM Model
fitAE4IPM    <- mxTryHard( modelAE4IPM, 20)
sumAE4IPM    <- summary( fitAE4IPM )
# <-
# <- page 29
# ############# ----------------------------------------------
# -> page 30
# ->
mxCompare(fitAE4V,fitAE4IPM) 
# <-
# <- page 30
#
# -> page 30
estLA=round(fitAE4IPM$LA$values,3)  # A loadings 
estTA=round(fitAE4IPM$TA$values,3)  # A residual variances
estLE=round(fitAE4IPM$LE$values,3)  # E loadings
estTE=round(fitAE4IPM$TE$values,3)  # E residual variances
#
t(estLA)     # A loadings 
diag(estTA) # A residual variances
t(estLE)    # E loadings
diag(estTE)  # E residual variances
SA=estLA%*%t(estLA) + estTA  # A cov matrix
SE=estLE%*%t(estLE) + estTE  # E cov matrix
SPh=SA+SE     # expected phenotypic cov matrix
round(SA/SPh,3)   # proportions
round(SE/SPh,3)  # proportions
# <- page 30
# ############# ----------------------------------------------
# -> page 31
LA_res<-fitAE4IPM$LA$values # common A factor loadings
cvarA<-LA_res%*%t(LA_res) # common A variance matrix
TA_res<-diag(fitAE4IPM$TA$values) # residual A variance vector
totalA<-cvarA+diag(TA_res) # total A variance matrix
#
LE_res<-fitAE4IPM$LE$values # common E factor loadings
TE_res<-diag(fitAE4IPM$TE$values) # common E variance matrix
cvarE<-LE_res%*%t(LE_res) # residual E variance vector
totalE<-cvarE+diag(TE_res) # total E variance matrix
#
totalPh<-totalA +totalE # total phenotypic variance matrix
# <- page 31
# -> page 31
diag(totalA)/diag(totalPh) # narrow sense h2
diag(totalE)/diag(totalPh) # Env
# <- page 31
# -> page 31
commonA<-floor(diag(cvarA/totalPh)*100)
uniqueA<-floor(diag(TA_res/totalPh)*100)
commonE<-floor(diag(cvarE/totalPh)*100)
uniqueE<-floor(diag(TE_res/totalPh)*100)
# <- page 31
# ############# ----------------------------------------------
# -> page 32
Trait<-as.factor(rep(c(1,2,3,4),4))
Source<-c(rep("Common A",nv),rep("Common E",nv),
          rep("Unique A",nv),rep("Unique E",nv))
VC<-rep("Total variance",16)
Perc<-c(commonA,commonE,uniqueA,uniqueE)
data<-data.frame(Trait,Source,VC,Perc)
#
# put the components in the right order
data$Source <- factor(data$Source, levels = c("Unique E","Common E","Unique A", "Common A"))   
# <- page 32
# -> page 32
library(ggplot2)
p1 <- ggplot(data=data, aes(x=Trait, y=Perc, fill=Source )) + 
  geom_bar(stat="identity")  +   
  facet_wrap(~factor(VC), scales = "free_x")    +                  
  xlab(" ") + ylab("% of variance") +
  scale_fill_manual(name = " ",values = c("Purple","Violet","Red","tomato")) +
  ylim(0,100) +
  theme_classic(base_size = 17, base_family = "sans") +
  theme(legend.position = "bottom",
   axis.text.x = element_text(angle = 45, vjust = 1,
   size = 14, hjust = 1),legend.text=element_text(size=18))
p1
# <- page 32
# ############# ----------------------------------------------
# page 35
# ->
# -> page 37
# CPM 2
nv=4
svb0=c(2,2,2,2)   
#
# Create regression model for expected Mean Matrices
B0_     <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=svb0, labels=c("b01","b02","b03","b04"), name="b0" )
#
ExpMean <- mxAlgebra(expression=cbind(b0,b0), name='expMean')
# Create Matrices for Variance Components
LF      <- mxMatrix( type="Full", nrow=nv, ncol=1, free=TRUE, values=c(.5,.5,.5,.5), label=c("Lf1","Lf2","Lf3","Lf4"), name="LF" )
SA      <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=c(.5), label=c("sA"), name="SA" )
SD     <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, values=c(.0), label=c("sD"), name="SD" )
SE      <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=c(.5), label=c("sE"), name="SE" )
# scaling constraint
con_stand <- mxConstraint(expression= (sA+sD+sE) == 1, name='standLF')
# residual variances
TA      <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=c(.5,.5,.5,.5), label=c("tA1","tA2","tA3","tA4"), name="TA" )
TD      <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=FALSE, values=0, label=c("tD1","tD2","tD3","tD4"), name="TD" )
TE      <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=c(.7,.7,.7,.7), label=c("tE1","tE2","tE3","tE4"), name="TE" )
# Create Algebra for expected Variance/Covariance Matrices in MZ & DZ twins
VA      <- mxAlgebra( expression= LF%*%(SA)%*%t(LF)  + TA, name="VA" )
VD      <- mxAlgebra( expression= LF%*%(SD)%*%t(LF) + TD, name="VD" )
VE      <- mxAlgebra( expression= LF%*%(SE)%*%t(LF) + TE, name="VE" )
#
covP      <- mxAlgebra( expression= VA+VD+VE, name="V" )
covMZ     <- mxAlgebra( expression= VA+VD, name="cMZ" )
covDZ     <- mxAlgebra( expression= 0.5%x%VA+.025%x%VD, 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
#ACDE model
dataMZ    <- mxData( observed=N4mz, type="raw" )
dataDZ    <- mxData( observed=N4dz, type="raw" )
#ACDE model
# Create Expectation Objects for Multiple Groups
expMZ     <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars )
expDZ     <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars )
funML     <- mxFitFunctionML()
# 
pars      <- list(B0_, LF, TA, TD, TE, SA, SD, SE, VA, VD,VE, covP)
modelMZ   <- mxModel( pars, ExpMean,  covMZ, expCovMZ, ExpMean, dataMZ, expMZ, funML, name="MZ" )
modelDZ   <- mxModel( pars, ExpMean,   covDZ, expCovDZ, ExpMean, dataDZ, expDZ, funML, name="DZ" )	
multi     <- mxFitFunctionMultigroup( c("MZ","DZ") )
# model
modelAE4CPM2  <- mxModel("AE_CPM2", pars, modelMZ, modelDZ, multi,con_stand )
# RUN MODEL
fitAE4CPM2    <- mxTryHard( modelAE4CPM2, 20)
sumAE4CPM2    <- summary( fitAE4CPM2 )
# <- page 37
#
# -> page 37
mxCompare(fitAE4IPM,fitAE4CPM2)
# <- page 37
# ############# ----------------------------------------------
# -> page 38
FL=t(round(fitAE4CPM2$LF$values,3))
TA=diag(round(fitAE4CPM2$TA$values,3))
TE=diag(round(fitAE4CPM2$TE$values,3))
sA2=round(fitAE4CPM2$SA$values,3)
sE2=round(fitAE4CPM2$SE$values,3)
#
print(c(sA2, sE2))
print(FL)
print(TA)
print(TE)
# <- page 38
# ############# ----------------------------------------------
 # -> page 39
Nmz=dim(N4mz)[1]
Ndz=dim(N4dz)[1]
sumNmz=matrix(0,Nmz,2)
sumNmz[,1]=N4mz[,1]+N4mz[,2]+N4mz[,3]+N4mz[,4]
sumNmz[,2]=N4mz[,5]+N4mz[,6]+N4mz[,7]+N4mz[,8]
sumNdz=matrix(0,Ndz,2)
sumNdz[,1]=N4dz[,1]+N4dz[,2]+N4dz[,3]+N4dz[,4]
sumNdz[,2]=N4dz[,5]+N4dz[,6]+N4dz[,7]+N4dz[,8]
#
colnames(sumNdz)=c('sumN1','sumN2')
sumNdz=as.data.frame(sumNdz)
colnames(sumNmz)=c('sumN1','sumN2')
sumNmz=as.data.frame(sumNmz)
#
describe(sumNmz, range=F, skew=F)
describe(sumNdz, range=F, skew=F)
# <- page 39
# ############# ----------------------------------------------
# -> page 40
# univariate model AE model
# starting values
selVars=c('sumN1','sumN2')
nv=1  # number of phenotypes
ntv=2   # times number of twins
svVA=10*.4   #  A variance
svVE=10*.6   #  E variance
svVD=10*.0    # D variance = 0
svb0=10   # mean approx 
# Create regression model for expected Mean Matrices
ExpMean    <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svb0, labels=c("b0","b0"), name="expMean" )
# ADE model
# Create Matrices for Variance Components
covA      <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svVA, label="VA11", name="VA" ) 
covD      <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=FALSE, values=svVD, label="VD11", name="VD" )
covE      <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svVE, label="VE11", name="VE" )
# Create Algebra for expected Variance/Covariance Matrices in MZ & DZ twins
covP      <- mxAlgebra( expression= VA+VD+VE, name="V" )
covMZ     <- mxAlgebra( expression= VA+VD, name="cMZ" )
covDZ     <- mxAlgebra( expression= 0.5%x%VA+.25%x%VD, 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=sumNmz, type="raw" )
dataDZ    <- mxData( observed=sumNdz, type="raw" )
# Create Expectation Objects for Multiple Groups
expMZ     <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars )
expDZ     <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars )
funML     <- mxFitFunctionML()
pars      <- list(covA, covD, covE, covP )
modelMZ   <- mxModel( pars, ExpMean, covMZ, expCovMZ, ExpMean, dataMZ, expMZ, funML, name="MZ" )
modelDZ   <- mxModel( pars, ExpMean,  covDZ, expCovDZ, ExpMean, dataDZ, expDZ, funML, name="DZ" )	
multi     <- mxFitFunctionMultigroup( c("MZ","DZ") )
# Create Algebra for Variance Components
rowVC     <- rep('VarComp',nv)
colVC     <- rep(c('VA','VD','VE','SA','SD','SE'),each=nv)
estVC     <- mxAlgebra( expression=cbind(VA,VD,VE,VA/V,VD/V,VE/V), name="VarComp", dimnames=list(rowVC,colVC) )
ciAE     <- mxCI( "VarComp[1,1:8]" )
# Build Model with Confidence Intervals
modelAE  <- mxModel('univAE', pars, modelMZ, modelDZ, multi, estVC, ciAE )
# Run ADE Model
fitAE    <- mxRun( modelAE, intervals=F)
sumAE    <- summary( fitAE )
# <- page 40 /41
# ############# ----------------------------------------------
# -> page 41
fitAE$VarComp
# <- page 41

# ############# ----------------------------------------------
# ############# THE END
# ############# ----------------------------------------------