require(OpenMx)
mydatafile=read.csv("mydataanthro.csv", header=T, na.strings=".")
summary(mydatafile)
mycols=colnames(mydatafile)
VarsA <- 'height1'
VarsB <- 'height2'
Strtpaths <- 0.7
StrtMean <- 0.7
selVars=c(VarsA,VarsB)
print(selVars)
cat("\n")
mzData<-data.frame(subset(mydatafile, Zygosity==0, c(VarsA, VarsB, 'Age1','Age2', 'Gender1', 'Gender2')))
dzData<-data.frame(subset(mydatafile, Zygosity==1, c(VarsA, VarsB, 'Age1','Age2', 'Gender1', 'Gender2')))
datasetname=c(mydatafile,mycols)
datasetname
colnames(mzData) = c(VarsA, VarsB, 'Age1', 'Age2', 'Gender1', 'Gender2')
colnames(dzData) = c(VarsA, VarsB, 'Age1', 'Age2', 'Gender1', 'Gender2')
summary(mzData)
summary(dzData)
colMeans(mzData,na.rm=TRUE)
colMeans(dzData,na.rm=TRUE)
cov(mzData,use="complete")
cov(dzData,use="complete")
MZlist=as.vector(as.matrix(mzData))
MZfinite=MZlist[is.finite(MZlist)]
DZlist=as.vector(as.matrix(dzData))
DZfinite=DZlist[is.finite(DZlist)]
Nobs=length(MZfinite)=length(DZfinite)
print(head(mzData))
print(head(dzData))

#################
#Saturated model
#################
mytwinSatModel <- mxModel("twinSat", 
mxModel("MZ", 
mxMatrix(type="Full", nrow=1,ncol= 2,free=TRUE,values= c(0,0),label= "mean", name="expMean"),
mxMatrix(type="Full", nrow=1, ncol=2, free=TRUE, values= 0.1, label=c("betaGender "), name="bGender"),
mxMatrix(type="Full", nrow=1, ncol=2, free=TRUE, values= 0.1, label=c("betaAge"), name="bAge"),
mxMatrix("Lower", nrow= 2, ncol=2,free=TRUE,values=.5,name="CholMZ"),
mxAlgebra(expression= CholMZ %*% t(CholMZ), name="expCovMZ"), 
mxMatrix(type="Full",nrow=2,ncol=2,free=F,label=c("data.Age1", "data.Age2"), name="Age"), 
mxMatrix(type="Full",nrow=2,ncol=2,free=F,label=c("data.Gender1","data.Gender2"), name="Gender"),
meanAge= mxAlgebra(bAge%*%Age, name= "AgeR"),
meanGender = mxAlgebra (bGender%*%Gender, name= "GenderR"),
mxAlgebra( expression=expMean +AgeR + GenderR , name="expMeanMZ"),
mxData(mzData, type="raw"),
mxExpectationNormal("expCovMZ", "expMeanMZ", dimnames=selVars) ,
mxFitFunctionML()),

mxModel("DZ", 
mxMatrix(type="Full", nrow=1,ncol= 2,free=TRUE,values= c(0,0),label= "mean", name="expMean"),
mxMatrix(type="Full", nrow=1, ncol=2, free=TRUE, values= 0.1, label=c("betaAge"), name="bAge"),
mxMatrix(type="Full", nrow=1, ncol=2, free=TRUE, values= 0.1, label=c("betaGender"), name="bGender"),
mxMatrix("Lower", nrow= 2, ncol=2,free=TRUE,values=.5,name="CholDZ"),
mxAlgebra(expression= CholDZ %*% t(CholDZ), name="expCovDZ"), 
mxMatrix(type="Full",nrow=2,ncol=2,free=F,label=c("data.Age1", "data.Age2"), name="Age"),
mxMatrix(type="Full",nrow=2,ncol=2,free=F,label=c("data.Gender1", "data.Gender2"), name="Gender"),
meanAge= mxAlgebra(bAge%*%Age, name="AgeR"),
meanGender = mxAlgebra (bGender%*%Gender, name="GenderR"),
mxAlgebra( expression=expMean + AgeR + GenderR, name="expMeanDZ"),
mxData(dzData, type="raw"),
mxExpectationNormal("expCovDZ", "expMeanDZ", dimnames=selVars) ,
mxFitFunctionML()),
mxFitFunctionMultigroup( c("MZ.fitfunction",  "DZ.fitfunction")))
mytwinSatFit <- mxRun(mytwinSatModel) #The mxRun command evaluates the model.
LL_Sat <- mxEval(objective, mytwinSatFit)
print(summary(mxRun(mytwinSatModel)))
DF_Sat=Nobs-nrow(mytwinSatFit@output$standardErrors)

################
#Run ACE model
################
#start value for paths (Variance) is (sqrt(variance/3))
twinACEModel <- mxModel("twinACE",
mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=Strtpaths, label="a", name="X" ), 
mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=Strtpaths, label="c", name="Y" ),
mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=Strtpaths, label="e", name="Z" ),
mxAlgebra( expression=X %*% t(X), name="A" ),
mxAlgebra( expression=Y %*% t(Y), name="C" ),
mxAlgebra( expression=Z %*% t(Z), name="E" ),
mxMatrix(type="Full", nrow=1,ncol= 2,free=TRUE,values= StrtMean,label= "mean", name="expMean"),
mxMatrix(type="Full", nrow=1, ncol=2, free=TRUE, values= 0.1, label=c("betaAge"), name="bAge"),
mxMatrix (type="Full", nrow=1, ncol=2, free=TRUE, values= 0.1, label=c("betaGender"), name="bGender"),
mxAlgebra(expression=rbind(cbind(A+C+E , A+C),cbind(A+C, A+C+E)),name="expCovMZ"),
mxAlgebra(expression=rbind(cbind(A+C+E,0.5%x%A+C),cbind(0.5%x%A+C,A+C+E)), name="expCovDZ"),

mxModel("MZ", 
mxData( observed=mzData, type="raw" ),
mxMatrix(type="Full",nrow=2,ncol=2,free=F, label=c("data.Age1", "data.Age2"), name="Age"),
mxMatrix(type="Full",nrow=2,ncol=2,free=F, label=c("data.Gender1", "data.Gender2"), name="Gender"),
meanAge= mxAlgebra(twinACE.bAge%*%Age, name="AgeR"),
meanGender = mxAlgebra (twinACE.bGender%*%Gender, name="GenderR"),
mxAlgebra( expression=twinACE.expMean + GenderR + AgeR, name="expMeanMZ"),
mxExpectationNormal( covariance="twinACE.expCovMZ", means="expMeanMZ", dimnames=selVars),
mxFitFunctionML()),

mxModel("DZ", 
mxData( observed=dzData, type="raw" ),
mxMatrix(type="Full",nrow=2,ncol=2,free=F, label=c("data.Age1", "data.Age2"), name="Age"),
mxMatrix(type="Full",nrow=2,ncol=2,free=F, label=c("data.Gender1", "data.Gender2"), name="Gender"),
meanAge= mxAlgebra(twinACE.bAge%*%Age, name="AgeR"),
meanGender = mxAlgebra (twinACE.bGender%*%Gender, name="GenderR"),
mxAlgebra( expression=twinACE.expMean + GenderR + AgeR, name="expMeanDZ"),
mxExpectationNormal( covariance="twinACE.expCovDZ", means="expMeanDZ", dimnames=selVars),
mxFitFunctionML()),
mxFitFunctionMultigroup( c("MZ.fitfunction",  "DZ.fitfunction")))
twinACEFit <- mxRun(twinACEModel)
print(summary(twinACEFit))

##########################################################
#Summary statistics
##########################################################
DF_ACE=Nobs-nrow(twinACEFit@output$standardErrors)
LL_ACE <- mxEval(objective, twinACEFit)
mychi_ACE= LL_ACE - LL_Sat  
mychi_DF_ACE=DF_ACE-DF_Sat  
mychi_p_ACE=1-pchisq(mychi_ACE,mychi_DF_ACE) 
expMZcov_ACE <- mxEval(expCovMZ, twinACEFit)                                 
expDZcov_ACE <- mxEval(expCovDZ, twinACEFit)                                  
expMeans_ACE <- mxEval(expMean, twinACEFit)     
A_ACE <- mxEval(a*a, twinACEFit)                          
C_ACE <- mxEval(c*c, twinACEFit)                          
E_ACE <- mxEval(e*e, twinACEFit)
V <- (A_ACE+C_ACE+E_ACE)                                
a2_ACE <- A_ACE/V                   
c2_ACE <- C_ACE/V                   
e2_ACE <- E_ACE/V                   
ACE_mySE=round(twinACEFit@output$standardErrors,3)
ACE_myest=round(twinACEFit@output$estimate,3)
ACE_mylower=round(ACE_myest-1.96*ACE_mySE,3)
ACE_myupper=round(ACE_myest+1.96*ACE_mySE,3)

################
#Run AE model
################
twinAEModel <- mxModel(twinACEModel, name = "twinACE",  
mxMatrix(type="Full", nrow=1, ncol=1, free=FALSE, values=0, label="c", name="Y"))
twinAEFit <- mxRun(twinAEModel)

##########################################################
#Summary statistics
##########################################################
LL_AE <- mxEval(objective, twinAEFit)
DF_AE=Nobs-nrow(twinAEFit@output$standardErrors)
mychi_AE= LL_AE - LL_Sat  #subtract LL for Saturated model from LL for AE
mychi_DF_AE=DF_AE-DF_Sat  #subtract DF for Saturated model from DF for AE
mychi_p_AE=1-pchisq(mychi_AE,mychi_DF_AE)# compute chi square probability
#expected values for covs and means can be found in mxEval(expCovMZ, twinAEFit)
A_AE <- mxEval(a*a, twinAEFit)
C_AE <- mxEval(c*c, twinAEFit)
E_AE <- mxEval(e*e, twinAEFit)
V <- (A_AE+C_AE+E_AE)
a2_AE <- A_AE/V
c2_AE <- C_AE/V
e2_AE <- E_AE/V
mychiAEdiff=mychi_AE-mychi_ACE
myDFAEdiff=mychi_DF_AE-mychi_DF_ACE
mychiAEdiff_p=1-pchisq(mychiAEdiff,myDFAEdiff)
AE_mySE=round(twinAEFit@output$standardErrors,3)
AE_myest=round(twinAEFit@output$estimate,3)
AE_mylower=round(AE_myest-1.96*AE_mySE,3)
AE_myupper=round(AE_myest+1.96*AE_mySE,3)
print(summary(mxRun(twinAEModel)))

################
#Run CE model
################
twinCEModel <- mxModel(twinACEModel, name = "twinACE",
mxMatrix(type="Full", nrow=1, ncol=1, free=FALSE, values=0, label="a", name="X"))
twinCEFit <- mxRun(twinCEModel)
DF_CE=Nobs-nrow(twinCEFit@output$standardErrors)
LL_CE <- mxEval(objective, twinCEFit)
mychi_CE= LL_CE - LL_Sat  
mychi_DF_CE=DF_CE-DF_Sat  
mychi_p_CE=1-pchisq(mychi_CE,mychi_DF_CE) 
A_CE <- mxEval(a*a, twinCEFit)
C_CE <- mxEval(c*c, twinCEFit)
E_CE <- mxEval(e*e, twinCEFit)
V <- (A_CE+C_CE+E_CE)
a2_CE <- A_CE/V
c2_CE <- C_CE/V
e2_CE <- E_CE/V
mychiCEdiff=mychi_CE-mychi_ACE
myDFCEdiff=mychi_DF_CE-mychi_DF_ACE
mychiCEdiff_p=1-pchisq(mychiCEdiff,myDFCEdiff)
CE_mySE=round(twinCEFit@output$standardErrors,3)
CE_myest=round(twinCEFit@output$estimate,3)
CE_mylower=round(CE_myest-1.96*CE_mySE,3)
CE_myupper=round(CE_myest+1.96*CE_mySE,3)
print(summary(mxRun(twinCEModel)))

################
#Run E model
################
twinEModel <- mxModel(twinAEModel, name = "twinACE",  mxMatrix(type="Full", nrow=1, ncol=1, free=FALSE, values=0, label="a", name="X"))
twinEFit <- mxRun(twinEModel) 
DF_E=Nobs-nrow(twinEFit@output$standardErrors)
LL_E <- mxEval(objective, twinEFit)
mychi_E= LL_E - LL_Sat  
mychi_DF_E=DF_E-DF_Sat  
mychi_p_E=1-pchisq(mychi_E,mychi_DF_E) 
A_E <- mxEval(a*a, twinEFit)
C_E <- mxEval(c*c, twinEFit)
E_E <- mxEval(e*e, twinEFit)
V <- (A_E+C_E+E_E)
a2_E <- A_E/V
c2_E <- C_E/V
e2_E <- E_E/V
mychiEdiff=mychi_E-mychi_ACE
myDFEdiff=mychi_DF_E-mychi_DF_ACE
mychiEdiff_p=1-pchisq(mychiEdiff,myDFEdiff)
E_mySE=round(twinEFit@output$standardErrors,3)
E_myest=round(twinEFit@output$estimate,3)
E_mylower=round(E_myest-1.96*E_mySE,3)
E_myupper=round(E_myest+1.96*E_mySE,3)
print(summary(mxRun(twinEModel)))
##########################
#Output to compare models
##########################
myoutput <- rbind(cbind("__________________________","_______________","_________________","______________"),

cbind("ACE model","A","C","E"), cbind("Unsquared path estimates", ACE_myest[1],ACE_myest[2],ACE_myest[3]),

cbind("Standard errors",ACE_mySE[1],ACE_mySE[2],ACE_mySE[3]),

cbind("Lower 95% CI",ACE_mylower[1], ACE_mylower[2], ACE_mylower[3]), cbind("Upper 95% CI",ACE_myupper[1], ACE_myupper[2], ACE_myupper[3]), cbind(".",".",".","."),

cbind("Unstandardized variance comps",round(A_ACE,3), round(C_ACE,3), round(E_ACE,3)), cbind("Standardized variance comps",round(a2_ACE,3),round(c2_ACE,3),round(e2_ACE,3)),

cbind(".","chisq","DF","p"), cbind("saturated vs. ACE",round(mychi_ACE,3), mychi_DF_ACE,round(mychi_p_ACE,3)), cbind("__________________________","_______________","_________________","______________"),

cbind("AE model","A","C","E"), cbind("Unsquared path estimates", AE_myest[1],".", AE_myest[2]), cbind("Standard errors",AE_mySE[1],".",AE_mySE[2]), cbind("Lower 95% CI",AE_mylower[1],".",AE_mylower[2]), cbind("Upper 95% CI",AE_myupper[1],".",AE_myupper[2]), cbind(".",".",".","."), cbind("Unstandardized variance comps",round(A_AE,3), round(C_AE,3), round(E_AE,3)), cbind("Standardized variance comps",round(a2_AE,3), round(c2_AE,3),round(e2_AE,3)), cbind(".","chisq","DF","p"),

cbind("saturated vs. AE",round(mychi_AE,3), mychi_DF_AE,round(mychi_p_AE,3)), cbind("ACE vs AE",round (mychiAEdiff,3),1,round(mychiAEdiff_p,3)), cbind(".",".",".","."), cbind ("__________________________","_______________","_________________","______________"),

cbind("CE model","A","C","E"), cbind("Unsquared path estimates",".", CE_myest[1],CE_myest[2]), cbind("Standard errors",".",CE_mySE[1],CE_mySE[2]),cbind("Lower 95% CI",".", CE_mylower[1], CE_mylower[2]), cbind("Upper 95% CI",".", CE_myupper[1], CE_myupper[2]), cbind(".",".",".","."), cbind("Unstandardized variance comps",round(A_CE,3),round(C_CE,3),round(E_CE,3)),

cbind("Standardized variance comps",round (a2_CE,3),round(c2_CE,3),round(e2_CE,3)), cbind(".","chisq","DF","p"),

cbind("saturated vs. CE",round(mychi_CE,3), mychi_DF_CE,round(mychi_p_CE,3)), cbind("ACE vs CE",round(mychiCEdiff,3), 1, round(mychiCEdiff_p,3)), cbind ("__________________________","_______________","_________________","______________"),

cbind("E model","A","C","E"), cbind("Unsquared path estimates",".", E_myest[1],E_myest[2]), cbind("Standard errors",".",E_mySE[1],E_mySE[2]),cbind("Lower 95% CI",".", E_mylower[1], E_mylower[2]), cbind("Upper 95% CI",".", E_myupper[1], E_myupper[2]), cbind(".",".",".","."), cbind("Unstandardized variance comps",round(A_E,3),round(C_E,3),round(E_E,3)),

cbind("Standardized variance comps",round (a2_E,3),round(c2_E,3),round(e2_E,3)), cbind(".","chisq","DF","p"),

cbind("saturated vs. E",round(mychi_E,3), mychi_DF_E,round(mychi_p_E,3)), cbind("ACE vs E",round(mychiEdiff,3), 1, round(mychiEdiff_p,3)), cbind ("__________________________","_______________","_________________","______________"),cbind(date(),".",".","."))

myoutput2 <- data.frame(myoutput)

names(myoutput2)=selVars

#write.table(myoutput2, "myfileout.txt",sep = "\t",quote=F)

write.table(myoutput2, sprintf("%s.txt", VarsA) ,sep = "\t",quote=F)

print(myoutput2)