inst/models/passing/ContinuousOnlyWLSTest.R

In OpenMx: Extended Structural Equation Modelling

#
#   Copyright 2007-2021 by the individuals mentioned in the source code history
#
#   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.


#------------------------------------------------------------------------------
# Author: Michael D. Hunter
# Date: 2012.10.02
# Filename: wlsContTest.R
# Purpose: Test the mxDataWLS function for continuous-only data in an
#  OpenMx model.
#------------------------------------------------------------------------------


#------------------------------------------------------------------------------
# Revision history
#  Wed 20 Aug 2014 13:25:04 Central Daylight Time -- Michael Hunter modified to use mxDataWLS.


#--------------------------------------
# Needed packages

library(testthat)
require(OpenMx)
data(Bollen)

suppressWarnings(RNGversion("3.5"))
set.seed(1)
got <- mxGenerateData(Bollen[, 1:8], nrows=10)
omxCheckEquals(nrow(got), 10)

#--------------------------------------
# Set up  model matrices

manvar <- names(Bollen[, 1:8])

lval <- matrix(
		c(1, 0,
		  1, 0,
		  1, 0,
		  1, 0,
		  0, 1,
		  0, 1,
		  0, 1,
		  0, 1),
		byrow=TRUE,
		ncol=2, nrow=8)
lfre <- matrix(as.logical(lval), ncol=2)
lfre[1, 1] <- FALSE
lfre[5, 2] <- FALSE
llab <- matrix(c(paste("lam", 1:4, sep=""), rep(NA, 8), paste("lam", 1:4, sep="")), ncol=2)


lx <- mxMatrix(name="Lam", values=lval, free=lfre, ncol=2, nrow=8, labels=llab, dimnames=list(manvar, c("F1", "F2")))


td <- mxMatrix(name="Theta", type="Symm", ncol=8,
	values=
	c(.8,  0,  0,  0, .2,  0,  0,  0,
	      .8,  0, .2,  0, .2,  0,  0,
	          .8,  0,  0,  0, .2,  0,
	              .8,  0,  0,  0, .2,
	                  .8,  0,  0,  0,
	                      .8,  0, .2,
	                          .8,  0,
	                              .8),
	free=c(T,F,F,F,T,F,F,F,
	         T,F,T,F,T,F,F,
	           T,F,F,F,T,F,
	             T,F,F,F,T,
	               T,F,F,F,
	                 T,F,T,
	                   T,F,
	                     T),
	dimnames=list(manvar, manvar)
)
diag(td$labels) <- paste("var", 1:8, sep="")
selMat <- matrix(
	  c(5,1,
		4,2,
		6,2,
		7,3,
		8,4,
		8,6), ncol=2, byrow=TRUE)
td$labels[selMat] <- paste("cov", c(51, 42, 62, 73, 84, 86), sep="")
td$labels[selMat[,2:1]] <- paste("cov", c(51, 42, 62, 73, 84, 86), sep="")

ph <- mxMatrix(name="Phi", type="Symm", ncol=2, free=T, values=c(.8, .2, .8), labels=paste("phi", c(1, 12, 2), sep=""), dimnames=list(c("F1", "F2"), c("F1", "F2")))


#--------------------------------------
# Set-up WLS model

wlsMod <- mxModel("Test case for WLS Objective function from Bollen 1989",
	lx, ph, td,
	mxExpectationLISREL(LX=lx$name, PH=ph$name, TD=td$name),
	mxFitFunctionWLS(),
	mxData(Bollen[, 1:8], 'raw')
)

dwlsMod <- mxModel(wlsMod, mxFitFunctionWLS("DWLS"))

ulsMod <- mxModel(wlsMod, mxFitFunctionWLS("ULS"))


# Run WLS model
wlsRun <- mxRun(wlsMod)
summary(wlsRun)
omxCheckTrue(is.null(wlsRun$output$calculatedHessian))

expect_equal(length(mxGetExpected(wlsRun, 'standvector')),
             summary(wlsRun)$observedStatistics)

dwlsRun <- mxRun(dwlsMod)

ulsRun <- mxRun(ulsMod)

expect_equal(median(log(diag(wlsRun$data$observedStats$acov))),
             median(log(diag(dwlsRun$data$observedStats$acov))), 1)

expect_equal(diag(wlsRun$data$observedStats$fullWeight),
             diag(dwlsRun$data$observedStats$fullWeight))

print(cbind(coef(wlsRun), coef(dwlsRun), coef(ulsRun)))

omxCheckCloseEnough(vechs(cor(cbind(coef(wlsRun), coef(dwlsRun), coef(ulsRun)))),
                    c(.931, .922, .996), 1e-3)

#TODO Fix summary for WLS data/fitfunctions
# Standard errors correct?
# observed statistics are not 0
# degrees of freedom should be correct when obs stats is fixed
# -2 log likelihood should be "Fit Function Value"?
# AIC, BIC for WLS?


# Compare parameter estimates
# WLS estimated parameters reported in Bollen (1989, p. 428, Table 9.4)
bollenParam <- c(l1=1.00, l2=1.11, l3=1.05, l4=1.16, e1=1.30, e2=7.12, e3=3.53, e4=2.72, e5=2.29, e6=3.77, e7=2.60, e8=3.45)

fitParam <- c(mxEval(Lam, wlsRun)[1:4,1], diag(mxEval(Theta, wlsRun)))

omxCheckCloseEnough(bollenParam, fitParam, epsilon=0.01)

j1 <- omxManifestModelByParameterJacobian(wlsRun, standardize = TRUE)
# Tedious to check all entries, but if variances match then
# other labels are probably correct.
expect_equivalent(j1[1:8, paste0('var',1:8)], diag(8))

#--------------------------------------
# Marginals should be fairly close to cumulants

wlsMO <- omxAugmentDataWithWLSSummary(mxData(Bollen[,1:8], 'raw'), allContinuousMethod = "marginals")
dwlsMO <- omxAugmentDataWithWLSSummary(mxData(Bollen[,1:8], 'raw'), "DWLS", allContinuousMethod = "marginals")

omxCheckCloseEnough(cor(vech(wlsMO$observedStats$cov),
                        vech(wlsRun$data$observedStats$cov)), 1, 5e-3)
omxCheckCloseEnough(cor(vech(wlsMO$observedStats$asymCov[-1:-8,-1:-8]),
                        vech(wlsRun$data$observedStats$asymCov)), 1, .31)
omxCheckCloseEnough(cor(diag(dwlsMO$observedStats$useWeight)[-1:-8],
                        diag(dwlsRun$data$observedStats$useWeight)), 1, .21)

#--------------------------------------
# Re-run with ML
mlMod <- mxModel(wlsMod, name="Rerun WLS model from Bollen as ML",
	mxFitFunctionML(),
	mxData(cov(Bollen[ , 1:8]), "cov", numObs=nrow(Bollen))
)
mlRun <- mxRun(mlMod)


#--------------------------------------
# Compare estimates

mlSum <- summary(mlRun)
wlsSum <- summary(wlsRun)

rms <- function(x, y){sqrt(mean((x-y)^2))}

# parameters are sort of close
expect_equal(rms(mlSum$parameters[,5], wlsSum$parameters[,5]), 0, 0.65)

# standard errors are close
expect_equal(rms(mlSum$parameters[,6], wlsSum$parameters[,6]), 0, 0.18)

# Chi square is on par
expect_equal(mlSum$Chi - wlsSum$Chi, 0, 11)

# Chi square df are the same
omxCheckEquals(mlSum$ChiDoF, wlsSum$ChiDoF)

# RMSEA confidence intervals overlap
omxCheckTrue( (wlsSum$RMSEA < mlSum$RMSEACI[2]) & (wlsSum$RMSEA > mlSum$RMSEACI[1]))
omxCheckTrue( (mlSum$RMSEA < wlsSum$RMSEACI[2]) & (mlSum$RMSEA >= wlsSum$RMSEACI[1]))


#--------------------------------------

wlsMod$data$observed[,1] <- 0.
expect_error(mxRun(wlsMod), "Test case for WLS Objective function from Bollen 1989.data: 'y1' has observed variance less than 1.49012e-08")

wlsMod <- mxModel(wlsMod, mxFitFunctionWLS(allContinuousMethod= 'marginals'))
expect_error(mxRun(wlsMod), "WLS Objective function from Bollen 1989.data: 'y1' has observed variance less than 1.49012e-08")