#
#   Copyright 2007-2009 The OpenMx Project
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
# 
#        http://www.apache.org/licenses/LICENSE-2.0
# 
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

require(OpenMx)

# Compute polychoric (or tetrachoric) correlation matrix

# Example Data
data <- read.table("mddndzf.dat", na.string=".", 
	col.names=c("t1neur1", "t1mddd4l", "t2neur1", "t2mddd4l"))

polychoricMatrix <- function(data) {
nvar <- dim(data)[[2]]
nthresh <- vector(mode="integer",nvar)

# Make variables into factors
for (i in 1:nvar) {
    nthresh[i] <- length(table(data[,i]))-1
    #nthresh[i] <- nthresh[i] - 0
    data[,i] <- mxFactor(data[,i], c(0:nthresh[i]))
}
maxnthresh <- max(nthresh)

# Populate matrix with threshold deviations, starting threshold 1 at -1 and putting maximum at +1
minthresh <- -1
maxthresh <- 1
thresholdDeviationValues <- matrix(0,nrow=maxnthresh, ncol=nvar)
thresholdDeviationValues[1,] <- minthresh
thresholdDeviationLbounds <- matrix(nrow=maxnthresh, ncol=nvar)
thresholdDeviationLabels <- matrix(nrow=maxnthresh, ncol=nvar)
thresholdDeviationLabels[1,] <- paste("ThresholdDeviation ", 1, 1:nvar)
thresholdDeviationFree <- matrix(F,nrow=maxnthresh, ncol=nvar)
thresholdDeviationFree[1,] <- TRUE
correlationLabels <- matrix(NA,nrow=nvar,ncol=nvar)
for (i in 1:nvar) { 
    if(nthresh[i]>1){
    for (j in 2:nthresh[i]) {
        thresholdDeviationValues[j,i] <- (maxthresh - minthresh) / nthresh[i]
        thresholdDeviationLbounds[j,i] <- .001
        thresholdDeviationLabels[j,i] <- paste("ThresholdDeviation ", j, i)
        thresholdDeviationFree[j,i] <- TRUE
}}
    for (k in 1:nvar) {
        if (i > k) {
        correlationLabels[i,k] <- paste("r",i,k)
        correlationLabels[k,i] <- paste("r",i,k)}
                }}
nameList <- names(data)
tnames <- paste("Threshold",1:maxnthresh,sep='')

# Define the model
model <- mxModel('model')
model <- mxModel(model, mxMatrix("Stand", name = "R", nrow = nvar, ncol = nvar, free=TRUE, labels=correlationLabels, dimnames=list(nameList, nameList)))
model <- mxModel(model, mxMatrix("Zero", name = "M", nrow = 1, ncol = nvar, free=FALSE, dimnames = list(NULL, nameList)))

# Threshold differences:
model <- mxModel(model, mxMatrix("Full", name="thresholdDeviations", ncol = nvar, nrow = maxnthresh, free=thresholdDeviationFree, values=thresholdDeviationValues, 
		lbound=thresholdDeviationLbounds, labels = thresholdDeviationLabels))

# For Multiplication
model <- mxModel(model, mxMatrix("Lower", name="UnitLower", nrow = maxnthresh, ncol = maxnthresh, free=F, values=1))

# Algebra to compute Threshold matrix
model$thresholds <- mxAlgebra(UnitLower %*% thresholdDeviations, dimnames=list(tnames,nameList))

# Define the objective function
objective <- mxFIMLObjective(covariance="R", means="M", thresholds="thresholds")

# Define the observed covariance matrix
dataMatrix <- mxData(data, type='raw')

# Add the objective function and the data to the model
model <- mxModel(model, objective, dataMatrix)

# Run the job
model <- mxRun(model)

# Populate seMatrix for return
seMatrix <- matrix(NA,nvar,nvar)
k<-0
for (i in 1:nvar){
    for (j in i:nvar){
        if(i != j) {
            k <- k+1
            seMatrix[i,j] <- model@output$standardErrors[k]
            seMatrix[j,i] <- model@output$standardErrors[k]
        }
    }
}
# Add dimnames to thresholds, which oddly are not in model$thresholds' output
thresholds <- matrix(model@output$algebras$model.thresholds, nrow=maxnthresh, ncol=nvar, dimnames=list(tnames,nameList))      

# Toss results back to front end
return(list(polychorics=model$R@values, thresholds=thresholds, polychoricStandardErrors=seMatrix,Minus2LogLikelihood=model@output$Minus2LogLikelihood))
}