All-continuous data do not include means for weights for the means.

This should add a switch to mxDataWLS to compute the weights in the ordinal/joint way vs in the all-continuous way.

Perhaps call the argument allContinuousMethod. It takes string values 'cumulants' and 'marginal'. It defaults to 'cumulants'.

The Holzinger & Swineford abilities dataset packaged with OpenMx at OpenMx/data/HS.ability.data.txt differs from other sources for this data. Specifically:

1. The numbers of observations for each school differ from other sources. In the file, Grant-White has 156 observations and Pasteur has 145. However, Joreksog (1969) has these dimensions reversed, with Grant-White having 145 observations (p. 193). This difference is confirmed in other sources (p. 318, table 3).
2. The data itself differs from other online sources, specifically it does not agree with the dataset packaged with the MBESS package or that reported on the Mplus website. The MBESS and Mplus datasets match one another.

The dataset discussed is also an identical match for the packaged OpenMx/data/HS.fake.data.txt, with the R help call returning the documentation for HS.ability.data.txt in both cases. I'm not aware how these datasets are intended to relate to one another but this might help clarify why the H&S dataset differs from other sources?

• Variation on cov2cor(), for MxAlgebras, that is "robust" to non-positive variances, either by masking them with zeroes or filtering out the necessary rows and columns.

• Function for MxAlgebras that assembles a covariance matrix from an upper left block, a lower left block, and a lower right block.

The change in -2logL between the point-estimate solution and the solution for a confidence limit is a pretty important piece of information for gauging the validity of a CI, but it doesn't appear in the 'CI details' table. Sure, the user can calculate it him/herself, but the change in fit is important enough that the user shouldn't have to.

During mxTryHard, it would be nice if progress could be kept to a single line that could be refreshed like mxRun console progress output.

As a final report, it might be helpful for mxTryHard to output the variance per parameter of the solutions that it found. For example, if you have a model with 10 parameters where 1 of them is poorly identified then larger variance could help indicate which parameter it is.

If the variance is fixed to zero (RAM's default) then we issue a non-pos-def covariance matrix error. It would be more helpful to complain that ordinary variable X has variance of exactly 0.

The StanHeaders 2.10.0 R package, which will be released on CRAN in a couple of days, moves a header file, which causes R CMD check to fail for OpenMx with

multi_normal_sufficient.hpp:15:54: fatal error: stan/math/prim/scal/meta/VectorViewMvt.hpp: No such file or directory

This can be fixed by changing

#include <stan/math/prim/scal/meta/VectorViewMvt.hpp>

in multi_normal_sufficient.hpp to

#include <stan/math/prim/mat/meta/VectorViewMvt.hpp>

See this forums thread. Apparently, the generic for generating data has no method defined for the MxExpectationMixture class. That means it's not possible to do the bootstrap LRT via mxCompare() for mixture models, which was the whole motivating use case for the bootstrap LRT implementation in the first place!

What an embarrassing oversight...

mxPath currently prints by showing their list of components (see below). It would be more informative to show the paths these will create. SO:

mxPath(letters[1:3], connect="unique.pairs")

Would be helpfully displayed something like:

	# Single arrow paths
	a -> a [label="NA"];
	a -> b [label="NA"];
	b -> b [label="NA"];
	a -> c [label="NA"];
	b -> c [label="NA"];
	c -> c [label="NA"];
	# Covariances...

current print method

mxPath 
$from:  'a', 'b', and 'c' 
$to:  NA 
$arrows:  1 
$values:  NA 
$free:  TRUE 
$labels:  NA 
$lbound:  NA 
$ubound:  NA 
$connect:  unique.pairs 
$joinKey:  NA 

Running a model ("m1") with an error in the algebra (attached as zipped m1.Rda)) model, get the following error:

Error in runHelper(model, frontendStart, intervals, silent, suppressWarnings,  : 
  Non-conformable matrices [(1 x 2) and (1 x 2)] in Matrix Multiply.

Be great to know which matrices... in which models...

m1.rda.zip

I'll start

Handy S3 methods

• logLik to get the log likelihood from a model
• confint for Wald-type confidence intervals on free parameters
• anova for comparing nested models
• coef to extract free parameters
• simulate to generate data from a model

More functions

• mxFitFunctionMultigroup for multiple group models.
• mxFitFunctionMixture for (independent) mixture models
• mxFitFunctionHiddenMarkov for Markov dependent mixture models
• mxFitFunctionWLS & mxDataWLS for weighted least squares estimation
• mxGenerateData to generate data from a model
• mxSE to get standard errors on arbitrary algebraic expressions
• mxGetExpected to extract expected means, covariances, and/or thresholds from models
• mxMI to get modification indicies
• mxTryHard & mxTryHardOrdinal to repeatedly fit a model until it behaves
• mxAutoStart for automatic starting values
• omxDefaultComputePlan to generate a default compute sequence, for editing.
• mxBootstrap for nonparametric bootstrapping to get SEs and CIs
• mxParametricBootstrap for parametric bootstrapping to get marginal p-values for parameters
• mxComputeNelderMead for a general-purpose, derivative-free alternative to the GD optimizers
• imxRobustSE for robust standard errors in non-multilevel FIML analysis of raw data
• mxCompare for model comparison, which is now smart about detecting invalid comparisons, and which can do the bootstrap LRT
• mxCompareMatrix is a variation of mxCompare
• mxFactorScores to estimate regression, ML, or weighted ML factor scores
• mxCheckIdentification to check if a model is locally identified
• mxRefModels for fitting saturated and independence models

The rowDiagnositcs argument to mxFitFunctionML returns a vector of all zeros in OpenMx 2.9.4. The feature was not heavily tested: only checking that the vectors of row likelihoods from two different models were the same in e.g. inst/models/passing/StateSpaceAlg.R. It is also checked that it does not throw an error in the man page for mxFitFunctionML. Below is a script, the last line of which should give the first several row likelihoods, but currently gives all zeros.

library(OpenMx)

# Simulate some data
x=rnorm(1000, mean=0, sd=1)
y= 0.5*x + rnorm(1000, mean=0, sd=1)
tmpFrame <- data.frame(x, y)
tmpNames <- names(tmpFrame)

# Define the matrices
M <- mxMatrix(type = "Full", nrow = 1, ncol = 2, values=c(0,0), 
              free=c(TRUE,TRUE), labels=c("Mx", "My"), name = "M")
S <- mxMatrix(type = "Full", nrow = 2, ncol = 2, values=c(1,0,0,1), 
              free=c(TRUE,FALSE,FALSE,TRUE), labels=c("Vx", NA, NA, "Vy"), name = "S")
A <- mxMatrix(type = "Full", nrow = 2, ncol = 2, values=c(0,1,0,0), 
              free=c(FALSE,TRUE,FALSE,FALSE), labels=c(NA, "b", NA, NA), name = "A")
I <- mxMatrix(type="Iden", nrow=2, ncol=2, name="I")

# Define the expectation
expCov <- mxAlgebra(solve(I-A) %*% S %*% t(solve(I-A)), name="expCov")
expFunction <- mxExpectationNormal(covariance="expCov", means="M", dimnames=tmpNames)

# Choose a fit function
fitFunction <- mxFitFunctionML(rowDiagnostics=TRUE)
# also return row likelihoods, even though the fit function
#  value is still 1x1

# Define the model
tmpModel <- mxModel(model="exampleModel", M, S, A, I, expCov, expFunction, fitFunction, 
                    mxData(observed=tmpFrame, type="raw"))

# Fit the model and print a summary
tmpModelOut <- mxRun(tmpModel)
summary(tmpModelOut)

fitResOnly <- mxEval(fitfunction, tmpModelOut)
attributes(fitResOnly) <- NULL
fitResOnly

# Look at the row likelihoods alone
fitLikeOnly <- attr(mxEval(fitfunction, tmpModelOut), 'likelihoods')
head(fitLikeOnly)

The code is all in R. Should be easy to fix.

Hi OpenMX team,

I am trying to simulate fake data based on an existing dataframe. I have converted all the columns to ordered factors, as required. However, I receive the follow error when trying to create simulated data using the mxGenerateData function. Below is the code. Many thanks in advance for your help.

simulatedData <- mxGenerateData(projectDataSubset_for_simulation,nrow=100)

Calculating asymptotic summary statistics ...
Error: you must specify 'nrow' and 'ncol' arguments in mxMatrix(values = wlsData$thresholds, name = "thresh")

Cheers,

Tim

summary() bug or deep mystery of the universe

http://openmx.ssri.psu.edu/thread/4199
Is the package (OpenMx) associated with MxRAMModel objects not being stored as part of the class definition?

Running a model with different starts for the same parameters leads to this Error

Error: The fixed variable '' has been assigned multiple starting values! See matrix 'A' at location (2, 6) and matrix 'A' at location (1, 6) If you want to randomly select one of these values, call model <- omxAssignFirstParameters(model) before running again.

As omxAssignFirstParameters only touches free parameters, when the variable is FIXED, the error should say:

Error: The fixed variable '' has been assigned multiple starting values! See matrix 'A' at location (2, 6) and matrix 'A' at location (1, 6) Equate the values manually, or free the parameters and call model <- omxAssignFirstParameters(model) before running again.

With StanHeaders 2.17.1, R CMD check throws an error (with g++ and clang++) like

In file included from /opt/r-devel/library/StanHeaders/include/stan/math/mix/mat/functor/finite_diff_grad_hessian.hpp:6:0,
                 from /opt/r-devel/library/StanHeaders/include/stan/math/mix/mat.hpp:13,
                 from omxMLFitFunction.cpp:19:
/opt/r-devel/library/StanHeaders/include/stan/math/mix/mat/functor/hessian.hpp: In instantiation of ‘void stan::math::hessian(const F&, const Eigen::Matrix<double, -1, 1>&, double&, Eigen::Matrix<double, -1, 1>&, Eigen::Matrix<double, -1, -1>&) [with F = multi_normal_deriv]’:
omxMLFitFunction.cpp:230:77:   required from here
/opt/r-devel/library/StanHeaders/include/stan/math/mix/mat/functor/hessian.hpp:51:15: error: no match for call to ‘(const multi_normal_deriv) (const Eigen::Matrix<double, -1, 1>&)’
         fx = f(x);
              ~^~~
omxMLFitFunction.cpp:109:4: note: candidate: template<class T> T multi_normal_deriv::operator()(Eigen::Matrix<Scalar, -1, 1>&) const
  T operator()(Eigen::Matrix<T, Eigen::Dynamic, 1>& x) const {
    ^~~~~~~~
omxMLFitFunction.cpp:109:4: note:   template argument deduction/substitution failed:
In file included from /opt/r-devel/library/StanHeaders/include/stan/math/mix/mat/functor/finite_diff_grad_hessian.hpp:6:0,
                 from /opt/r-devel/library/StanHeaders/include/stan/math/mix/mat.hpp:13,
                 from omxMLFitFunction.cpp:19:
/opt/r-devel/library/StanHeaders/include/stan/math/mix/mat/functor/hessian.hpp:51:15: note:   types ‘Eigen::Matrix<Scalar, -1, 1>’ and ‘const Eigen::Matrix<double, -1, 1>’ have incompatible cv-qualifiers
         fx = f(x);
              ~^~~
/opt/r-devel/etc/Makeconf:168: recipe for target 'omxMLFitFunction.o' failed
make: *** [omxMLFitFunction.o] Error 1
ERROR: compilation failed for package ‘OpenMx’

@bob-carpenter : What changed here?

IIRC, we used to check for NAs in definition variables and halt. But this isn't happening now. We should.

Background: I was getting odd results building a new model with algebras and definition variables. Everything looked OK, but the model was returning start values and code RED).

Turns out the issue is that the definition variables contained NAs, but the model didn't report this... no idea what it did when it hit them, in fact.

Here's a quick model which I think should trigger a stop("definition vars 'data.wt2' can't have missing values...") error.

NOTE: It turns out this toy model doesn't cause evaluation, so isn't a test case. See below for an actual test case.

require(umx)
data(twinData) # ?twinData from Australian twins.
fit1 = mxModel("fit1",
	mxMatrix(name="x", "Full", nrow=1, ncol=1, free=F, values=1),
	mxMatrix(name="def", "Full", nrow=1, ncol=1, free=F, labels= "data.wt2"),
	mxAlgebra(x %*% def, name="prod"),
	mxFitFunctionAlgebra("prod"),
	mxData(twinData[,"wt2", drop=F], type="raw")
)
# ====================
# = This should fail =
# ====================
fit1 = mxRun(fit1) 
# Because...
sum(complete.cases(twinData$wt2))
[1] 86

require(umx); data(twinData) 
# Help optimizer by measuring height in cm
twinData$ht1 = twinData$ht1*100
twinData$ht2 = twinData$ht2*100

mzData = subset(twinData, zygosity == "MZFF")
dzData = subset(twinData, zygosity == "DZFF")
umx_set_optimizer("SLSQP") # preferably NPSOL: CSOLNP needs setup to run this model.
m1 = umxCP(selDVs = c("ht", "wt"), sep = "", dzData = dzData, mzData = mzData)

mxCheckIdentification(m1)

Error in mxCheckIdentification(m1) :
Whoa Nelly. I found an MxConstraint in your model. I just cannot work under these conditions. I will be in my trailer until you reparameterize your model without using mxConstraint().

This is currently a warning, Means requested, but model has no means. Add appropriate TX, TY, KA, and/or KA matrices to get real means.

It should probably be an error, but there is some other code that would need some slight reorganization to make this change.

Can we define the correct independence model for a twin model, and ensure mxRefmodels tests that?

Currently (2.7.12.45) , anova() now accepts a model as input, but outputs an inscrutable error if that's all that's provided

require(umx)
data(twinData) ;  selDVs = c("ht1", "ht2")
mzData <- twinData[twinData$zygosity %in% "MZFF", ]
dzData <- twinData[twinData$zygosity %in% "DZFF", ]
m1 = umxACE(selDVs = selDVs, dzData = dzData, mzData = mzData)

anova(m1)
Error in if (stats[i - 1, "m2ll"] > stats[i, "m2ll"]) { : 
  missing value where TRUE/FALSE needed

mxRefModels fails with univariate models.

# example rda file
load("https://github.com/OpenMx/OpenMx/files/674322/x.txt", verbose = TRUE)

refModels = mxRefModels(x, run = TRUE)

Error in mxData(observed = obsdata, type = datatype, numObs = numObs) :
Observed argument is neither a data frame nor a matrix

mxCompare is the preferred way to see an AIC comparison of models. But mxCompare now warns users to use AIC comparison... Can we remove this warning? Or establish who wants to be warned about this?

mxCompare(m2, m1) # ADE is better
  base comparison ep minus2LL   df      AIC    diffLL diffdf  p
1  ADE       <NA>  4 9658.230 3851 1956.230        NA     NA NA
2  ADE        ACE  4 9659.215 3851 1957.215 0.9849158      0 NA
Warning message:
In collectStatistics(nextBaseSummary, nextCompareSummary) :
  Model 'ADE' has the same degrees of freedom as ACE which means the models are not nested and should not be compared with the likelihood ratio test. Compare these models using the information criteria statistics

mxRename appears to be checking objects that are not models in deciding if a new name is legal.
Here's a model, called ADE, which I want to rename to ACE. Should be legal, IMHO, but isn't:
m1$name
[1] "ADE"
> m2 = mxRename(m1, "ACE")
Error in mxRename(m1, "ACE") :
The name 'ACE' is already used in the model
The issue turns out to be that there's an algebra called ACE in a submodel of ADE.

1. I wouldn't have thought that algebras, matrices etc would be checked.
2. Inconsistently, there's no problem creating a model called "ACE" in the same circumstances. And it runs fine.

It would be helpful to change the mxRename error to include where the first instance of the clash was found. So "The name 'ACE' was already used in the model (m1$top$algebras$ACE)"

BUT, should this even be illegal? If so, should it also be illegal to use the name of an algebra as a model name? I can't see why.

Previously, using the ML fitfunction with the GREML expectation would fit a GREML-type model by ML rather than REML. The GREML expectation knew it was talking to an ML fitfunction, and knew how to present the backend FIML fitfunction with a single row of data, a covariance matrix, and a mean vector (the regression yhats). It was not the recommended way to fit a GREML-type model, but a user who really wanted to do it that way could.

The recent overhaul of the backend fitfunction API was sorely needed and will go a long way toward making the code easier to maintain. But, it changed the behavior so that now using the ML fitfunction with a GREML expectation ends up using the GREML fitfunction in the backend.

Hi,

There appears to be a build error on the stable branch when trying to install this on a linux machine, as follows:

omxFIMLFitFunction.cpp: In member function ‘virtual void omxFIMLFitFunction::invalidateCache()’:
omxFIMLFitFunction.cpp:820:48: error: ‘typeof’ was not declared in this scope
omxFIMLFitFunction.cpp:820:49: error: template argument 1 is invalid
omxFIMLFitFunction.cpp: In member function ‘virtual void omxFIMLFitFunction::init()’:
omxFIMLFitFunction.cpp:1050:56: error: ‘typeof’ was not declared in this scope
omxFIMLFitFunction.cpp:1050:57: error: template argument 1 is invalid
make[1]: *** [omxFIMLFitFunction.o] Error 1
make[1]: Leaving directory `/panfs/panasas01/sscm/zr15964/OpenMx/src'
ERROR: compilation failed for package ‘OpenMx’
* removing ‘/tmp/Rtmpps0oZm/devtools_install_4338205c1876/OpenMx’
Error: Command failed (1)
Execution halted
make: *** [install] Error 1

I used the folllowing commands to try and do this:

git clone https://github.com/OpenMx/OpenMx.git
cd OpenMx/
git checkout stable
make install

Many thanks

Load up the attached RAM model (a supermodel with 2 submodels)

load("https://github.com/OpenMx/OpenMx/files/1704191/eqmeans.rda.txt")

running mxStandardizeRAMpaths yields the following error:

mxStandardizeRAMpaths(eqmeans@submodels[[2]], SE = TRUE)
Error in ParamsCov[paramnames, paramnames] : subscript out of bounds

It works fine if SE is not requested:

mxStandardizeRAMpaths(eqmeans@submodels[[2]], SE = FALSE)

        name             label matrix  row  col  Raw.Value        Raw.SE  Std.Value    Std.SE
 1 DZ.S[1,1] bmi1_with_bmi1_DZ      S bmi1 bmi1 0.88451032 not requested 1.00      not requested
 2 DZ.S[2,1] bmi1_with_bmi2_DZ      S bmi2 bmi1 0.46340859 not requested 0.52      not requested
 3 DZ.S[2,2] bmi2_with_bmi2_DZ      S bmi2 bmi2 0.88395528 not requested 1.00      not requested

Also works fine on the other submodel

mxStandardizeRAMpaths(eqmeans@submodels[[2]], SE = TRUE)

       name             label matrix  row  col  Raw.Value      Raw.SE  Std.Value                     Std.SE
1 MZ.S[1,1] bmi1_with_bmi1_MZ      S bmi1 bmi1 0.88238409 0.020542628 1.000    3.9984e-13
2 MZ.S[2,1] bmi1_with_bmi2_MZ      S bmi2 bmi1 0.46613570 0.016449657 0.528   1.1986e-02
3 MZ.S[2,2] bmi2_with_bmi2_MZ      S bmi2 bmi2 0.88174717 0.020505953 1.000   1.5254e-15

OpenMx version: 2.8.3.71 [GIT v2.8.3-71-g0aa1fc0]
R version: R version 3.4.3 (2017-11-30)
Platform: x86_64-apple-darwin15.6.0
MacOS: 10.13.4
Default optimiser: NPSOL

eqmeans.rda.txt

Allow mxGenerateData to take mxData as input
Would be nice to implement data simulation using mxDataWLS + saturated model + mxGenerateData

CSOLNP can enter a closed loop, such as when running the second and third MxModels in the attached script.
180209--CSOLNP.R.txt

(see also the < 2.8.3 version of this document)

Handy S3 methods

• xx to ...

More functions

• print for mxPath enhanced: mxPaths now displays a compact list of the actual paths created

Add to the CITATION file information for the SEM state space modeling paper

This is a step toward full WLS confidence intervals. First, we need the backend to know if it's full WLS. If it is, then profile CIs are fine.

In the TwinAnalysis_Multivariate_Matrix_Raw.R file,
the current version of the code throws an error.

fit = mxRun(model)
Running ACE with 11 parameters
Warning message:
In model 'ACE' Optimizer returned a non-zero status code 6. The model does not satisfy the first-order optimality conditions to the required accuracy, and no improved point for the merit function could be found during the final linesearch (Mx status RED)

In addition, the estimate no longer returns the correct result

print(ACEest);
#          A1    A2 C1 C2     E1    E2
# ht1   0.938 0.000  0  0  0.345 0.000
# bmi1 -0.030 0.883  0  0 -0.145 0.445

is now

print(ACEest);
#          A1    A2 C1 C2     E1    E2
# ht1   **-0.938** 0.000  0  0  0.345 0.000
# bmi1 -0.030 0.883  0  0 -0.145 0.445

A new StanHeaders is about to be uploaded to CRAN. In src/omxMLFitFunction.cpp , you need to

#include <stan/math/fwd/mat/vectorize/apply_scalar_unary.hpp>

to get R CMD check to pass again. CRAN will probably want a new OpenMx ASAP. Sorry for the non-existent notice but the default MCMC sampler in StanHeaders 2.10 was broken so we have to upload 2.11 now.

Just to store some record of the .version field bug building mxModels.
Up to 2.7.9.48, and with umx 1.5.5 (CRAN as of this post), it can be elicited in RStudio (not sure which version) but also on a Mac running R.app, not RStudio.

Can't elicit the behavior with umx 1.6.0 but also completely unclear why not - not changes made in the umxACE code that I can see. I'm wondering if the bug was about building umx on a version of OpenMx earlier than the one it is now running on. Anyhow will close as of historical interest

data(twinData)
selDVs = c("ht", "wt") # umx will add suffix (in this case "") + "1" or '2'
mzData <- twinData[twinData$zygosity %in% c("MZFF", "MZMM"),]
dzData <- twinData[twinData$zygosity %in% c("DZFF", "DZMM", "DZOS"), ]
m1 = umxACE(selDVs = selDVs, dzData = dzData, mzData = mzData, suffix = ''); 

Error in (function (cl, name, valueClass)  : 
  assignment of an object of class “character” is not valid for @‘.version’ in an object of class “MxModel.ACE”; is(value, "package_version") is not TRUE

mxVersion()
OpenMx version: 2.7.9.48 [GIT v2.7.9-48-gac93eda]
R version: R version 3.3.3 (2017-03-06)
Platform: x86_64-apple-darwin13.4.0 
MacOS: 10.12.5
Default optimiser: CSOLNP

omxParallelCI() provides a way for users to get CIs without redoing the primary optimization to find the point estimates, and without using a custom compute plan. Unfortunately, getting the CI details table in verbose summary output isn't possible with omxParallelCI(). With argument run=FALSE, the user would have to run the returned MxModel object and look at each independent submodel's CI details table, and with run=TRUE, there doesn't seem to be a way to get CI details at all. The problem is that the CI details are stored in the output slot of the MxComputeConfidenceInterval object. With run=FALSE, each independent submodel has its own compute plan containing a CI step, and with run=TRUE, there is no CI step in the compute sequence of the returned MxModel object.

I propose that CI details instead be stored in the MxModel's output slot, and that omxParallelCI() be modified to assemble the CI details table in the returned MxModel's output slot when run=TRUE.

omxParallelCI() was apparently thus named based on the idea that OpenMx would run independent submodels in parallel (via 'snowfall'?), which isn't implemented, and isn't really a priority to be implemented. Perhaps we should create an alias function with a less misleading name--say, omxRunCI()?--and use that as the recommended interface (while of course retaining omxParallelCI() for backward compatibility).

mxRun doesn't seem to check that class(model == "MxModel") so the error is a bit confusing: Says the input was a model, just an old one.

mxRun(1)

in lapply(model@matrices, findSquareBrackets) :
trying to get slot "matrices" from an object of a basic class ("numeric") with no slots
In addition: Warning message:
In mxRun(1) :
You are using OpenMx version 2.7.16 with a model created by an old version of OpenMx. This
may work fine (fingers crossed), but if you run into trouble then please recreate your model with the current version of OpenMx.

twinWLS.R script specifies a model that is not identified without a constraint. WLS gives a terrible error message.

> aceWLS <- mxRun(aceWLS)
Running aceML with 5 parameters
Error in solve.default(t(d) %*% V %*% d) : 
  system is computationally singular: reciprocal condition number = 7.45013e-18
> traceback()
5: solve.default(t(d) %*% V %*% d)
4: solve(t(d) %*% V %*% d) at MxFitFunctionWLS.R#248
3: imxWlsStandardErrors(model) at MxRun.R#257
2: runHelper(model, frontendStart, intervals, silent, suppressWarnings, 
       unsafe, checkpoint, useSocket, onlyFrontend, useOptimizer) at MxRun.R#56
1: mxRun(aceWLS)

Proposed error message: "Standard error matrix is not inveritble. Is your model not identified or in need of a constraint? Try re-running the model with standard errors turned off."

Could make this a warning and just set all SEs to NA.

If the user puts the model first, and expression second, they get an error which suggests some component of the expression was not found in the model.

mxEval(m1, solve(sqrt(I * Vtot)) + ) Error in solve(sqrt(I * Vtot)) : object 'Vtot' not found

Be nice to is(expression, "MxModel") and stop if input 1 is an mxModel, or directly test if input1 is a valid expression.

if(is(expression, "MxModel")){
  stop("Input 1 should be an expression, you gave me a ", class(expression)
}

A model is not a valid expression:

mxEval(m1$top, m1)
Cannot evaluate the object 'top' in the model 'ACE'

The attached script (adapted from
models/passing/IntroSEM-OneFactorCov.R ) fits a common-factor model to
covariance-matrix data by ML, then fits a nested submodel to the same
data by WLS, and compares the two fitted MxModel objects via
mxCompare() . The comparison is of course not valid, but mxCompare()
does not complain.

Similarly, comparisons involving the GREML fitfunction should use -2
times the unrestricted loglikelihood (which is also returned to the
frontend in the fitfunction object), not -2 times the restricted
loglikelihood.

See 170818--bad_mxCompare.txt

It would be useful if mxStandardizeRAMpaths also reported standardized means (M paths). It currently ignores them.

re commit 7a12549

Be nice to tell people how to estimate the hessian when vcov() fails?

if(mxOption(model, "Calculate Hessian")== "No"){
	stop(paste("Can't get vcov, as you have Hessian calculation turned off:",
	    	"Turn on with mxOption(model, 'Calculate Hessian', 'Yes')")
} else {
	stop(paste("Parameter variance covariance matrix is not available:",
	    	"Perhaps the Hessian could not be calculated?")
}

At @tbates suggestion, consider using diagonally weighted least squares in mxAutoStart. This would especially help for problems with large differences in the scales of variables (e.g. inst/models/passing/TwinAnalysis_Multivariate_Matrix_Raw.R ... height and bmi are on very different scales and autostart gives negative variances for this one). The big possible drawback is speed. If DWLS can be used without big differences in speed, then this is a pretty easy change and a good fix.

Love to see weighted least squares in addition to ML!

Hi all,
We are planning to update RcppEigen to include the recent version (3.3.x) of the Eigen library. When checking the reverse dependencies, I found that OpenMx compiled well, but caused a runtime error given below.

* checking examples ... ERROR
Running examples in ‘OpenMx-Ex.R’ failed
The error most likely occurred in:

> ### Name: mxDataWLS
> ### Title: Create MxData Object for Least Squares (WLS, DLS, ULS) Analyses
> ### Aliases: mxDataWLS
> 
> ### ** Examples
> 
> 
> # Create and fit a model using mxMatrix, mxAlgebra, mxExpectationNormal, and mxFitFunctionWLS
> 
> library(OpenMx)
> 
> # Simulate some data
> 
> x=rnorm(1000, mean=0, sd=1)
> y= 0.5*x + rnorm(1000, mean=0, sd=1)
> tmpFrame <- data.frame(x, y)
> tmpNames <- names(tmpFrame)
> wdata <- mxDataWLS(tmpFrame)
Calculating asymptotic summary statistics for 2 continuous and 0 ordinal variables ...
> 
> # Define the matrices
> 
> 
> S <- mxMatrix(type = "Full", nrow = 2, ncol = 2, values=c(1,0,0,1), 
+               free=c(TRUE,FALSE,FALSE,TRUE), labels=c("Vx", NA, NA, "Vy"), name = "S")
> A <- mxMatrix(type = "Full", nrow = 2, ncol = 2, values=c(0,1,0,0), 
+               free=c(FALSE,TRUE,FALSE,FALSE), labels=c(NA, "b", NA, NA), name = "A")
> I <- mxMatrix(type="Iden", nrow=2, ncol=2, name="I")
> 
> # Define the expectation
> 
> expCov <- mxAlgebra(solve(I-A) %*% S %*% t(solve(I-A)), name="expCov")
> expFunction <- mxExpectationNormal(covariance="expCov", dimnames=tmpNames)
> 
> # Choose a fit function
> 
> fitFunction <- mxFitFunctionWLS()
> 
> # Define the model
> 
> tmpModel <- mxModel(model="exampleModel", S, A, I, expCov, expFunction, fitFunction, 
+                     wdata)
> 
> # Fit the model and print a summary
> 
> tmpModelOut <- mxRun(tmpModel)
Running exampleModel with 3 parameters
Error: The following error occurred while evaluating the subexpression 'solve(exampleModel.I - exampleModel.A)' during the evaluation of 'expCov' in model 'exampleModel' : system is computationally singular: reciprocal condition number = 4.08701e-32
Execution halted

It should be a numerical problem with singular matrices, but since I'm not familiar with OpenMx internals, it's pretty hard for me to debug. I wonder if you have any insights on this issue, and a patch for it will help us moving forward on the new version of RcppEigen.

The development version of the RcppEigen package can be retrieved here.

Thank you.

If VC is 1x1 and we ask for the CIs of mxCI( "VC[1,1:3]" ) we get the unhelpful error message e.g.

Error in runHelper(model, frontendStart, intervals, silent, suppressWarnings,  : 
  Requested improper value (1, 2) from (1, 1) matrix 'oneACEvc.VC'

Below is a complete script that reproduces the error

# ----------------------------------------------------------------------------------------------------------------------
# Program: oneACEvc.R  
#  Author: Hermine Maes
#    Date: 01 04 2018 
#
# Twin Univariate ACE model to estimate causes of variation across multiple groups
# Matrix style model - Raw data - Continuous data
# -------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|

# Load Libraries & Options
library(OpenMx)
library(psych); library(polycor)
source("miFunctions.R")
# download this file from
# https://hermine-maes.squarespace.com/#/help/
# use 4openmx as requested

# ----------------------------------------------------------------------------------------------------------------------
# PREPARE DATA

# Load Data
data(twinData)
dim(twinData)
describe(twinData[,1:12], skew=F)

# Select Variables for Analysis
vars      <- 'bmi'                     # list of variables names
nv        <- 1                         # number of variables
ntv       <- nv*2                      # number of total variables
selVars   <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="")

# Select Data for Analysis
mzData    <- subset(twinData, zyg==1, selVars)
dzData    <- subset(twinData, zyg==3, selVars)

# Generate Descriptive Statistics
colMeans(mzData,na.rm=TRUE)
colMeans(dzData,na.rm=TRUE)
cov(mzData,use="complete")
cov(dzData,use="complete")

# Set Starting Values
svMe      <- 20                        # start value for means
svPa      <- .2                        # start value for path coefficient
svPe      <- .5                        # start value for path coefficient for e

# ----------------------------------------------------------------------------------------------------------------------
# PREPARE MODEL

# Create Algebra for expected Mean Matrices
meanG     <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labVars("mean",vars), name="meanG" )

# Create Matrices for Variance Components
covA      <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="VA11", name="VA" ) 
covC      <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="VC11", name="VC" )
covE      <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="VE11", name="VE" )

# Create Algebra for expected Variance/Covariance Matrices in MZ & DZ twins
covP      <- mxAlgebra( expression= VA+VC+VE, name="V" )
covMZ     <- mxAlgebra( expression= VA+VC, name="cMZ" )
covDZ     <- mxAlgebra( expression= 0.5%x%VA+ VC, name="cDZ" )
expCovMZ  <- mxAlgebra( expression= rbind( cbind(V, cMZ), cbind(t(cMZ), V)), name="expCovMZ" )
expCovDZ  <- mxAlgebra( expression= rbind( cbind(V, cDZ), cbind(t(cDZ), V)), name="expCovDZ" )

# Create Data Objects for Multiple Groups
dataMZ    <- mxData( observed=mzData, type="raw" )
dataDZ    <- mxData( observed=dzData, type="raw" )

# Create Expectation Objects for Multiple Groups
expMZ     <- mxExpectationNormal( covariance="expCovMZ", means="meanG", dimnames=selVars )
expDZ     <- mxExpectationNormal( covariance="expCovDZ", means="meanG", dimnames=selVars )
funML     <- mxFitFunctionML()

# Create Model Objects for Multiple Groups
pars      <- list( meanG, covA, covC, covE, covP )
modelMZ   <- mxModel( pars, covMZ, expCovMZ, dataMZ, expMZ, funML, name="MZ" )
modelDZ   <- mxModel( pars, covDZ, expCovDZ, dataDZ, expDZ, funML, name="DZ" )
multi     <- mxFitFunctionMultigroup( c("MZ","DZ") )

# Create Algebra for Variance Components
rowVC     <- rep('VC',nv)
colVC     <- rep(c('VA','VC','VE','SA','SC','SE'),each=nv)
estVC     <- mxAlgebra( expression=cbind(VA,VC,VE,VA/V,VC/V,VE/V), name="VarC", dimnames=list(rowVC,colVC) )

# Create Confidence Interval Objects
ciACE     <- mxCI( "VC[1,1:3]" )  # OFFENDING LINE

# Build Model with Confidence Intervals
modelACE  <- mxModel( "oneACEvc", pars, modelMZ, modelDZ, multi, estVC, ciACE )

# ----------------------------------------------------------------------------------------------------------------------
# RUN MODEL

# Run ACE Model
fitACE    <- mxRun( modelACE, intervals=T ) # LINE THAT GIVES ERROR

Basically all of the umx functions do this, and thus we cannot generate data for them. The mxGetExpected functions handle this. It shouldn't be too hard to fix it for data generation too.

When using SLSQP to optimize an MxModel that contains MxConstraints, it is possible for it to fail to find a local minimum but not warn the user. The reason is that we are using NLOPT's stopping criteria of change in fit less than tolerance; NLOPT does not automatically check the KKT optimality conditions. For unconstrained problems, our backend uses our own finite-difference derivative routine to check the objective's gradient and Hessian, but we do not do this for constrained problems (and correctly so, as the optimality conditions for constrained problems involve more than just the objective's derivatives).