Wed, 09/14/2011 - 18:25

Hi,

I am not sure how to specify thresholds when > 1, I was hoping for some help. I have run the ACE models with one threshold (binary response data) successfully and have tried to incorporate the correct code for thresholds >1 into my script from several different scripts but have had no luck.

It seems that I am not putting the right information into dimnames? My data contains one variable with 6 levels, for each twin. I'm also not sure how the specified matrices combine and any tips on how that occurs would be greatly appreciated.

I'll attach my script for the saturated model (which I'll then adapt for the model with parameters specified).

Thanking you,

Jane

p.s. the script may be a little disorganised - I have been playing around with it

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

ord5ThrModel.R | 1.88 KB |

Hi

With thresholds>1 we usually arrange things with an algebra to ensure that threshold 3 > threshold 2 > threshold 1 during optimization. The usual approach is to estimate a matrix in which the first row contains threshold 1 for each variable and is free and unbounded. Then rows 2 onwards are deviations from the row above, and these are constrained to be strictly positive (lower bound say .001). The algebraic trick is to pre-multiply this matrix with a lower triangular matrix of 1's.

There's an example script at http://www.vipbg.vcu.edu/~vipbg/HGEN619/MultiTwinRawOrdEx.R - the key lines are:

mxMatrix( type="Full", nrow=nth, ncol=ntv, free=TRUE, values=thValues, lbound=thLBound, name="Thre" ),

mxMatrix( type="Lower", nrow=nth, ncol=nth, free=FALSE, values=1, name="Inc" ),

mxAlgebra( expression= Inc %*% Thre, name="ThreInc"),

mxAlgebra( expression= ThreInc, dimnames=list(thRows,selVars), name="expThreDZ" ),

mxFIMLObjective( covariance="expCovDZ", means="expMeanDZ", dimnames=selVars, thresholds="expThreDZ")

- note that the mxAlgebra gets dimnames.

HTH

re the comment in this thread

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

"I'll also add that OpenMx does not require the "lower matrix trick" to keep thresholds strictly increasing."

This matrix lock the rows to the same threshold. But the comment implies that OpenMx is smart enough to keep the thresholds monotonic on its own?

nthresh = 3; nVar = 4, nSib=2

mxMatrix("Full", nthresh, nVar*nSib, free=T, name="expThre",

values = seq(0,1, length = nthresh),

labels = glue("th",(1:nthresh))

)

@labels

[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]

[1,] "th1" "th1" "th1" "th1" "th1" "th1" "th1" "th1"

[2,] "th2" "th2" "th2" "th2" "th2" "th2" "th2" "th2"

[3,] "th3" "th3" "th3" "th3" "th3" "th3" "th3" "th3"

@values

[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]

[1,] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0

[2,] 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5

[3,] 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0

To my knowledge, there is nothing in OpenMx that would automatically prevent thresholds from being out of order *during optimization* if the lower matrix trick was not used. Yes, it is certainly possible to start the matrix with thresholds nicely ordered, as Ryne points out in thread http://openmx.psyc.virginia.edu/thread/978 but that does not guarantee that the optimizer will never try a set of parameter values that violates that initial (and desired) ordering.

It would be difficult to do this - it would require openmx to detect that a matrix was used for the thresholds directly (i.e., not the result of an algebra), to internally apply the lower matrix trick, then undo it when communicating back to the user. Reworking the Hessian to obtain correct standard errors (or CI's if requested) would be necessary also.

I've talked with Tim a few times about this, and he claims that it shouldn't be an issue. Line 210, 250 and several other places in omxFIMLObjective.c calls 'checkIncreasing()' on the threshold matrix prior to integration, which at least checks to make sure the thresholds are in the right order. I don't know what checks this function runs, or what happens when it finds a failure, but that's a question for Tim (Brick), as I think he wrote it.

To my knowledge, we have no test suite scripts for which threshold reversal has been a problem, and I haven't seen a script that leads to this error in OpenMx. I'd be very happy to see one, if only to see what the behavior of the FIML objective.

I'll also add that I think that (provided its not necessary for a particular problem), the lower matrix algebra just adds another layer of abstraction. I don't doubt that there are classic Mx scripts that require this algebra, but I'd rather see it as an advanced technique that's taught only after thresholds are fully understood.

OpenMx does check for non-increasing thresholds during optimization, but the system isn't perfect.

The check is performed in the FIML objective function every time the thresholds are recalculated. In the case of a failure, the fallback strategy is the same one used for the case of an objective function returning an infinite value: OpenMx asks the optimizer to back up and try again. If the optimizer doesn't know what to do, the error that comes back says that the thresholds aren't strictly increasing.

In practice, this policy seems to do quite well at keeping thresholds strictly increasing, and we haven't seen it fail in our test suite, but there's no reason to believe that it always works. If anybody's gotten the error about their thresholds being not strictly increasing in a real optimization, please let us know.

I can think of one way we might be able to let people specify an ordering on their free parameters (and therefore on their thresholds), but I'm not sure how helpful it'd be. If anyone has thoughts on that, let us know.

Hi Ryne,

As a learner with OpenMx I found that consideration of the use of the lower triangular matrix (to provide the cumulative displacement across thresholds) useful. It gave me an idea of how matrices can and are used creatively in OpenMx and a more integrated way of thinking about coding. So I guess I'm saying I would prefer they were not excluded as a form of simplification.

Thanks for the input! Your comment helped me rethink my position. While I still think we should be clearer about the actual need for the algebra (i.e., don't do it with binary data), it is a great teaching tool about the flexibility of OpenMx. We really should run some tests to see in what situations the deviations (algebra) approach performs better than just a flat threshold matrix.

That's my main question: if the optimiser has to back out every time it picks a non-monotonic value, that sounds like a lot of backing out, potentially.

I guess just stripping out the lower multi from an ordinal two-group twin script and comparing 100 run times would tell us if it matters... next time I do something ordinal, I'll try and remember to play

You're right, but I suspect that this is more a theoretical problem than a practical one. Thresholds would have to be fairly close together for this to be an issue, and I'd guess that the high interdependence between the thresholds for any particular variable would keep step sizes both low and coincident and prevent most crossings. In cases where this is an issue, I think the algebra approach will work better, as box constraints like the lbounds used in a deviation matrix are a lot easier to satisfy than inequalities.

If you do experiment with this (and please do!), don't forget to remove the lbounds on the deviations matrix.

thanks for confirming

Thanks the model ran and I am learning :)

J