I get the following error, which is one of a series of virtually incomprehensible ones encountered trying to get this script to work.

Error in as.character(x) :

cannot coerce type 'closure' to vector of type 'character'

This problem is designed to fit a 2-class latent class model to binary data, using a mixture distribution.

Attachment | Size |
---|---|

lazarsfeld.txt | 175 bytes |

lazarsfeld.R | 3.25 KB |

I haven't run your code, but I see two problems, either one of which could cause a weird error.

First, I think you should delete "Class1 <-" and "Class2 <-". When you try to assign those objects, you're currently in an mxModel function call.

Second, change "mxRun(lca)" to "mxRun(lcamodel)". You're using the OpenMx namespace when you should be using the R namespace.

Hi Mike,

I fixed up the model. It doesn't run yet, the current problem is that the threshold matrices need column names. The script contained a few typos of minor importance. In mxMatrix(), the argument is called "labels" instead of "names". The function is mxBounds() instead of mxBound(). The function mxConstraint() requires the names of algebras for the left-hand and right-hand sides. The script had the expressions themselves inside quote strings. Just remember the rule: algebra expressions can never appear inside quotes.

But really the mysterious error was caused by the fact that your model was named "lca" and one of your algebras was also named "lca". This caused the library to replace the entire model with a single algebra, and then subsequently it tried to insert the objective function inside an algebra object instead of a model object and then the badness happens.

I'll add a bugfix after I no longer have a fever. A separate bug that your script happened to trigger. It turns out that mxAlgebra(1, name = 'foo') causes an error. Just like mxAlgebra((A), name = 'foo') causes an error for similar reasons. I'll fix that bug too when I feel better.

Oh sorry. I forgot to mention that the "Iden" matrix type is implemented to forbid any free parameters. So does "Unit" and "Zero". Maybe that's too strict, and we should allow free parameters for these types in the mxMatrix() constructor? OTOH, I would expect to see a forum post from some of our beta testers asking us why we allowed them to accidentally set free parameters on these types. I have no opinion on this, I'll do whatever is the consensus on the matter.

Umm, wait. Don't we want the expected covariance matrix to be symmetric? Replace:

mxMatrix("Full", name = "R", values = mxEval(R, class1), free=TRUE)

with:

mxMatrix("Symm", name = "R", values = mxEval(R, class1), free=TRUE)

I think I implemented so that it will generate exactly 10 free parameters if nvar = 4.

IMHO, Iden Unit and Zero should be free=F

As a shortcut for filling the values on free matrices, perhaps recognising matrix types as key words for values and free would be good. So then this would be legal:

But I am guessing others might be averse to that.

Nasty naming effect!

I try and prefix names with a type marker, like modelLCA and algLCA to avoid that.

Hope you're feeling better soon - rotten viruses trying to run their code on our machinery not good...

Woohoo! It now works. I have uploaded, svn 1001, in trunk/models/passing/LCAlazarsfeld.R with Mx1 script & output in trunk/models/passing/LCAlazarsfeld.mx and trunk/models/passing/LCAlazarsfeld.mxo

In general, it was a bit of a fight to get this to work, and this is something I am experiencing generally with OpenMx. Granted, this may be due to my incompetence in R. However, I can usually work my way around error messages if they are reasonably clear, and even when they're a bit cryptic. Seems to me that we need a REALLY GOOD LOOK AT THE ERROR MESSAGES from OpenMx. Otherwise, there is a danger that it will become the province of the expert programmer only, and will not be used by the general clinician/researcher we wish to form a large part of the user base.

Sorry r1002 may not work on your platform. I didn't notice that ThresholdClass1 and ThresholdClass2 are interchangeable, with respect to correctness. This has been fixed in r1005.

echo mike's comment about the esoteric nature of error/warning messages. here is one that i just got:

In model 'thisOMxModel' NPSOL returned a non-zero status code 6. 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)

i've been doing minimization for a while now and i have never encountered the phrase "first-order optimality conditions." (please tell me what the "second order" and "third order" optimality conditions are.) there are two "optimality conditions:" (1) all gradient elements very close to 0 and (2) all eigenvalues of the Hessian (NOT the updated Hessian) are greater than 0. this is NOT my opinion--it is basic calculus.

imho, if the line search is flat, then there is a problem with either the updated hessian or step-size restrictions in the linesearch. in either case, a "clinician" is likely to be completely befuddled. given the vast literature on human factors research, OpenMx can be adjusted to avoid such circumstances.

Ok, optimization issues can borrow Mx1's error messages. For IFAIL = 6 we had:

********************************************

I am not sure I have found a solution that satisfies

Kuhn-Tucker conditions for a minimum.

NAG's IFAIL parameter is 6

Looks like I got stuck here. Check the following:

1. The model is correctly specified

2. Starting values are good

3. You are not already at the solution

The error can arise if the Hessian is ill-conditioned

You can try resetting it to an identity matrix and fit

from the solution by putting TH=-n on the OU line

where n is the number of refits that you want to do

If all else fails try putting NAG=30 on the OU line

and examine the file NAGDUMP.OUT and the NAG manual

********************************************

Obviously, the above needs to be changed because there is no option th in OpenMx, afaik.

It seems like there isn't a full list of mxOptions in the documentation, so that at least a link to the NPSOL manual, which can be found on the manuals link on the page http://www.sbsi-sol-optimize.com/asp/sol_product_npsol.htm would be a start.

I would modify the text of the IFAIL = 6 result to read:

********************************************

Optimization may not have been successful, NPSOL IFAIL= 6. For a full technical description of what this means, please see http://www.sbsi-sol-optimize.com/asp/sol_product_npsol.htm and read the notes below.

There are several things to note about IFAIL = 6

1. Often, but not always, IFAIL = 6 is returned when optimization has completed successfully.

2. It is still the responsibility of the user to check that the solution is a global minimum. Unfortunately, optimization is not an exact science.

3. Checking that the solution is a global minimum can be done several ways:

i) Fit the model again from the present solution several times. If the fit improves, and keeps improving over the successive optimizations then a global solution probably has not been found, and the number of refits should be increased until the solution stabilizes. Here is some example R code that can be used to refit the model:

refits <- 3

for(i in 1:refits) {

model <- mxOption(model, "Cold Start")

model <- mxRun(model)

summary(model)

}

ii) Try starting optimization from a different place by changing the starting values of the parameters. These could be manually changed, but note that changing all the starting values from say .6 to .7 may not be sufficient. Random starting values could be assigned using, e.g., rnorm(), but be careful to keep the range of possible starting values reasonable.

iii) Consider whether the model is identified. IFAIL = 6 does occur more often with underidentified, or nearly-identified models. If two different sets of parameter estimates yield the same -2lnL (or whichever fit function is being used) then the model can be said to be underidentified. Sometimes this condition arises because of characteristics of the data, particularly if the dataset does not provide much information to estimate every parameter of the model.

iv) Examine technical output of optimization returned by NPSOL. At present this is not captured by OpenMx and this will be corrected soon. Under Mx1 it was optionally stored in the file NAGDUMP.OUT

********************************************

And I welcome input as to how to make this text most easily understood

thanks, mike, but i think you misunderstood my meaning. i had solved the problem within a few seconds of getting the message, several days before i posted. the point i was trying to make was that the messages--given that they are taken literally from the NPSOL value of "inform"--make sense to the illuminati in numerical methods, but will not make sense to casual users who use OpenMx as a tool. how many clinicians who google "first-order optimality" will be able to understand http://plato.asu.edu/papers/paper94/node3.html?

here are two suggestions:

(1) for the relevant values of inform, calculate the four conditions listed under inform=0 in the NPSOL manual. if all four conditions hold, then say that the model has converged.

(2) i always look at: (a) max(abs(gradient element)); (b) number of gradient elements with an absolute value greater than .001; (c) number of gradient elements exactly equal to 0--useful in isolating an unidentified parameter; and (4) the smallest eigenvalue of a numerical estimate of the hessian (NOT the hessian that NPSOL spits back--that is almost guaranteed to be pos def even if the model has redundant parms).