# ------------------------------------------------------------------------------
# Program: twinAceSibs.R  
# Author: Sarah Medland 
# Date: 13.04.2012 
#
# Univariate Twin Analysis model to estimate causes of variation
# Matrix style model - Raw data - Continuous data
# -------|---------|---------|---------|---------|---------|---------|---------|

# Load Library
require(OpenMx)
mxOption(NULL, "Default optimizer", "NPSOL") 
# ------------------------------------------------------------------------------
# PREPARE DATA

# Load Data

mzData <- read.csv("mzData.csv")
dzData <- read.csv("dzData.csv")



# Select Variables for Analysis
nv        <- 2       # number of variables
nind      <- 2
ntv       <- nv*nind    # number of total variables
selVars   <- c("Score_t1", "Disease_t1", "Score_t2", "Disease_t2")


defVars   <- c("ageTwins")

useVars   <- c(selVars,defVars)

# Set the Binary variable up for OpenMx
mzData$Disease_t1 <- mxFactor(mzData$Disease_t1, levels=c(0:1) )
mzData$Disease_t2 <- mxFactor(mzData$Disease_t2, levels=c(0:1) )

dzData$Disease_t1 <- mxFactor(dzData$Disease_t1, levels=c(0:1) )
dzData$Disease_t2 <- mxFactor(dzData$Disease_t2, levels=c(0:1) )

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

# Age correction
# Regression effects
Beta1        <- mxMatrix( type="Full", nrow=1, ncol=2, free=TRUE, 
                          values=.1, labels=c('bAgeCon', 'bAgeBin'), name="beta1" )
# Independent variable
obsAge       <- mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, 
                          labels=c("data.ageTwins", "data.ageTwins"), name="Age")

Beta2        <- mxMatrix( type="Full", nrow=1, ncol=2, free=TRUE, 
                          values=.1, labels=c('bSexCon', 'bSexBin'),
                          name="beta2" )
# Independent variable
obsSex       <- mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, 
                          labels=c("data.sex_t1", "data.sex_t2"),
                          name="Sex")

#setting up the regression
intercept      	<- mxMatrix( type="Full", nrow=1, ncol=ntv, free=c(T,F),
                             values=c(3,0), 
                           labels=c("meanP","binary"), name="intercept" )

threshold       <-mxMatrix( type="Full", nrow=1, ncol=nind, free=T, 
                            values=1, 
                            labels="thresh", name="Threshold" )

expMean  	<- mxAlgebra( intercept + (Age %x% beta1) + (Sex %x% beta2) , name="expMean")


inclusions   <- list (Beta1, Beta2,  obsAge, obsSex, intercept, expMean, threshold)

# ------------------------------------------------------------------------------

# ACE Model
# Matrices declared to store A, C, and E Variances
A     <- mxMatrix( type="Symm", nrow=nv, ncol=nv, 
                       free=TRUE, values=.6, label=c("a11","a21","a22"), 
                       name="A" ) 
C     <- mxMatrix( type="Symm", nrow=nv, ncol=nv, 
                       free=TRUE, values=.2, label=c("c11","c21","c22"), 
                       name="C" )
E     <- mxMatrix( type="Symm", nrow=nv, ncol=nv, 
                       free=c(T,T,F), values=1, label=c("e11","e21","e22"),
                       name="E" )
	

# Algebra for expected Variance/Covariance Matrices in MZ & DZ twins
covMZ     <- mxAlgebra( expression= rbind( 
                              	    cbind(A+C+E , 	     A+C),
                                        cbind(A+C,	       A+C+E)), name="expCovMZ" )
covDZ     <- mxAlgebra( expression= rbind( 
					          cbind(A+C+E , 	.5%x%A+C),
                                        cbind(.5%x%A+C,	   A+C+E)), name="expCovDZ" )

Imat        <- mxMatrix(name= "I", type="Iden", nrow = nv, ncol = nv)
# Data objects for Multiple Groups
dataMZ    <- mxData( observed=mzData, type="raw" )
dataDZ    <- mxData( observed=dzData, type="raw" )

# Objective objects for Multiple Groups
expMZ     <- mxExpectationNormal( covariance="expCovMZ", means="expMean",
                              dimnames=selVars,
                              thresholds="Threshold",
                              threshnames=c('Disease_t1','Disease_t2')  )


expDZ     <- mxExpectationNormal( covariance="expCovDZ", means="expMean", 
                              dimnames=selVars ,
                              thresholds="Threshold", 
                              threshnames=c('Disease_t1','Disease_t2')  )

# Combine Groups
pars      <- list( A,C,E, Imat)
funML     	<- 	mxFitFunctionML()

modelMZ   <- mxModel( pars, inclusions, covMZ, dataMZ, expMZ, funML, name="MZ" )
modelDZ   <- mxModel( pars, inclusions, covDZ, dataDZ, expDZ, funML, name="DZ" )
multi     	<- mxFitFunctionMultigroup( c("MZ","DZ") )
AceModel  <- mxModel( "ACE", pars, modelMZ, modelDZ, funML, multi )
AceModel  <-  mxTryHard(AceModel)
# ------------------------------------------------------------------------------
# RUN MODEL


# Run ACE model
AceFit    <- mxRun(AceModel)
AceSumm   <- summary(AceFit)
AceSumm
round(AceFit@output$estimate,4)


# ------------------------------------------------------------------------------
# FIT SUBMODELS

# Run AE model
AeModel   <- mxModel( AceFit, name="AE" )
AeModel   <- omxSetParameters( AeModel, labels="c11", free=FALSE, values=0 )
AeModel   <- omxSetParameters( AeModel, labels="c21", free=FALSE, values=0 )
AeModel   <- omxSetParameters( AeModel, labels="c22", free=FALSE, values=0 )
AeFit     <- mxRun(AeModel)
round(AeFit@output$estimate,4)