# Program: threshold.R  

#Aim of this script is to estimate the threshold dividing the ever vs never smokers
# Read data from tab delimited file (two_cat.dat) with 'NA' as missing values
# Variabels: Zyg	twin1	twin2	Age	Sex
# zyg: 1=mz, 2=dz 
# Sex: 1=male, 2=female 

require(OpenMx)
Canabis <- read.table ('two_cat.dat', header=T )

# Print Descriptive Statistics
# ---------------------------------------------------------------------
summary(Canabis$twin1)
table(Canabis$twin1)

# Select data
# ---------------------------------------------------------------------
Canabis1 <-data.frame(Canabis$twin1, Canabis$Age, Canabis$Sex)
names(Canabis1) <- c("twin1", "Age", "Sex")

# Specify and Run Saturated Model (Tetrachoric correlations) with RawData
# -----------------------------------------------------------------------
nvar <- 1
nthresh <- 1
Vars <-('twin1')
selVars <- ('twin1')
defVars <- c('age','sex')

AgeSexRegressionModel<- mxModel("AgeSexRegression",
#Bris<-
mxModel("BrisbaneData",
  # Matrix & Algebra for expected means vector (SND), Thresholds and correlation
	mxMatrix(type="Zero", nrow=nvar, ncol=nvar, name="expMean"),
	mxMatrix(type="Full", nrow=nvar, ncol=nthresh, free=TRUE, values=0,
		 name="ucThresh", label="Uncorrected_Threshold"),
	mxMatrix(type="Stand", nrow=nvar, ncol=nvar, name="expCor"),

#Matrices for age and sex corrections
      	mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels="data.Age", name="Age"),
      	mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels="data.Sex", name="Sex"),
      	mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=0.5, labels="AgeBeta", name="AgeB"),
      	mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=0.5, labels="SexBeta", name="SexB"),
	mxAlgebra( expression= AgeB%*%Age + SexB%*%Sex, name="AgeSexR"),
	mxAlgebra( expression= (ucThresh + AgeSexR), name="expThresh", dimnames=list(NA,selVars) ), 

 	mxData(Canabis1, type="raw"),
 	mxFIMLObjective( covariance="expCor", means="expMean", thresholds="expThresh",  dimnames=selVars )),
 	mxAlgebra(BrisbaneData.objective, name="-2LL" ),
 	mxAlgebraObjective("-2LL") 
 
 ) 
AgeSexRegressionFit <- mxRun(AgeSexRegressionModel)
AgeSexRegressionSumm <- summary(AgeSexRegressionFit)
AgeSexRegressionSumm


# Generate Model Output
# -----------------------------------------------------------------------
ExpThresholds <- mxEval(BrisbaneData.expThresh, AgeSexRegressionFit )
ExpThresholds
ExpTetraCor <-mxEval(BrisbaneData.expCor, AgeSexRegressionFit )
ExpTetraCor
Sex_Age_Betas <-mxEval(c(BrisbaneData.SexB,BrisbaneData.AgeB), AgeSexRegressionFit )
Sex_Age_Betas
minusTwiceLL <- mxEval(objective, AgeSexRegressionFit )
minusTwiceLL

AgeSexRegressionFit@submodels$BrisbaneData@matrices$Age