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

# Estimate threshold from 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")
exp     	<- 		mxExpectationNormal( covariance="expCor", means="expMean", dimnames=selVars,
				thresholds="expThresh" )
funML     	<- 	mxFitFunctionML()

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

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