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


#Data objects

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

#ACE model

pathAm <- mxMatrix(type="Full", nrow=nv, ncol=nv, free=TRUE, values=7, 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=7, 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=7, label="em11", name="em")  #creates Mx matrix for specific environmental factors for males (1*1 matrix)

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

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

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

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

# Objective objects for Multiple Groups

objMZM  <- mxExpectationNormal( covariance="expCovMZM", means="expMeanMZM", dimnames=selVars )
objDZM  <- mxExpectationNormal( covariance="expCovDZM", means="expMeanDZM", dimnames=selVars )

# Choose a fit function

fitFunction <- mxFitFunctionML()

pars      <- list( pathAm, pathCm, pathEm, covAm, covCm, covEm, VarM) #creates a vector including these elements

ACEestimates   <- mxAlgebra( expression=cbind(Am/Vm,Cm/Vm,Em/Vm), name="VarsACE",) #divide genetic and environmental variances by total variance 
                     
ciVars    <- mxCI("VarsACE") 
 
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" )



minus2ll  <- mxAlgebra( expression=MZM.objective + DZM.objective, name="m2LL" )

  
  obj       <- mxFitFunctionAlgebra( "m2LL" )  
  
AceModel  <- mxModel( "ACE", pars, modelMZM, modelDZM, minus2ll, obj, ACEestimates, ciVars )



AceFit    <- mxRun(AceModel, intervals=F) #run the model
summary(AceFit) #presents the model fit results

#str(AceFit)

round(AceFit@output$estimate,4) # path estimates and means
pathEstimatesACE <- round(mxEval(cbind(Am/Vm,Cm/Vm,Em/Vm), AceFit),4) # ACE estimates
rownames(pathEstimatesACE) <- 'pathEstimates'                                           # labeling, for clarity
colnames(pathEstimatesACE) <- c('Am/Vm','Cm/Vm','Em/Vm')
pathEstimatesACE
variances <- mxEval(cbind(Am,Cm,Em), AceFit)
variances


#FIT AE SuBMODEL 
#AEModel   <- mxModel( AceFit, name="AE" )
#AEModel   <- omxSetParameters( AEModel, labels="cm11", free=FALSE, values=0 )
#AEFit     <- mxRun(AEModel, intervals=F )
#summary(AEFit)
#mxCompare(AceFit, AEFit)