#-------------------------------------------------------------------------------------
# smoking initiation and adiposity(wc defined)
# Bivariate Twin model for ordinal data
# Phenotypes
# Variable 1: 	whtr_g categorised in 2 classes: 0,1
# Variable 2: 	smk_cur categorised in 2 classes:0,1
# Type of data:	Raw ordinal data
# zygosity: 		0=MZ 1=DZ
# -------------------------------------------------------------------------------------
# Note: in reality you would rather analyse the continuous distributions
 
rm(list=ls())
#ls()

require(OpenMx)
require(psych)
library(foreign)
require(readstata13)

# change optimizer ########
#mxOption(NULL, "Default optimizer", "NPSOL")
#mxOption( NULL, 'Default optimizer' ,'SLSQP' )
mxOption(NULL, "Default optimizer", "CSOLNP")

#Bivdata <- read.table ('ADHDdiaage.csv', header=T, sep=",", na.strings="NA")
getwd("C:/Users/GGWSXY/Desktop/lyn_data")
#Bivdata <- as.data.frame(read.spss("shortsleep176pair.sav")) ##176¶Ô
Bivdata <-  read.dta13("input/merge/total_pair_matched5745.dta")

#Ãô¸ÐÐÔ·ÖÎöÅųýѪÌÇ¿ØÖƽÏÀíÏëºÍ½ÏºÃ
#Bivdata <- Bivdata[Bivdata$diacon1!=1,]
#Bivdata <- Bivdata[Bivdata$diacon1!=2,]
#Bivdata <- Bivdata[Bivdata$diacon2!=1,]
#Bivdata <- Bivdata[Bivdata$diacon2!=2,]


names(Bivdata)
summary(Bivdata)
str(Bivdata)
describe(Bivdata)
table(Bivdata$smk_cur1)
table(Bivdata$smk_cur2)
#hist(Bivdata$dia1); hist(Bivdata$dia2)

#Make the integer variables ordered factors (Note:because of the 'cuting' above, this is not necessary here)
Bivdata$smk_cur1 	<-mxFactor(Bivdata$smk_cur1, levels=c(0:1) )
Bivdata$smk_cur2 	<-mxFactor(Bivdata$smk_cur2, levels=c(0:1) )
Bivdata$whtr_g1 <- mxFactor(Bivdata$whtr_g1, levels= c(0:1))
Bivdata$whtr_g2 <- mxFactor(Bivdata$whtr_g2, levels= c(0:1))

mzData <- subset(Bivdata, zygosity2==0, c(whtr_g1,whtr_g2,smk_cur1,smk_cur2))
dzData <- subset(Bivdata, zygosity2==1, c(whtr_g1,whtr_g2,smk_cur1,smk_cur2))

vars		<-c('whtr_g','smk_cur')
selVars	<-c('whtr_g1', 'smk_cur1', 'whtr_g2', 'smk_cur2' )
useVars	<-c('whtr_g1', 'smk_cur1', 'whtr_g2', 'smk_cur2' )

#Make the integer variables ordered factors (Note:because of the 'cuting' above, this is not necessary here)

describe(mzData)
describe(dzData)

sink('result/03_SEM_results/s3b_biv_whtr_smkcur_cat22_byage.doc')

table(mzData$smk_cur1, mzData$smk_cur2)
table(dzData$smk_cur1, dzData$smk_cur2)
# To get the CTs for ADHD
table(mzData$smk_cur1, mzData$smk_cur2)
table(dzData$smk_cur1, dzData$smk_cur2)

sum(is.na(mzData))

sum(is.na(dzData))

mzData <- na.omit(mzData)

nv 	<- 2		# number of variables per twin
ntv <- nv*2	# number of variables per pair
nth	<- 1		# number of max thresholds

# 1) Fits a constrained Polychoric correlation model 
# TH same across twins but different across zyg groups
# Age effect is different acros variables, but same across thresholds within variables (if c>2)
# There is one overall rPH between var1-2 and the x-trait x-twin correlations are symmetric
# ------------------------------------------------------------------------------------------------------------------------------
library( polycor )
polychor( table(Bivdata$smk_cur1, Bivdata$whtr_g1) , std.err = T , ML=T )
polychor( table(mzData$whtr_g1, mzData$whtr_g2) , std.err = T , ML=T )
polychor( table(dzData$whtr_g1, dzData$whtr_g2) , std.err = T , ML=T )
polychor( table(mzData$smk_cur1, mzData$smk_cur2) , std.err = T , ML=T )
polychor( table(dzData$smk_cur1, dzData$smk_cur2) , std.err = T , ML=T )
polychor( table(mzData$whtr_g1, mzData$smk_cur1) , std.err = T , ML=T )
polychor( table(dzData$whtr_g1, dzData$smk_cur1) , std.err = T , ML=T )


# CREATE LABELS & START VALUES as objects(to ease specification)
LabThMZ	<-c('Tmz_11','Tmz_21')			# THs for var 1 and 2 for a twin individual (mz)
LabCorMZ	<-c('r21','rMZ1','MZxtxt','MZxtxt','rMZ2','r21')
LabThDZ	<-c('Tdz_11','Tdz_21')			# THs for var 1 and 2 for a twin individual (dz)
LabCorDZ	<-c('r21','rDZ1','DZxtxt','DZxtxt','rDZ2','r21')
LabCov	<-c('BageThdia', 'BageThdia', 'BageThadhd', 'BageThadhd') 
ThPat		<-c(T,F,T,F)
StCorMZ	<-c(0.04, .9, 0.04,0.04, 0.87, 0.04)
StCorDZ	<-c(0.04, .7, 0.03,0.03, .73, 0.04)
StTH	<-c(0.2, 0, 0.07, 0 )

# Define definition variables to hold the Covariates

# Matrix & Algebra for expected means (SND), Thresholds, effect of Age on Th and correlations
Mean		<-mxMatrix( type="Zero", nrow=1, ncol=ntv, name="M" )

Tmz		<-mxMatrix( type="Full", nrow=nth, ncol=nv, free=ThPat, values=StTH, lbound=c(-3,.001), ubound=3, labels=LabThMZ, name="ThMZ")
ThresMZ	<-mxAlgebra( expression= cbind(ThMZ, ThMZ), name="expThresMZ")
CorMZ		<-mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=StCorMZ, labels=LabCorMZ, lbound=-.99, ubound=.99, name="expCorMZ") 

Tdz		<-mxMatrix( type="Full", nrow=nth, ncol=nv, free=ThPat, values=StTH, lbound=c(-3,.001), ubound=3, labels=LabThDZ, name="ThDZ")
ThresDZ	<-mxAlgebra( expression= cbind(ThDZ , ThDZ), name="expThresDZ")
CorDZ		<-mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=StCorDZ, labels=LabCorDZ, lbound=-.99, ubound=.99, name="expCorDZ")

# Data objects for Multiple Groups
dataMZ   	<- mxData( observed=mzData, type="raw" )
dataDZ   	<- mxData( observed=dzData, type="raw" )

# Objective objects for Multiple Groups
objMZ    	<- mxExpectationNormal( covariance="expCorMZ", means="M", dimnames=selVars, thresholds="expThresMZ" )
objDZ    	<- mxExpectationNormal( covariance="expCorDZ", means="M", dimnames=selVars, thresholds="expThresDZ" )
fitFunction<-mxFitFunctionML()
# Combine Groups
modelMZ	<- mxModel( Mean, Tmz,  ThresMZ, CorMZ, dataMZ, objMZ, name="MZ",fitFunction )
modelDZ	<- mxModel( Mean, Tdz,  ThresDZ, CorDZ, dataDZ, objDZ, name="DZ",fitFunction )
minus2ll	<- mxAlgebra( expression=MZ.objective + DZ.objective, name="m2LL" )
obj		<- mxFitFunctionAlgebra(algebra = "m2LL")
Conf		<- mxCI (c ('MZ.expCorMZ[2,1]') )
SatModel 	<- mxModel( "Sat", modelMZ, modelDZ, minus2ll, obj, Conf )

# -------------------------------------------------------------------------------------------------------------------------------
# 1) RUN Saturated Model

SatFit	<- mxRun(SatModel, intervals=TRUE)
(SatSumm	<- summary(SatFit))

# Generate SatModel Output

round(SatFit@output$estimate,4)

SatFit$MZ.expCorMZ
SatFit$DZ.expCorDZ

(ExpTetraCorMZ	<-mxEval(MZ.expCorMZ, SatFit))
(ExpTetraCorDZ	<-mxEval(DZ.expCorDZ, SatFit))

# ---------------------------------------------------------------------------------------------------
# 2) Specify and Run Sub1Model:  equating Thresholds across MZ and DZ group

Sub1Model	<- mxModel(SatModel, name="sub1")
Sub1Model	<- omxSetParameters(Sub1Model, labels=LabThDZ, newlabels=LabThMZ, free=ThPat, values=StTH, lbound=c(-3,.001), ubound=3 )
Sub1Fit	<- mxRun(Sub1Model, intervals=TRUE)
(Sub1Summ	<- summary(Sub1Fit))
round(Sub1Fit@output$estimate,4)

# Generate Sub1Model Output
(ExpTetraCorMZ	<-mxEval(MZ.expCorMZ, Sub1Fit))
(ExpTetraCorDZ	<-mxEval(DZ.expCorDZ, Sub1Fit))
mxCompare(SatFit,Sub1Fit)

# change optimizer ########
#mxOption(NULL, "Default optimizer", "NPSOL")
mxOption( NULL, 'Default optimizer' ,'SLSQP' )
#mxOption(NULL, "Default optimizer", "CSOLNP")

# CREATE LABELS for Cholesky decomposition with a function
laLower	<- function(la,nv) { paste(la,rev(nv+1-sequence(1:nv)),rep(1:nv,nv:1),sep="_") }
LabTh		<-c('T_11','T_12')

# 3) Specify ACE Model, with ONE overall set of Thresholds with Age effects on thresholds
# ------------------------------------------------------------------------------------------

# CREATE LABELS for Cholesky decomposition with a function
laLower	<- function(la,nv) { paste(la,rev(nv+1-sequence(1:nv)),rep(1:nv,nv:1),sep="_") }
LabTh		<-c('T_11','T_12')


# Matrix & Algebra for expected means (SND), Thresholds, effect of Age on Th and correlations
Mean		<-mxMatrix( type="Zero", nrow=1, ncol=ntv, name="M" )

Tr		<-mxMatrix( type="Full", nrow=nth, ncol=nv, free=ThPat, values=StTH, labels=LabTh, lbound=c(-3,.001), ubound=3, name="Th")
 
Thres		<-mxAlgebra( expression= cbind(Th , Th ), name="expThres")


# Matrices declared to store a, c, and e Path Coefficients
pathA	<- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=c(sqrt(0.08),-.3,sqrt(0.08)), label=laLower("a",nv), name="a" ) 
pathC	<- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=c(sqrt(0.08),-.1,sqrt(0.08)), label=laLower("c",nv), name="c" )
pathE	<- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=c(sqrt(0.08),-.3,sqrt(0.08)), label=laLower("e",nv), name="e" )
	
# Matrices generated to hold A, C, and E components and total Variance
covA	<- mxAlgebra( expression=a %*% t(a), name="A" )
covC	<- mxAlgebra( expression=c %*% t(c), name="C" ) 
covE	<- mxAlgebra( expression=e %*% t(e), name="E" )
covP	<- mxAlgebra( expression=A+C+E, name="V" )


# 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     , 0.5%x%A+C),
                                       cbind(0.5%x%A+C , A+C+E)), name="expCovDZ" )

# Algebra to compute total variances and standard deviations (whtr_ggonal only)

matI      <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I")
invSD     <- mxAlgebra( expression=solve(sqrt(I*V)), name="iSD")
rph	<- mxAlgebra( expression= solve(sqrt(I*V)) %*% V %*% solve(sqrt(I*V)), name="Rph")
rA	<- mxAlgebra( expression= solve(sqrt(I*A)) %*% A %*% solve(sqrt(I*A)), name="Ra" )
rC	<- mxAlgebra( expression= solve(sqrt(I*C)) %*% C %*% solve(sqrt(I*C)), name="Rc" )
rE	<- mxAlgebra( expression= solve(sqrt(I*E)) %*% E %*% solve(sqrt(I*E)), name="Re" )


# Algebra to compute standardized variance components
covP	<- mxAlgebra( expression=A+C+E, name="V" )
StA	<- mxAlgebra( expression=A/V, name="h2")
StC	<- mxAlgebra( expression=C/V, name="c2")
StE	<- mxAlgebra( expression=E/V, name="e2")


#std <- mxAlgebra(expression=cbind(iSD %*% a,iSD %*% c,iSD %*% e),name="std")
stdA <- mxAlgebra(expression=iSD %*% a,name="stda")
stdC <- mxAlgebra(expression=iSD %*% c,name="stdc")
stdE <- mxAlgebra(expression=iSD %*% e,name="stde")

# Constraint on total variance of Ordinal variables (A+C+E=1)
mUnv	<- mxMatrix( type="Unit", nrow=nv, ncol=1, name="Unv" )
varL	<- mxConstraint( expression=diag2vec(V)==Unv, name="VarL" )

# Data objects for Multiple Groups
dataMZ   <- mxData( observed=mzData, type="raw" )
dataDZ   <- mxData( observed=dzData, type="raw" )

# Objective objects for Multiple Groups
objMZ	<- mxExpectationNormal( covariance="expCovMZ", means="M", dimnames=selVars, thresholds="expThres" )
objDZ	<- mxExpectationNormal( covariance="expCovDZ", means="M", dimnames=selVars, thresholds="expThres" )
fitFunction<-mxFitFunctionML()
# Combine Groups

pars1		<-list(Mean,Tr,Thres,mUnv,varL)
pars2		<-list(pathA, pathC, pathE, covA, covC, covE, covP, StA, StC, StE, matI, rph, rA, rC, rE,invSD,stdA, stdC,stdE )
modelMZ	<- mxModel( pars1, pars2, covMZ, dataMZ, objMZ, name="MZ",fitFunction )
modelDZ	<- mxModel( pars1, pars2, covDZ, dataDZ, objDZ, name="DZ",fitFunction )
minus2ll	<- mxAlgebra( expression=MZ.objective + DZ.objective, name="m2LL" )

obj		<- mxFitFunctionAlgebra( algebra="m2LL" )
Conf1		<- mxCI (c ('h2[1,1]', 'h2[2,2]', 'c2[1,1]', 'c2[2,2]', 'e2[1,1]', 'e2[2,2]','stda','stdc','stde'))
Conf2		<- mxCI (c ('Rph[2,1]', 'Ra[2,1]', 'Rc[2,1]', 'Re[2,1]') )
AceModel	<- mxModel( "ACE", pars2, modelMZ, modelDZ, minus2ll, obj, Conf1, Conf2)


# ------------------------------------------------------------------------------
# 4) RUN AceModel

AceFit   <- mxRun(AceModel, intervals=TRUE)
(AceSumm  <- summary(AceFit))
round(AceFit@output$estimate,4)

(h2	<-mxEval(MZ.A, AceFit))
(c2	<-mxEval(MZ.C, AceFit))
(e2	<-mxEval(MZ.E, AceFit))

AceFit$MZ.V

mxCompare(SatFit,AceFit)
mxCompare(Sub1Fit,AceFit)



# 4) Specify ACE Model, with ONE overall set of Thresholds with Age effects on thresholds
# ------------------------------------------------------------------------------------------

# Fit ACE_noa1 model
# -----------------------------------------------------------------------

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

CE_CE_Model <- omxSetParameters(model=AceFit,labels=c('a_1_1','a_2_1','a_2_2'),free=F,values=0,name="CE_CE")
CE_CE_Model@intervals <- list()
CE_CE_Model <- mxModel(CE_CE_Model, Conf1, mxCI (c ('Rph[2,1]', 'Rc[2,1]', 'Re[2,1]') ))


CE_CE_Fit <- mxRun(CE_CE_Model, intervals=TRUE)
CE_CE_Summ <- summary(CE_CE_Fit,verbose=T)
CE_CE_Summ

mxCompare(AceFit,CE_CE_Fit)


# Fit ACE_noc model
# -----------------------------------------------------------------------

AE_AE_Model <- omxSetParameters(model=AceFit,labels=c('c_1_1','c_2_1','c_2_2'),free=F,values=0,name="AE_AE")
AE_AE_Model@intervals <- list()
AE_AE_Model <- mxModel(AE_AE_Model, Conf1, mxCI (c ('Rph[2,1]', 'Ra[2,1]', 'Re[2,1]') ))

AE_AE_Fit <- mxRun(AE_AE_Model, intervals=TRUE)
AE_AE_Summ <- summary(AE_AE_Fit)
AE_AE_Summ                         

mxCompare(AceFit,AE_AE_Fit)

# Fit ACE_noa22 model
# -----------------------------------------------------------------------


AE_CE_Model <- omxSetParameters(model=AceFit,labels=c('c_1_1','c_2_1','a_2_2'),free=F,values=0,name="AE_CE")
AE_CE_Model@intervals <- list()
AE_CE_Model <- mxModel(AE_CE_Model, Conf1, mxCI (c ('Rph[2,1]', 'Ra[2,1]', 'Re[2,1]') ))

AE_CE_Fit <- mxRun(AE_CE_Model, intervals=TRUE)
AE_CE_Summ <- summary(AE_CE_Fit)
AE_CE_Summ                         

mxCompare(AceFit,AE_CE_Fit)


# Fit ACE_noca21 model
# -----------------------------------------------------------------------

AE_AE2_Model <- omxSetParameters(model=AceFit,labels=c('c_1_1','c_2_1','c_2_2','a_2_1'),free=F,values=0,name="AE_AE2")
AE_AE2_Model@intervals <- list()
AE_AE2_Model <- mxModel(AE_AE2_Model, Conf1, mxCI (c ('Rph[2,1]', 'Ra[2,1]', 'Re[2,1]') ))

AE_AE2_Fit <- mxRun(AE_AE2_Model, intervals=TRUE)
AE_AE2_Summ <- summary(AE_AE2_Fit)
AE_AE2_Summ                         

mxCompare(AceFit,AE_AE2_Fit)


# Fit ACE_noa21 model
# -----------------------------------------------------------------------

ACE_ACE_Model <- omxSetParameters(model=AceFit,labels="a_2_1",free=F,values=0,name="ACE_ACE")

ACE_ACE_Fit <- mxRun(ACE_ACE_Model, intervals=TRUE)
ACE_ACE_Summ <- summary(ACE_ACE_Fit)
ACE_ACE_Summ

mxCompare(AceFit,ACE_ACE_Fit)
#tableFitStatistics(AceFit,ACE_ACE_Fit)

CholeskyACENestedc <- list(CE_CE_Fit,AE_AE_Fit,AE_CE_Fit,AE_AE2_Fit,ACE_ACE_Fit)
mxCompare(AceFit,CholeskyACENestedc)

#mxCompare(AceFit,CholeskyACENesteda,CholeskyACENestedc)
AE_AE3_Model <- omxSetParameters(model=AceFit,labels=c('c_1_1','c_2_1','c_2_2','a_2_1','e_2_1'),free=F,values=0,name="AE_AE3")
AE_AE3_Model@intervals <- list()
AE_AE3_Model <- mxModel(AE_AE3_Model, Conf1, mxCI (c ('Rph[2,1]', 'Ra[2,1]', 'Re[2,1]') ))
AE_AE3_Fit <- mxRun(AE_AE3_Model, intervals=TRUE)
AE_AE3_Summ <- summary(AE_AE3_Fit,verbose=T)
AE_AE3_Summ    

mxCompare(AceFit,AE_AE3_Fit)

CholeskyACENestedc <- list(CE_CE_Fit,AE_AE_Fit,AE_CE_Fit,AE_AE2_Fit,ACE_ACE_Fit,AE_AE3_Fit)
mxCompare(AceFit,CholeskyACENestedc)

sink()