# 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)
print( "Note no subset command because I want to use all the data")
head(Canabis1)
print( "This won't work because data names cannot contain '.'")
names(Canabis1) <- "twin1"
head(Canabis1)

# Specify and Run Saturated Model (Tetrachoric correlations) with RawData
# -----------------------------------------------------------------------
nvar <- 1
nthresh <- 1
Vars <-('twin1')
selVars <- ('twin1')
Canabis1$twin1 <- mxFactor(Canabis1$twin1, levels=c(0:1) )


Means      <-   mxMatrix(type="Zero", nrow=nvar, ncol=nvar, name="expMean" )
Thresholds <-	mxMatrix(type="Full", nrow=nvar, ncol=nthresh, free=TRUE, values=0,
		 name="expThresh", label="threshold1",dimnames=list('th1',selVars) )
Corelation <-   mxMatrix(type="Stand", nrow=nvar, ncol=nvar, name="expCor")
data 	   <- 	mxData(Canabis1, type="raw")
obj        <-   mxFIMLObjective( covariance="expCor", means="expMean", dimnames=selVars,
		 thresholds="expThresh" )

ThresholdModel <- mxModel("checkThreshold", Means, Thresholds, Corelation, data, obj)

checkThresholdFit <- mxRun(ThresholdModel)
checkThresholdSumm <- summary(checkThresholdFit)
checkThresholdSumm


# Generate Model Output
# -----------------------------------------------------------------------
ExpThresholds <- mxEval(BrisbaneData.expThresh, checkThresholdFit )
ExpThresholds
ExpTetraCor <-mxEval(BrisbaneData.expCor, checkThresholdFit )
ExpTetraCor
minusTwiceLL <- mxEval(objective, checkThresholdFit )
minusTwiceLL