umxSaturated.R

umxSaturated <- function(m1, evaluate = T) {
	# Use case
	# m1_sat = umxSaturated(m1)
	# summary(m1, SaturatedLikelihood=m1_sat$SaturatedLikelihood, IndependenceLikelihood=m1_sat$IndependenceLikelihood)
	manifests           = m1@manifestVars
	nVar                = length(manifests)
	theData             = m1@data@observed
	dataMeans           = colMeans(theData)
	meansLabels         = paste("mean", 1:nVar, sep="")
	loadingsLabels      = paste("F", 1:nVar, "loading", sep="")
	factorLoadingStarts = t(chol(cov(theData, use = "pairwise.complete.obs")))
	independenceStarts  = diag(cov(theData, use = "pairwise.complete.obs"))

	# Set latents to a new set of 1 per manifest
	# Set S matrix to an Identity matrix (i.e., variance fixed@1)
	# Set A matrix to a Cholesky, manifests by manifests in size, free to be estimated 
	# TODO: start the cholesky at the cov values
	m2 <- mxModel("sat",
    	# variances set at 1,,,
		# mxMatrix(name = "factorVariances", type="Iden" , nrow = nVar, ncol = nVar), # Bunch of Ones on the diagonal
	    mxMatrix(name = "factorMeans"    , type="Zero" , nrow = 1   , ncol = nVar), # Bunch of Zeros
	    mxMatrix(name = "factorLoadings" , type="Lower", nrow = nVar, ncol = nVar, free=T, values = factorLoadingStarts), # labels = loadingsLabels),
	    mxAlgebra(name = "expCov"        , expression = factorLoadings %*% t(factorLoadings)),

	    mxMatrix(name = "expMean", type="Full", nrow = 1, ncol = nVar, values = dataMeans, free = T, labels = meansLabels),
	    mxFIMLObjective(covariance="expCov", means="expMean", dimnames = manifests),
	    mxData(theData, type="raw")
	)
	m3 <- mxModel("independence",
	    # TODO: slightly inefficient, as this has an analytic solution
	    mxMatrix(name = "variableLoadings" , type="Diag", nrow = nVar, ncol = nVar, free=T, values = independenceStarts), # labels = loadingsLabels),
	    mxAlgebra(variableLoadings %*% t(variableLoadings), name = "expCov"),

	    mxMatrix(name = "expMean", type="Full", nrow = 1, ncol = nVar, values = dataMeans, free = T, labels = meansLabels),
	    mxFIMLObjective(covariance="expCov", means="expMean", dimnames = manifests),
	    mxData(theData, type="raw")
	)

	if(evaluate){
		message("I am going to run the saturated and independence models: this may take some time")
		m2 = mxRun(m2)
		m3 = mxRun(m3)
	}
	message("you can use this result in the summary function like this:
	summary(m1, SaturatedLikelihood=m1_sat$SaturatedLikelihood, IndependenceLikelihood=m1_sat$IndependenceLikelihood)
	or use 
	umxReportFit(m1, saturatedModels = m1_sat)")
	return(list(SaturatedLikelihood = m2, IndependenceLikelihood = m3))
}