#
#   Copyright 2007-2010 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.

# -----------------------------------------------------------------------------
# Program: testProd.R  
# Author: Michael Neale
# Date: 2019.03.28
#
# ModelType: Marginal Maximum Likelihood with test of non-normality, allowing for product latent variables
# DataType: Ordinal (Binary)
# Field: None
#
# Purpose: Test for product interactions between latent variables
#       Matrix style model input - Raw data input
#       Based on
# Schmitt et al 2005: 
#   Semi-Nonparametric Methods for Detecting Latent Non-normality: 
#   A Fusion of Latent Trait and Ordered Latent Class Modeling. Multivariate Behavioral Research 2006; 41(4),427-443
# Neale (1998):
# Modeling interaction and nonlinear effects with Mx: A general approach. In Marcoulides, G. and Schumacker, R., editors, Interaction and Non-linear Effects in Structural Equation Modeling, pages 43–61. Lawrence Erlbaum Associates, (1998).
#
# RevisionHistory:
# -----------------------------------------------------------------------------

# 1 Load Libraries & Options
require(OpenMx)
require(statmod)
source("prodLatent.R")
set.seed(778921)

# 2 Prepare Data
mxOption(NULL, "Default Optimizer", "NPSOL")
mxOption(NULL, key="Number of Threads", value= (omxDetectCores() - 1))
nVariables<-10
nFactors<-1
nThresholds<-1
nSubjects<-500

loadings <- matrix(.7,nrow=nVariables,ncol=nFactors)
residuals <- 1 - (loadings * loadings)
sigma <- loadings %*% t(loadings) + vec2diag(residuals)
mu <- matrix(0,nrow=nVariables,ncol=1)

nSim <- 500
# Set up matrices to store results
nPar <- nVariables * (nThresholds + nFactors) - nFactors*(nFactors-1)/2
nParNL <- nPar + (nFactors*(nFactors+1)/2) * nVariables
nStat <- 2 #  Increase if more statistics being held
results   <- matrix(NA, nSim, nPar + nStat)
resultsNL <- matrix(NA, nSim, nParNL + nStat)
# Begin Simulation Loop
# Step 3: simulate multivariate normal data
for (iSim in 1:nSim)
{
	continuousData <- mvtnorm::rmvnorm(n=nSubjects,mu,sigma)
	
	# Step 4: chop continuous variables into ordinal data 
	# with nThresholds+1 approximately equal categories, based on 1st variable
	quants<-quantile(continuousData[,1],  probs = c((1:nThresholds)/(nThresholds+1)))
	ordinalData<-matrix(0,nrow=nSubjects,ncol=nVariables)
	for(i in 1:nVariables)
	{
	ordinalData[,i] <- cut(as.vector(continuousData[,i]),c(-Inf,quants,Inf))
	}
	
	# Step 5: make the ordinal variables into R factors
	ordinalData <- mxFactor(as.data.frame(ordinalData),levels=c(1:(nThresholds+1)))
	
	# Step 6: name the variables
	selvars<-paste("Item",1:nVariables,sep="")
	names(ordinalData)<-selvars
	
	
	# -----------------------------------------------------------------------------
	#  Load the product latent function
	
	fitNL <- prodLatent (ordinalData, nQuadPoints=14, nFactors=1, oblique=FALSE, calcSEs="No", calcHessian="No", extraTries=9, run=F)
	fit <- omxSetParameters(fitNL, labels=paste(rep(paste("prodloading",1:nVariables,sep="_"),1), rep(1, each=nVariables, len=nVariables),sep="_"), free=F, values=0, name="Linear")
	fit  <- mxTryHardOrdinal(fit)
#	fitNL <- omxSetParameters(fitNL, labels=paste(rep(paste("prodloading",1:nVariables,sep="_"),1), rep(1, each=nVariables, len=nVariables),sep="_"), free=T, values=0, name="F1 and F1*F1")
	parameters <- omxGetParameters(fit)
	fitNL <- omxSetParameters(fit, names(parameters), values=omxGetParameters(fit), name="F1 and F1sq")
	fitNL <- omxSetParameters(fitNL, labels=paste(rep(paste("prodloading",1:nVariables,sep="_"),1), rep(1, each=nVariables, len=nVariables),sep="_"), free=T, values=0, name="F1 and F1xF1")
	fitNL  <- mxTryHardOrdinal(fitNL)

	results[iSim,] <- c(fit$output$Minus2LogLikelihood, fit$output$status$code, fit$output$estimate)
	resultsNL[iSim,] <- c(fitNL$output$Minus2LogLikelihood, fitNL$output$status$code,  fitNL$output$estimate)
}
save.image("simProd10item.RData")