I've noticed this a while ago but forgot to post it to you developers.
Something goes wrong when you include definition variables in a model where both continuous and categorical variables are included, and use the definition variables for regression on the mean. For example running this code yields on my machine:
library(OpenMx) N <- 2000 u <- rbinom(N,1,.5) x <- .5*u+rnorm(N) y <- mxFactor( rbinom(N,1,pnorm(-2+u)) , levels=c(0,1) ) model <- mxModel( 'BinCont', mxMatrix('Full',nrow=1,ncol=2,free=c(T,T),name='Betas'), mxMatrix('Full',nrow=1,ncol=1,free=F,label='data.u',name='U'), mxMatrix('Full',nrow=1,ncol=2,free=c(T,F),name='Means'), mxAlgebra( Means + Betas%x%U , name='eMean'), mxMatrix('Full',nrow=1,ncol=1,free=T,name='Thresh'), mxMatrix('Symm',nrow=2,ncol=2,free=c(T,T,F),values=c(1,0,1),name='Cov'), mxData( data.frame(x,y,u) , type='raw'), mxFIMLObjective( means='eMean', covariance='Cov',thresholds='Thresh',threshnames='y',dimnames=c('x','y') ) ) mxRun( model )@output$estimate[1:2] > mxRun( model )@output$estimate[1:2] Running BinCont BinCont.Betas[1,1] BinCont.Betas[1,2] 0.1 0.1
So nothing has ahppened with the regression coefficients. This is true on my Windows machine, as well as on my Linux server. To make things even more complicated, this happen when changing the number of thresds used for fitting model:
> mxOption( NULL , 'Number of Threads' , 2) > mxRun( model )@output$estimate[1:2] Running BinCont BinCont.Betas[1,1] BinCont.Betas[1,2] 0.5015814 0.9228213 > mxOption( NULL , 'Number of Threads' , 3) > mxRun( model )@output$estimate[1:2] Running BinCont BinCont.Betas[1,1] BinCont.Betas[1,2] 1.1945377 0.5031962 > mxOption( NULL , 'Number of Threads' , 4) > mxRun( model )@output$estimate[1:2] Running BinCont BinCont.Betas[1,1] BinCont.Betas[1,2] 0.5015814 0.9228213
I.e. the result seems to be correct when stating the use of an even number of threads, but not when using an odd number of threads! To make things even more confusing, this is true on the Linux server (where the number of threads actually is changed) as well as on Windows (where only one thread is used).
Am I doing something wrong? Is the use of FIMLObjective causing something of this?
Printer-friendly version
I ran the same code on a Mac (10.8.4, quad core i7, R 2.15.3, OpenMx 999.0.0-2883) and got the correct answer with 1, 2 and 4 threads, but unrealistic betas for 3 threads. I have no idea what's causing this.
This model is also reporting indepndence and saturated M2LLs of -1.9xxx and -2, when the model M2LL is around 6700. The "bad" solution with 3 threads has a lower M2LL, which is weird. All models status 0: 1, 2 & 4 have 13 npsol iterations and 145 evaluations, while 3 has 11 and 138. If I rerun the model with model 3's start values (model1a <- res3), then 2 and 4 return to their previous solutions, where 1 thread gets stuck at a local solution with a M2Ll of 6900.
I compared the raw likelihoods of model 2 and model 3, and did nothing productive beyond make a cool graph.
Note on models below. res1 has the run model with 1 thread, res2 has 2 threads, etc. As 3 threads had the problem, I reran each model with res3's final values, and called everything res1a through res4a.
Confirm this as a bug. With repository version updated as of today, 999.0.0-2914 I get betas stuck at starting values with 1 thread, and inconsistent results with 2-6 threads
@values [,1] [,2] [1,] 0.5973846 0.9912441 [,1] [,2] [1,] 1.131837 0.6336549 [,1] [,2] [1,] 1.415889 0.4357451 [,1] [,2] [1,] 0.9577509 0.6713813 [,1] [,2] [1,] 1.169111 0.5255151Here's code to work on
Please try to reproduce with SVN 3717.
I now get consistent results with R 3.1.1 and rev3171 under OSX Mavericks. Nice going Joshua! Do you think this change will affect other scripts with ordinal data?
Joint continuous-ordinal was completely broken with respect to definition variables. Pure ordinal or pure continuous models are not affected.
Major congratulations on finding and fixing this bug, Josh! A joint ordinal continuous script that had been giving weird SE's and which had parameter estimates of regressions-on-definition-variables stuck at the starting values is now giving very sensible results. Great job!