require(OpenMx)

# ----------------------------------
# Read the data and print descriptive statistics.

data(factorExample1)

# ----------------------------------
# Build an OpenMx single factor covariance model with fixed variance

indicators <- names(factorExample1)
latents <- c("F1")
loadingLabels <- paste("b_", indicators, sep="")
uniqueLabels <- paste("U_", indicators, sep="")
meanLabels <- paste("M_", indicators, sep="")
factorVarLabels <- paste("Var_", latents, sep="")

oneFactorCov1 <- mxModel("Single Factor Covariance Model with Fixed Variance",
												 type="RAM",
												 manifestVars=indicators,
												 latentVars=latents,
												 mxPath(from=latents, to=indicators, 
												 			 #           arrows=1, all=TRUE, 
												 			 arrows=1, connect="all.pairs", 
												 			 free=TRUE, values=.2, 
												 			 labels=loadingLabels),
												 mxPath(from=indicators, 
												 			 arrows=2, 
												 			 free=TRUE, values=.8, 
												 			 labels=uniqueLabels),
												 mxPath(from=latents,
												 			 arrows=2, 
												 			 free=FALSE, values=1, 
												 			 labels=factorVarLabels),
												 mxData(observed=cov(factorExample1), type="cov", numObs=500)
)

oneFactorCov1Out <- mxRun(oneFactorCov1, suppressWarnings=TRUE)

summary(oneFactorCov1Out)

omxCheckError(mxBootstrap(oneFactorCov1Out),
							"MxComputeBootstrap: data 'Single Factor Covariance Model with Fixed Variance.data' of type 'cov' cannot have row weights")

# ----------------------------------
# check for correct values

expectVal <- c(0.68396, 0.32482, 0.10887, 0.47441, 0.60181, 1.12064,
							 1.25934, 0.64739, 0.71873, 0.3528, 0.17619, 0.19354, 0.79988, 0.63306,
							 0.36763, 0.34024, 0.23404, 0.85441)

expectSE <- c(0.035244, 0.022431, 0.020808, 0.044662, 0.042295, 0.045789, 
							0.048931, 0.03064, 0.049368, 0.02492, 0.01197, 0.01234, 0.052165, 
							0.042853, 0.032175, 0.034939, 0.017353, 0.057948)

# cat(deparse(round(oneFactorCov1Out$output$estimate, 5)))
omxCheckCloseEnough(expectVal, oneFactorCov1Out$output$estimate, 0.001)

omxCheckCloseEnough(expectSE, 
										as.vector(oneFactorCov1Out$output[['standardErrors']]), 0.001)

omxCheckCloseEnough(1435.94, oneFactorCov1Out$output$minimum, 0.001)


# ----------------------------------
# Check that WLS can be swapped into a path model just by
#  adding a WLS data set.
oneFactorCovWLS <- mxModel(oneFactorCov1Out, name='WLS',
													 mxDataWLS(factorExample1)
)
oneFactorCovWLS <- omxSetParameters(model=oneFactorCovWLS,labels="b_x1",free=F,values=0)

oneFactorCovWLSOut <- mxRun(oneFactorCovWLS)

mxCompare(oneFactorCov1Out,oneFactorCovWLSOut)