#install.packages("psych")
setwd("\\\\ATKK/home/s/silvento/data/opetus/OpenMx/omat esitykset")
require(OpenMx)
require(psych)

#Prepare data
BMIdata <- read.table ("bmidata.dat", header=T, na=-99)
head(BMIdata, n=10)
#dim(BMIdata)
#describe(BMIdata)
vars <- 'bmi16'
nv <- 1
ntv <- nv*2
selVars <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="")
mzmData <- subset(BMIdata, zyg==2 & sex1==0, selVars)
dzmData <- subset(BMIdata, zyg==4 & sex1==0, selVars)
#str(dzmfData)


# ------------------------------------------------------------------------------
# PREPARE SATURATED MODEL

### 1. MEANS ###
# Algebra for expected Mean Matrices in the five zygosity groups
meanMZM    <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=round(colMeans(mzmData,na.rm=TRUE),digits=0),# [mean twin1, mean twin2]
                        labels=c("meanT1mzm","meanT2mzm"), name="expMeanMZM" )                                        # use labels to set parameters to be equal (same name=same estimate)
meanDZM    <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=round(colMeans(dzmData,na.rm=TRUE),digits=0),# but: starting values must be the same too!
                        labels=c("meanT1dzm","meanT2dzm"), name="expMeanDZM" )



### 2. (CO-)VARIANCES ###
# Determine starting values, (co-)variances
startCovMZM     <- round(vech(cov(mzmData,use="pairwise.complete.obs")),2) #extracts vartwin1-covtwin12-vartwin2
startCovDZM     <- round(vech(cov(dzmData,use="pairwise.complete.obs")),2)

print(startCovMZM)

# Algebra for expected Variance/Covariance Matrices in the five zygosity groups
covMZM  <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, label=c("varMZT1","covMZ","varMZT2"), values=c(4,3,4), name="expCovMZM" ) 
covDZM  <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, label=c("varDZT1","covDZ","varDZT2"), values=c(4,3,4), name="expCovDZM" ) 

print(covMZM)

### 3. CORRELATIONS ### (= standardised results)
Imat <- mxMatrix( type="Iden", nrow=ntv, ncol=ntv, name="I")

isdMZM <- mxAlgebra( expression = solve(sqrt(I*expCovMZM)), name="iSDmzm") # "solve" inverts the matrix [so it becomes e.g. 1/sqrt(var1)]
isdDZM <- mxAlgebra( expression = solve(sqrt(I*expCovDZM)), name="iSDdzm")

corMZM <- mxAlgebra( expression = iSDmzm%*%expCovMZM%*%iSDmzm, name="expCorMZM")
corDZM <- mxAlgebra( expression = iSDdzm%*%expCovDZM%*%iSDdzm, name="expCorDZM")

# Data objects for the five zygosity groups
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" )

# Objective functions for the five zygosity groups
objMZM     <- mxExpectationNormal( covariance="expCovMZM", means="expMeanMZM", dimnames=selVars ) # uses full-information maximum likelihood
objDZM     <- mxExpectationNormal( covariance="expCovDZM", means="expMeanDZM", dimnames=selVars )

fitFunction <- mxFitFunctionML()

# Specify model
modelMZM   <- mxModel( "MZM", meanMZM, covMZM, Imat, isdMZM, corMZM, dataMZM, objMZM, fitFunction ) # put everything into the model
modelDZM   <- mxModel( "DZM", meanDZM, covDZM, Imat, isdDZM, corDZM, dataDZM, objDZM, fitFunction )
 
minus2ll       <- mxAlgebra(expression= MZM.objective + DZM.objective, name="minus2sumloglikelihood" ) # combine the five groups
obj            <- mxFitFunctionAlgebra( "minus2sumloglikelihood" )
ciCor          <- mxCI( c('MZM.expCorMZM','DZM.expCorDZM' ))

twinSatModel   <- mxModel( "twinSat", modelMZM, modelDZM, minus2ll, obj, ciCor)

# ------------------------------------------------------------------------------
# RUN SATURATED MODELS

twinSatFit     <- mxRun( twinSatModel, intervals=F ) # it usually takes long to calculate CIs -> first do "intervals=F" to see whether the script is working at all
summary( twinSatFit )


# Fix variances to be same for 1. and 2. twins  

eqSDTwinModel    <- mxModel(twinSatFit, name="eqSDTwin" )
eqSDTwinModel     <- omxSetParameters( eqSDTwinModel, label=c("varMZT1","varMZT2"), free=TRUE, # label=old labels -> newlabels
                                         values=c(5,5), newlabels=c("varMZ","varMZ") )                # starting values have to be identical too, choose based on estimates!

eqSDTwinModel     <- omxSetParameters( eqSDTwinModel, label=c("varDZT1","varDZT2"), free=TRUE, # label=old labels -> newlabels
                                         values=c(5,5), newlabels=c("varDZ","varDZ") )                # starting values have to be identical too, choose based on estimates!               

eqSDTwinFit      <- mxRun( eqSDTwinModel, intervals=F)
summary( eqSDTwinFit)

mxCompare(twinSatFit,eqSDTwinFit)


# Fix additionally variances to be same for MZ and DZ twins 

eqSDTwin2Model    <- mxModel(eqSDTwinFit, name="eqSD2Twin" )
eqSDTwin2Model     <- omxSetParameters( eqSDTwin2Model, label=c("varMZ","varDZ"), free=TRUE, # label=old labels -> newlabels
                                         values=c(5,5), newlabels=c("var","var") )                # starting values have to be identical too, choose based on estimates!
   

eqSDTwin2Fit      <- mxRun(eqSDTwin2Model, intervals=F)
summary( eqSDTwin2Fit)

mxCompare(twinSatFit,eqSDTwin2Fit)



# Fix means to be same for 1. and 2. twins 
 
eqMEANTwinModel    <- mxModel(twinSatFit, name="eqMEANTwin" )
eqMEANTwinModel     <- omxSetParameters( eqMEANTwinModel, label=c("meanT1mzm","meanT2mzm"), free=TRUE, # label=old labels -> newlabels
                                         values=c(20,20), newlabels=c("meanmzm","meanmzm") )                # starting values have to be identical too, choose based on estimates!

eqMEANTwinModel      <- omxSetParameters( eqMEANTwinModel, label=c("meanT1dzm","meanT2dzm"), free=TRUE, 
                                         values=c(20,20), newlabels=c("meandzm","meandzm") )                

eqMEANTwinFit      <- mxRun( eqMEANTwinModel, intervals=F)
summary(eqMEANTwinFit)

mxCompare(twinSatFit,eqMEANTwinFit)


# Fix additionally means to be same for MZ and DZ twins 

eqMEANTwin2Model    <- mxModel(eqMEANTwinFit, name="eqMEAN2Twin" )
eqMEANTwin2Model     <- omxSetParameters( eqMEANTwinModel, label=c("meanmzm","meandzm"), free=TRUE, # label=old labels -> newlabels
                                         values=c(20,20), newlabels=c("mean","mean") )                # starting values have to be identical too, choose based on estimates!

eqMEANTwin2Fit      <- mxRun( eqMEANTwin2Model, intervals=F)
summary( eqMEANTwin2Fit)

mxCompare(twinSatFit,eqMEANTwin2Fit)