# 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, na.strings="." )

# 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')

checkThresholdModel<- mxModel("checkThreshold",
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="expThresh", label="threshold1",dimnames=list('th1',selVars) ),
	mxMatrix(type="Stand", nrow=nvar, ncol=nvar, name="expCor"), 
 	mxData(Canabis1, type="raw"),

 	mxFIMLObjective( covariance="expCor", means="expMean", dimnames=selVars,
		 thresholds="expThresh" )),
 	mxAlgebra(BrisbaneData.objective, name="-2LL" ),
 	mxAlgebraObjective("-2LL") 
 
 ) 
checkThresholdFit <- mxRun(checkThresholdModel)
checkThresholdSumm <- summary(checkThresholdFit)
checkThresholdSumm


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