Mon, 09/21/2009 - 19:16

Hi everyone,

Here's an interesting problem I am running into. I was working on improving convergence rates on a script that was giving me some problems. I finally got better convergence on my laptop after switching to FIML and changing the scaling of coefficients. (I think the variances of x and derivatives I was estimating were differing by several orders of magnitude.) Then it was time to get things working on Xgrid. After a while I noticed that one of my Xgrid computers (my desktop) was almost never converging, while my other computers were converging almost all the time.

There were two differences that I could find: my desktop had R 2.9.1 (other computer 2.9.2) and my desktop had a new version of OpenMx (other computer was updated 2 weeks ago). I updated R on the desktop, no improvement. I updated OpenMx on the other machine --- almost none of the models converge now on both machines.

Something has changed with OpenMx in the new version. Is there something new that I might not be specifying in a function which was not previously required? Could there be something different with mxFIMLObjective? Other ideas?

Question unrelated to the previous topic, but related to changes in the new version:

I am running a large number of models using a for loop. Before (v0.1.3) when running multiple models, R didn't tend to give errors that would cause it halt execution. Consequently, I could run through hundreds of models without problem. (Naturally I am storing the status codes so that I can check that models converged.) Now I get the following: "Error: The job for model 'LDE_Model_1' exited abnormally with the error message: Non-positive-definite at free parameters: 0.793 -1.439 0.072 0.054...." This error terminates the current loop, canceling the fitting of numerous additional models. Is there any way for OpenMx to register this in the status codes, rather than halting execution of the program?

You could use result <- try(mxRun(model)). Keep in mind that if an error occurs, you still can't get the result of the optimization under the current system. If an error occurs, then the value stored by result is an invisible object of class '"try-error"' containing the error message.

Somehow I never ran into try(). I'll give a try. Thanks!

Here's some sample code. Sorry about the .R extensions on everything, the forum would not let me use .dat. Run Test.R, it will load information from the other two files.

server: 2 code 4's, 23 code 6's.

desktop 2: 2 code 4's, 23 code 6's.

laptop: 23 code 0's, 2 code 1's.

desktop 1: 23 code 0's, 2 code 1's.

All are running R 2.9.2. The server and desktop 2 are using v0.1.4 of OpenMx. The laptop and desktop 1 are using v0.1.3 of OpenMx.

Did we change the effects of byrow=TRUE?

Try removing byrow=TRUE from the line about matrix L (line 42, I think).

This fixed it on my machine.

Looks good initially --- I'll run a larger test tomorrow.

For those who might be interested: The values of matrix L were set by a E x 3 matrix of the same name. It looks like the L matrix was correctly specified in v0.1.3 regardless of whether the "byrow" statement was included. (Even setting it to FALSE resulted in the same matrix being produced.) In v0.1.4 the inclusion/removal of the byrow statement makes a difference.

The test went well, aside from not getting to it until Friday. The change in the L matrix more than doubled the number of converging models. Thanks!

Poop. There is a bug in 0.1.4-827 and the bug is that when 'values', 'lbound', or 'ubound' are specified as matrix arguments to the mxMatrix() function AND when byrow = TRUE then the matrices are read by columns. Let me explain how on earth this happened. A few weeks ago I noticed one of our scripts had the call mxMatrix(values = 1:9) and this causing OpenMx to report an error because 1:9 are of integer type and values must be of float type. So I added the line to mxMatrix: values <- as.numeric(values). And this line is AFTER I test that is.numeric(values) returns TRUE. * Quirk #1 of the example: is.numeric() returns true on integers or floating point numbers, and as.numeric() always returns floating point numbers. What I did not realise is that as.numeric(values) will convert a matrix into a vector (!?!) Furthermore it converts the matrix into a vector in column-major order, such that when byrow = TRUE the column-major vector is read in row-major order and the badness happens.

The fix has been checked into subversion, along with a test case to make sure this doesn't happen again. We'll release 0.1.5 with this fix on Friday.

Glad to hear this one was found. Nice catch.

Ah - that clarifies my observation (written below).

I never did like that as.numeric() would change matrices into vectors.

so that's bad for me here too

SVN Mx, "R version 2.9.2 (2009-08-24)"

Earlier you mentioned downgrading. Have you tried that? Does it fix the problem?

If the symptoms help: The laptop (OpenMx 0.1.3) with the weaker processor has no problem outpacing the desktop (0.1.4). I just ran 10 models, the desktop produced 10 code-reds, the laptop produces 10 zeros (i.e. all was fine).

So you downgraded the laptop again (you mentioned you tried it on the new version and got failures to converge)?

What happens when you downgrade the desktop?

Do you have the script to share with a url for the data?

This is based on the script in the following forum:

http://openmx.psyc.virginia.edu/thread/144

I haven't tried that script, but it should produce similar problems.

hi,

Running fine here: returns in about 1 second, no errors with openmx from SVN and R 2.9.1. Updated to 2.9.2, and still molto bene.

ldeModel1Fit <- mxRun(ldeModel1) A 3 1 0.041294837 20.0094740 A 3 2 -1.206272186 4.8355851 S 1 1 0.134677338 3.2759996 S 1 2 -0.231808468 6.1104767 S 2 2 0.000001000 4.4628945 S 3 3 0.000001000 3.2374304 U 1 1 0.018200208 0.1093465 U 2 2 0.003239532 0.1314263 U 3 3 0.016667088 6.8067205 U 4 4 0.129030656 0.7500015 U 5 5 0.149651825 3.4870282 U 6 6 0.116827792 2.2173529 U 7 7 0.145345958 1.7590130

Running LDE_Model_1

name matrix row col parameter estimate error estimate

1

2

3

4

5

6

7

8

9

10

11

12

13

Observed statistics: 29

Estimated parameters: 13

Degrees of freedom: 16

-2 log likelihood: -1204.045

Saturated -2 log likelihood: -1765.218

Chi-Square: 561.1724

p: 3.871965e-109

AIC (Mx): 529.1724

BIC (Mx): 244.2398

adjusted BIC:

RMSEA: 0.6020645

running under 2.9.2

Observed statistics: 29

Estimated parameters: 13

Degrees of freedom: 16

-2 log likelihood: -1204.045

Saturated -2 log likelihood: -1765.218

Chi-Square: 561.1724

p: 3.871965e-109

AIC (Mx): 529.1724

BIC (Mx): 244.2398

RMSEA: 0.6020645

I'll post a script that recreates the problem later today. (First I'm off to prep & teach a class.)

I hadn't gotten around to updating the laptop. I've been using source() to get the updates, so I'm not sure how I would go about downgrading the desktop.

The code includes a lot of other pieces, but here are the pieces related to the mxModel:

ldeModel <- mxModel("LDE_Model_1",

mxMatrix("Full", 3, 3,

values=c( 0, 0, 0,

0, 0, 0,

-.01,0, 0),

free=c(FALSE,FALSE,FALSE,

FALSE,FALSE,FALSE,

TRUE, TRUE,FALSE),

name="A", byrow=TRUE),

mxMatrix("Symm", 3, 3,

values=c(.1,0, 0,

0,.1, 0,

0, 0,.1),

free=c( TRUE, TRUE,FALSE,

TRUE, TRUE,FALSE,

FALSE,FALSE, TRUE),

lbound=c(0,NA, NA,

NA,0, NA,

NA,NA,0),

name="S", byrow=TRUE),

mxMatrix("Iden", 3, name="I"),

mxFIMLObjective("R","Means")

)

RunModel <- function(toAnalyze,loadings,manifestVars,type) {

if(type=="Orig") L <- loadings$LDE.Orig

if(type=="GOLD") L <- loadings$LDE.GOLD

len <- length(manifestVars)

fit <- mxRun(mxModel(ldeModel,toAnalyze, mxMatrix("Full",free=FALSE,values=L,name="L",byrow=TRUE), mxMatrix("Diag", len, len, values=.05,free=TRUE, name="U", lbound=0), mxMatrix("Full", free=TRUE,nrow=1, ncol=len,values=0, name="Means", dimnames=list(NULL,manifestVars)), mxAlgebra(L %*% solve(I-A) %*% S %*% t(solve(I-A)) %*% t(L) + U, name="R", dimnames = list(manifestVars,manifestVars))))

Out <- c( fit@matrices$A@values[3,1:2], #eta/zeta

diag(fit@matrices$S@values), #derivative var

fit@matrices$S@values[2,1], #derivative cov

#fit@output$SaturatedLikelihood, #saturated likelihood

#fit@output$Minus2LogLikelihood, #-2LL

fit@output$status[[1]]) #Error Codes

return(Out)

}

Oh, and this:

mxData(Embedded,type="raw", numObs=dim(Embedded)[1])