#install.packages("psych") 
rm(list=ls(all=TRUE))

setwd("\\\\ATKK/home/s/silvento/data/opetus/OpenMx/omat esitykset")
require(OpenMx)
require(psych)#download OpenMx and psych packages is already installed


#Prepare data, data needs to be in wide format (one twin pair per line)

BMIdata <- read.table ("bmidata.dat", header=T, na=-99)
#head(BMIdata, n=10)
#dim(BMIdata)
#describe(BMIdata)

vars <- 'bmi16' #identify a variable used in the analyses (BMI18z)
nv <- 1 #number of variables used in the analyses (BMI18z)
ntv <- nv*2 #variables multiplied by 2 because one variable for each co-twin

selVars <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="") #creates a vector to select variables for both co-twins (BMI18z1 and BMI18z2)


mzmData <- subset(BMIdata, zyg==2, selVars) #creates a subset of data including MZM twins and variables BMI18z1 and BMI18z2 
dzmData <- subset(BMIdata, zyg==4, selVars)
mzfData <- subset(BMIdata, zyg==1, selVars)
dzfData <- subset(BMIdata, zyg==3, selVars)
dosmfData <- subset(BMIdata, zyg==5 & sex1==0 & sex2==1 , selVars)
dosfmData <- subset(BMIdata, zyg==5 & sex1==1 & sex2==0 , selVars)


dataMZM    <- mxData( observed=mzmData, type="raw" ) # creates MxData objects, will be used as argument in MxModel ("raw data in OpenMx format")
dataDZM    <- mxData( observed=dzmData, type="raw" )
dataMZF    <- mxData( observed=mzfData, type="raw" )
dataDZF    <- mxData( observed=dzfData, type="raw" )
dataDOSMF    <- mxData( observed=dosmfData, type="raw" )
dataDOSFM    <- mxData( observed=dosfmData, type="raw" )


#ACE model

pathAm <- mxMatrix(type="Full", nrow=nv, ncol=nv, free=TRUE, values=1, label="am11", name="am") #creates Mx matrix for additive genetic factors for males (1*1 matrix)
pathCm <- mxMatrix(type="Full", nrow=nv, ncol=nv, free=TRUE, values=1, label="cm11", name="cm")  #creates Mx matrix for common environmental factors for males (1*1 matrix)
pathEm <- mxMatrix(type="Full", nrow=nv, ncol=nv, free=TRUE, values=1, label="em11", name="em")  #creates Mx matrix for specific environmental factors for males (1*1 matrix)

pathAf <- mxMatrix(type="Full", nrow=nv, ncol=nv, free=TRUE, values=1, label="af11", name="af")
pathCf <- mxMatrix(type="Full", nrow=nv, ncol=nv, free=TRUE, values=1, label="cf11", name="cf")
pathEf <- mxMatrix(type="Full", nrow=nv, ncol=nv, free=TRUE, values=1, label="ef11", name="ef")

rados <- mxMatrix(type="Full", nrow=1, ncol=1, free=F, values=0.5, label="rados", lbound=-1, ubound=1, name="ra") #creates a freely estimated correlation parameter for OSDZ twins

covAm <- mxAlgebra (expression=am %*% t(am), name="Am") #calculate a new matrix: am times transpose(am); in principle takes the square of matrix
covCm <- mxAlgebra (expression=cm %*% t(cm), name="Cm")
covEm <- mxAlgebra (expression=em %*% t(em), name="Em")

covAf <- mxAlgebra (expression=af %*% t(af), name="Af")
covCf <- mxAlgebra (expression=cf %*% t(cf), name="Cf")
covEf <- mxAlgebra (expression=ef %*% t(ef), name="Ef")

VarM <- mxAlgebra (expression= Am+Cm+Em, name="Vm") #calculate total variance for males (additive genetic+common environmental + specific environmental variance)
VarF <- mxAlgebra (expression= Af+Cf+Ef, name="Vf")

CovMZM   <- mxAlgebra( expression= rbind( cbind(Vm, Am+Cm),                             
                                             cbind(Am+Cm, Vm)),       name="expCovMZM" ) #calculate expected variance-co-variance matrix for male MZ twins
CovDZM   <- mxAlgebra( expression= rbind( cbind(Vm, 0.5%x%Am+Cm),
                                             cbind(0.5%x%Am+Cm, Vm)), name="expCovDZM" )
CovMZF   <- mxAlgebra( expression= rbind( cbind(Vf, Af+Cf),       
                                             cbind(Af+Cf, Vf)),       name="expCovMZF" )
CovDZF   <- mxAlgebra( expression= rbind( cbind(Vf, 0.5%x%Af+Cf),
                                             cbind(0.5%x%Af+Cf, Vf)), name="expCovDZF" )
CovDOSMF   <- mxAlgebra( expression= rbind( cbind(Vm, ((ra%x%(am%*%t(af)))+(cm%*%t(cf)))),
                                             cbind((ra%x%(af%*%t(am))+(cf%*%t(cm))), Vf)), name="expCovDOSMF" )
CovDOSFM   <- mxAlgebra( expression= rbind( cbind(Vf, ((ra%x%(af%*%t(am)))+(cf%*%t(cm)))),
                                             cbind((ra%x%(am%*%t(af))+(cm%*%t(cf))), Vm)), name="expCovDOSFM" )


meanMZM      <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=c(0,0), labels=c("meanM","meanM"), name="expMeanMZM" ) # creates mean model for male MZ twins
meanDZM      <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=c(0,0), labels=c("meanM","meanM"), name="expMeanDZM" )
meanMZF      <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=c(0,0), labels=c("meanF","meanF"), name="expMeanMZF" )
meanDZF      <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=c(0,0), labels=c("meanF","meanF"), name="expMeanDZF" )
meanOSDZMF      <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=c(0,0), labels=c("meanM","meanF"), name="expMeanOSDZMF" )
meanOSDZFM      <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=c(0,0), labels=c("meanF","meanM"), name="expMeanOSDZFM" )



relativeAm <- mxAlgebra(expression=Am/Vm,name="relAm")
relativeCm <- mxAlgebra(expression=Cm/Vm,name="relCm")
relativeEm <- mxAlgebra(expression=Em/Vm,name="relEm")
relativeAf <- mxAlgebra(expression=Af/Vf,name="relAf")
relativeCf <- mxAlgebra(expression=Cf/Vf,name="relCf")
relativeEf <- mxAlgebra(expression=Ef/Vf,name="relEf")



# Objective objects for Multiple Groups
objMZM     <- mxExpectationNormal( covariance="expCovMZM", means="expMeanMZM", dimnames=selVars ) #estimate free parameter values using maximum likelihood estimator
objDZM     <- mxExpectationNormal( covariance="expCovDZM", means="expMeanDZM", dimnames=selVars )
objMZF     <- mxExpectationNormal( covariance="expCovMZF", means="expMeanMZF", dimnames=selVars )
objDZF     <- mxExpectationNormal( covariance="expCovDZF", means="expMeanDZF", dimnames=selVars )
objDOSMF     <- mxExpectationNormal( covariance="expCovDOSMF", means="expMeanOSDZMF", dimnames=selVars )
objDOSFM     <- mxExpectationNormal( covariance="expCovDOSFM", means="expMeanOSDZFM", dimnames=selVars )

fitFunction <- mxFitFunctionML()

pars      <- list( pathAm, pathCm, pathEm, pathAf, pathCf, pathEf, rados, covAm, covCm, covEm, covAf, covCf, covEf, VarM, VarF,
                   relativeAm, relativeCm, relativeEm,relativeAf, relativeCf, relativeEf ) #creates a vector including these elements


ACEestimates   <- mxAlgebra( expression=cbind(Am/Vm,Cm/Vm,Em/Vm,Af/Vf,Cf/Vf,Ef/Vf), name="VarsACE",) #divide genetic and environmental variances by total variance
                     

 ciVars    <- mxCI("VarsACE") #calculate 95% confidence intervals for standardized estimates

modelMZM   <- mxModel( pars, meanMZM, CovMZM, dataMZM, objMZM, fitFunction, name="MZM" )  #creates a model for male MZ twins
modelDZM   <- mxModel( pars, meanDZM, CovDZM, dataDZM, objDZM, fitFunction, name="DZM" )
modelMZF   <- mxModel( pars, meanMZF, CovMZF, dataMZF, objMZF, fitFunction, name="MZF" )
modelDZF   <- mxModel( pars, meanDZF, CovDZF, dataDZF, objDZF, fitFunction, name="DZF" )
modelDOSMF   <- mxModel( pars, meanOSDZMF, CovDOSMF, dataDOSMF, objDOSMF, fitFunction, name="DOSMF" )
modelDOSFM   <- mxModel( pars, meanOSDZFM, CovDOSFM, dataDOSFM, objDOSFM, fitFunction, name="DOSFM" )

minus2ll  <- mxAlgebra( expression=MZM.objective + DZM.objective + MZF.objective + DZF.objective + DOSMF.objective + DOSFM.objective , name="m2LL" )
obj       <- mxFitFunctionAlgebra( "m2LL" )
AceModel  <- mxModel( "ACE", pars, modelMZM, modelDZM, modelMZF, modelDZF, modelDOSMF, modelDOSFM, minus2ll, obj,ACEestimates,ciVars )


#str(modelMZM)
AceFit    <- mxRun(AceModel, intervals=F) #run the model
summary(AceFit) #presents the model fit results
AceFit$algebras$VarsACE

# FIT a SUBMODEL WIT NO SEX  DIFFERENCES IN STANDARDISED ESTIMATES

eqSexStanAceModel <- mxModel(AceModel, name="eqSexStanAce",
                        mxConstraint(expression=relAm==relAf, name="relAsame"))

eqSexStanAceModel <- mxModel(eqSexStanAceModel, name="eqSexStanAce",
                        mxConstraint(expression=relCm==relCf, name="relCsame"))


eqSexStanAceFit   <- mxRun(eqSexStanAceModel, intervals=F)

summary(eqSexStanAceFit)

eqSexStanAceFit$algebras$VarsACE