You are here

Failure to obtain CI's for a umx model

14 posts / 0 new
Last post
noamm's picture
Offline
Joined: 10/10/2018 - 17:20
Failure to obtain CI's for a umx model

Hi all,

I have used umx for a bivariate sex-limitation model, using the umxSexLim function, which have been a great help. However, I've failed in trying to compute CI's for my estimates, or any other indication of significance. Attached below are all the options I've tried. I would appreciate any help in how to compute any indication of significance.

Thanks!

umxCI(m2, run = "yes") 
Running Scalar with 19 parameters
CSOLNP 64815751 -1.00592 -1.93e-09 Error in runHelper(model, frontendStart, intervals, silent, suppressWarnings,  : 
 BLAS/LAPACK routine 'DGEBAL' gave error code -3
 
umxCI_boot(m2, type = "par.observed")
Error in if (model$data$type == "raw") { : argument is of length zero
 
umxConfint(m2, parm = "all",run = TRUE)
Running Scalar with 19 parameters
CSOLNP 50077024 503487 5.035e+05   Error in runHelper(model, frontendStart, intervals, silent, suppressWarnings,  : 
 BLAS/LAPACK routine 'DGEBAL' gave error code -3
 
 
x <- umxRun(m2, intervals = TRUE)
Running Scalar with 19 parameters
 
summary(x)
Summary of Scalar 
 
free parameters:
            name        matrix row       col      Estimate lbound ubound
1         am_r1c1        top.am   1         1  0.7895632830  1e-04       
2         am_r2c2        top.am   2         2  0.0556255100  1e-04       
3         cm_r1c1        top.cm   1         1  0.0001000004 1e-04!       
4         cm_r2c2        top.cm   2         2  0.5058710146  1e-04       
5         em_r1c1        top.em   1         1  0.6265102953  1e-04       
6         em_r2c2        top.em   2         2  0.7983627259  1e-04       
7         af_r1c1        top.af   1         1  1.2427580390  1e-04       
8         af_r2c2        top.af   2         2  0.0001000354 1e-04!       
9         cf_r1c1        top.cf   1         1  0.0001000000 1e-04!       
10        cf_r2c2        top.cf   2         2  0.1961368709  1e-04       
11        ef_r1c1        top.ef   1         1  0.3057593395  1e-04       
12        ef_r2c2        top.ef   2         2  0.9729683133  1e-04       
13        Ra_r2c1       top.Ram   1         2 -0.2442409485     -1      1
14        Rc_r2c1       top.Rcm   1         2  0.0629173480     -1      1
15        Re_r2c1       top.Rem   1         2  0.0883929800     -1      1
16 x_mean_m top.expMeanGm   1 parents_1 -0.0529413841              
17   y_mean_m top.expMeanGm   1   peers_1  0.1286657693              
18 x_mean_f top.expMeanGf   1 parents_1  0.0457339934              
19   y_mean_f top.expMeanGf   1   peers_1 -0.0864138755              
 
Model Statistics: 
              |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
      Model:             19                   2259               7060.65
  Saturated:             NA                     NA                    NA
Independence:             NA                     NA                    NA
Number of observations/statistics: 618/2278
 
Constraint 'top.Keep_it_Positive' contributes 0 observed statistics. 
Information Criteria: 
     |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
AIC:       2542.650               7098.650                 7099.921
BIC:      -7456.787               7182.754                 7122.432
CFI: NA 
TLI: 1   (also known as NNFI) 
RMSEA:  0  [95% CI (NA, NA)]
Prob(RMSEA <= 0.05): NA
To get additional fit indices, see help(mxRefModels)
timestamp: 2019-06-27 15:12:35 
Wall clock time: 0.2890499 secs 
optimizer:  CSOLNP 
OpenMx version number: 2.12.2 
Need help?  See help(mxSummary) 
 
umxSummary(m2, SE = TRUE)
umxSummarySexLim is a beta feature. Some things are broken. If any desired stats are not presented, let me know what's missing
Scalar -2 × log(Likelihood) = 7060.65
Scalar: 2-variable sex limitation analysis

Genetic Factor Correlations (male (Ram) lower triangle, female (Raf) upper)

x y
x 1.00 -0.24
y -0.24 1.00

Opposite sex A Correlations Rao

x y
x 1.00 -0.24
y -0.24 1.00

C Factor Correlations (male (Rcm) lower triangle, female (Rcf) upper)

x y
x 1.00 0.06
y 0.06 1.00

Opposite sex C Correlations Rco

x y
x 1.00 0.06
y 0.06 1.00

E Factor Correlations (male (Rem) lower triangle, female (Ref) upper)

x y
x 1.00 0.09
y 0.09 1.00

Standardized solution (top.[ACE][mf]Std + R[ac]o matrices). Use std=F for raw variance.

am af cm cf em ef Rao Rco
x 0.56 0.80 . . 0.44 0.20 1 1
y 0.07 0.06 0.36 0.16 0.57 0.78 1 1

note: no plots for umxPlotSexLim as yet.

tbates's picture
Offline
Joined: 07/31/2009 - 14:25
thoughts to help with CIs

The BLAS error code looks like a bug for OpenMx or its CSOLNP component? I've not seen it before

64815751 -1.00592 -1.93e-09 Error in runHelper(model, frontendStart, intervals, silent, suppressWarnings,  : 
BLAS/LAPACK routine 'DGEBAL' gave error code -3

But for profile confidence intervals, which can be tricky to compute, switching optimizer when a CI isn't returned can often help.

umx_set_optimizer("NPSOL") # or
umx_set_optimizer("SLSQP")

For what you are calling "any significance in general", I would be testing if parameters can be dropped or not, with

m2 = umxModify(m1, update = "parametertotest", compare = TRUE).

if you call umxModify(m1) (i.e., don't set any paths to update) it will return all the labels of the free parameters, which can help think through what paths you want to test. plot(m2) will confirm what path it is that you've dropped.

Also, mxSE() now works for models with constraints, and SE-based CIs are much quicker than profile-based.

call

mxSE("parameter_to_check", model)
or
mxSE(all[r,c, model) # row and column to pick out the cell(s) from the algebra-computed matrix that you want CIs on.
jpritikin's picture
Offline
Joined: 05/24/2012 - 00:35
DGEBAL

That's strange. We don't even call DGEBAL directly. We'd need a stack trace to diagnose further.

AdminRobK's picture
Online
Joined: 01/24/2014 - 12:15
two remarks
Also, mxSE() now works for models with constraints, and SE-based CIs are much quicker than profile-based.

Just wanted to point out that OP would need to update OpenMx to v2.13 for that...which might be a good idea anyhow. Also, profile-likelihood intervals have superior theoretical properties relative to Wald-style intervals.

AdminRobK's picture
Online
Joined: 01/24/2014 - 12:15
mxVersion?

What do you get from mxVersion()?

AdminRobK's picture
Online
Joined: 01/24/2014 - 12:15
umx bug?

Incidentally, this,

umxCI_boot(m2, type = "par.observed")
Error in if (model$data$type == "raw") { : argument is of length zero

, looks like a umx bug. In a multigroup model, the "container" MxModel generally doesn't have an MxData object.

tbates's picture
Offline
Joined: 07/31/2009 - 14:25
speleology :-)

FYI umxCI_boot function is from here https://openmx.ssri.psu.edu/thread/2598

Now superseded by OpenMx's built-in boot and summary features. All the CI and SE handling in umx is scheduled for a good clean out once #95 is closed for the version 3.0 umx release this Fall.

noamm's picture
Offline
Joined: 10/10/2018 - 17:20
Some responses

Hi all, and thanks for all the responses.

Some answers:

Regarding

 > mxVersion()
OpenMx version: 2.13.2 [GIT v2.13.2]
R version: R version 3.5.2 (2018-12-20)
Platform: x86_64-apple-darwin15.6.0 
MacOS: 10.14.1
Default optimizer: SLSQP
NPSOL-enabled?: No
OpenMP-enabled?: Yes   

Iv'e tried to use both optimizers but the first one have failed:

> umx_set_optimizer("NPSOL") # or
Error in umx_set_optimizer("NPSOL") : 
  The Optimizer 'NPSOL' is not legal. Legal values (from mxAvailableOptimizers() ) are:'CSOLNP' and 'SLSQP'
> umx_set_optimizer("SLSQP")
>    

However, when using the SLSQP optimizer, the CI's computation have failed:

> umxConfint(m2, parm = "all",run = TRUE)
Running Scalar with 19 parameters
               lbound estimate ubound lbound Code ubound Code                         
am_r1c1            NA    0.790     NA          NA          NA
am_r2c2            NA    0.056     NA          NA          NA
cm_r1c1            NA    0.000     NA          NA          NA
cm_r2c2            NA    0.506     NA          NA          NA
em_r1c1            NA    0.627     NA          NA          NA
em_r2c2            NA    0.798     NA          NA          NA
af_r1c1            NA    1.243     NA          NA          NA
af_r2c2            NA    0.000     NA          NA          NA
cf_r1c1            NA    0.000     NA          NA          NA
cf_r2c2            NA    0.196     NA          NA          NA
ef_r1c1            NA    0.306     NA          NA          NA
ef_r2c2            NA    0.973     NA          NA          NA
Ra_r2c1            NA   -0.244     NA          NA          NA
Rc_r2c1            NA    0.063     NA          NA          NA
Re_r2c1            NA    0.088     NA          NA          NA
x_mean_m     NA   -0.053     NA          NA          NA
y_mean_m       NA    0.129     NA          NA          NA
x_mean_f     NA    0.046     NA          NA          NA
y_mean_f       NA   -0.086     NA          NA          NA    

I have succeeded in computing the SE's for the unstandardised parameters. Is it best to just use them? Do you have any other recommendations?

Thank you again,
Noam.

AdminRobK's picture
Online
Joined: 01/24/2014 - 12:15
OK, so you don't have an

OK, so you don't have an NPSOL-enabled build of OpenMx installed. If you want one, follow the instructions on our "Download" page. NPSOL is closed-source/proprietary, so CRAN won't distribute NPSOL-enabled builds of OpenMx.

However, when using the SLSQP optimizer, the CI's computation have failed:

I don't know the inner workings of umxConfint(). Could you run your model with intervals=TRUE, and then run summary() on it with verbose=TRUE? That will show diagnostics for why none of the confidence limits was acceptable.

 > mxVersion()
OpenMx version: 2.13.2 [GIT v2.13.2]
R version: R version 3.5.2 (2018-12-20)
Platform: x86_64-apple-darwin15.6.0 
MacOS: 10.14.1
Default optimizer: SLSQP
NPSOL-enabled?: No
OpenMP-enabled?: Yes   

Thanks for posting your mxVersion() output. The summary() output in your opening post says it was from OpenMx version 2.12, but you are now running only version 2.13, correct? If so, would you mind trying again with CSOLNP? I want to see if the DGEBAL error occurs with the newest release or not.

noamm's picture
Offline
Joined: 10/10/2018 - 17:20
Thanks again,

Thanks again,

I've tried again with CSOLNP and running the model with  intervals=TRUE and verbose=TRUE:

 > umx_set_optimizer("CSOLNP")
> umxConfint(m2, parm = "all",run = TRUE)
Running Scalar with 19 parameters
               lbound estimate ubound lbound Code ubound Code                           
am_r1c1            NA    0.790     NA          NA          NA
am_r2c2            NA    0.056     NA          NA          NA
cm_r1c1            NA    0.000     NA          NA          NA
cm_r2c2            NA    0.506     NA          NA          NA
em_r1c1            NA    0.627     NA          NA          NA
em_r2c2            NA    0.798     NA          NA          NA
af_r1c1            NA    1.243     NA          NA          NA
af_r2c2            NA    0.000     NA          NA          NA
cf_r1c1            NA    0.000     NA          NA          NA
cf_r2c2            NA    0.196     NA          NA          NA
ef_r1c1            NA    0.306     NA          NA          NA
ef_r2c2            NA    0.973     NA          NA          NA
Ra_r2c1            NA   -0.244     NA          NA          NA
Rc_r2c1            NA    0.063     NA          NA          NA
Re_r2c1            NA    0.088     NA          NA          NA
x_mean_m     NA   -0.053     NA          NA          NA
y_mean_m       NA    0.129     NA          NA          NA
x_mean_f     NA    0.046     NA          NA          NA
y_mean_f       NA   -0.086     NA          NA          NA
> x <- umxRun(m2, intervals = TRUE)
Running Scalar with 19 parameters
> summary(x, verbose = TRUE)
Summary of Scalar 
 
data:
$MZm.data
   x_1          y_1          x_2          y_2       
 Min.   :-3.1845   Min.   :-2.1074   Min.   :-3.7694   Min.   :-2.0170  
 1st Qu.:-0.4972   1st Qu.:-0.8199   1st Qu.:-0.6840   1st Qu.:-0.6743  
 Median : 0.5081   Median : 0.0937   Median : 0.3831   Median : 0.1446  
 Mean   : 0.2552   Mean   : 0.1910   Mean   : 0.1665   Mean   : 0.3132  
 3rd Qu.: 1.2049   3rd Qu.: 1.1513   3rd Qu.: 1.2673   3rd Qu.: 1.4081  
 Max.   : 1.7867   Max.   : 3.4629   Max.   : 1.8530   Max.   : 3.8534  
 NA's   :3         NA's   :5         NA's   :4         NA's   :5        
 
$DZm.data
   x_1          y_1            x_2           y_2        
 Min.   :-4.3267   Min.   :-2.017027   Min.   :-5.65610   Min.   :-1.88303  
 1st Qu.:-0.9082   1st Qu.:-0.871336   1st Qu.:-0.77663   1st Qu.:-0.78236  
 Median :-0.3148   Median :-0.058217   Median : 0.04700   Median :-0.14555  
 Mean   :-0.3185   Mean   : 0.003414   Mean   :-0.09487   Mean   : 0.04082  
 3rd Qu.: 0.4589   3rd Qu.: 0.813996   3rd Qu.: 0.86106   3rd Qu.: 0.78922  
 Max.   : 1.6973   Max.   : 3.109265   Max.   : 1.85297   Max.   : 3.12115  
 NA's   :8         NA's   :13          NA's   :5          NA's   :8         
 
$MZf.data
   x_1          y_1           x_2          y_2        
 Min.   :-2.8581   Min.   :-2.01703   Min.   :-3.3053   Min.   :-1.88912  
 1st Qu.:-0.5959   1st Qu.:-0.88836   1st Qu.:-0.5585   1st Qu.:-0.77906  
 Median : 0.4997   Median :-0.30778   Median : 0.4965   Median :-0.20729  
 Mean   : 0.1566   Mean   :-0.04333   Mean   : 0.1902   Mean   :-0.08421  
 3rd Qu.: 1.0123   3rd Qu.: 0.86064   3rd Qu.: 1.1663   3rd Qu.: 0.70066  
 Max.   : 1.9847   Max.   : 2.80406   Max.   : 1.7287   Max.   : 2.72372  
 NA's   :2         NA's   :3          NA's   :3         NA's   :3         
 
$DZf.data
   x_1           y_1          x_2          y_2        
 Min.   :-3.90926   Min.   :-2.0170   Min.   :-3.5580   Min.   :-1.88912  
 1st Qu.:-0.75487   1st Qu.:-1.1508   1st Qu.:-0.5387   1st Qu.:-0.95107  
 Median : 0.31497   Median :-0.5643   Median : 0.2649   Median :-0.22216  
 Mean   : 0.05894   Mean   :-0.3270   Mean   : 0.1463   Mean   :-0.03517  
 3rd Qu.: 1.05503   3rd Qu.: 0.1913   3rd Qu.: 1.0077   3rd Qu.: 0.66856  
 Max.   : 1.85297   Max.   : 3.2326   Max.   : 1.9847   Max.   : 3.83830  
 NA's   :7          NA's   :9         NA's   :16        NA's   :21        
 
$DZo.data
   x_1           y_1           x_2           y_2         
 Min.   :-3.53320   Min.   :-2.01703   Min.   :-5.03291   Min.   :-2.017027  
 1st Qu.:-0.72178   1st Qu.:-0.74939   1st Qu.:-0.88953   1st Qu.:-0.792964  
 Median : 0.03384   Median : 0.05771   Median : 0.14311   Median :-0.085926  
 Mean   :-0.02784   Mean   : 0.14687   Mean   :-0.01495   Mean   : 0.004626  
 3rd Qu.: 0.78246   3rd Qu.: 0.88383   3rd Qu.: 1.02975   3rd Qu.: 0.753380  
 Max.   : 1.85297   Max.   : 4.42853   Max.   : 1.85297   Max.   : 2.951411  
 NA's   :25         NA's   :27         NA's   :12         NA's   :15         
 
free parameters:
             name        matrix row       col      Estimate  Std.Error A lbound
1         am_r1c1        top.am   1         1  0.7895632830 0.16331118    1e-04
2         am_r2c2        top.am   2         2  0.0556255100 0.26090790    1e-04
3         cm_r1c1        top.cm   1         1  0.0001000004 0.15521944 ! 1e-04!
4         cm_r2c2        top.cm   2         2  0.5058710146 0.21668825    1e-04
5         em_r1c1        top.em   1         1  0.6265102953 0.08777096    1e-04
6         em_r2c2        top.em   2         2  0.7983627259 0.11031195    1e-04
7         af_r1c1        top.af   1         1  1.2427580390 0.20415855    1e-04
8         af_r2c2        top.af   2         2  0.0001000354 0.28206959 ! 1e-04!
9         cf_r1c1        top.cf   1         1  0.0001000000 0.18284594 ! 1e-04!
10        cf_r2c2        top.cf   2         2  0.1961368709 0.15562563    1e-04
11        ef_r1c1        top.ef   1         1  0.3057593395 0.05631552    1e-04
12        ef_r2c2        top.ef   2         2  0.9729683133 0.13836033    1e-04
13        Ra_r2c1       top.Ram   1         2 -0.2442409485 0.11593499       -1
14        Rc_r2c1       top.Rcm   1         2  0.0629173480 0.24414338       -1
15        Re_r2c1       top.Rem   1         2  0.0883929800 0.03553627       -1
16 x_mean_m top.expMeanGm   1 x_1 -0.0529413841 0.05435899         
17   y_mean_m top.expMeanGm   1   y_1  0.1286657693 0.05590730         
18 x_mean_f top.expMeanGf   1 x_1  0.0457339934 0.05631440         
19   y_mean_f top.expMeanGf   1   y_1 -0.0864138755 0.04856195         
   ubound
1        
2        
3        
4        
5        
6        
7        
8        
9        
10       
11       
12       
13      1
14      1
15      1
16       
17       
18       
19       
 
Model Statistics: 
               |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
       Model:             19                   2259               7060.65
   Saturated:             NA                     NA                    NA
Independence:             NA                     NA                    NA
Number of observations/statistics: 618/2278
 
Constraint 'top.Keep_it_Positive' contributes 0 observed statistics. 
 
condition number of the information matrix:  66.96394 
maximum absolute gradient:  3.378536  ( cm_r1c1 )
chi-square:  χ² ( df=NA ) = NA,  p = NA
Information Criteria: 
      |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
AIC:       2542.650               7098.650                 7099.921
BIC:      -7456.787               7182.754                 7122.432
CFI: NA 
TLI: 1   (also known as NNFI) 
RMSEA:  0  [95% CI (NA, NA)]
Prob(RMSEA <= 0.05): NA
To get additional fit indices, see help(mxRefModels)
timestamp: 2019-08-11 13:59:03 
frontend time: 0.1872709 secs 
backend time: 0.179646 secs 
independent submodels time: 1.40667e-05 secs 
cpu time: 0.366931 secs 
Wall clock time: 0.366931 secs 
OpenMx version number: 2.13.2 
Need help?  See help(mxSummary) 

I can't understand from that why the confidence limits were NA.

In addition, I think I've tried to install what is needed for NPSOL, but have failed.

> source('https://openmx.ssri.psu.edu/software/getOpenMx.R')
You are now installing the latest version of OpenMx.Error in install.packages : Updating loaded packages
> install.packages(pkgs = c("OpenMx"), contriburl = contrib.url(repos, type = type), dependencies = NA, verbose = TRUE)
trying URL 'https://openmx.ssri.psu.edu/software/bin/macosx/el-capitan/contrib/3.5/OpenMx_2.12.2.tgz'
Content type 'application/x-gzip' length 31603003 bytes (30.1 MB)
==================================================
downloaded 30.1 MB
 
 
The downloaded binary packages are in
    /var/folders/s6/m437k24d69d4wb_cb0bvrzsr0000gn/T//RtmpfZrLMP/downloaded_packages
> umx_set_optimizer("NPSOL")
Error in umx_set_optimizer("NPSOL") : 
  The Optimizer 'NPSOL' is not legal. Legal values (from mxAvailableOptimizers() ) are:'CSOLNP' and 'SLSQP'

.

I would really appreciate any other suggestions.
Thank you,
Noam.

AdminRobK's picture
Online
Joined: 01/24/2014 - 12:15
some tips
I've tried again with CSOLNP and running the model with intervals=TRUE and verbose=TRUE:
[snip]
I can't understand from that why the confidence limits were NA.

summary() with verbose=TRUE isn't printing the CI details table I was expecting. In fact, it isn't printing anything at all about the CIs. I bet that's a umx issue. Could you try again, but run your model with mxRun() instead? Also, it looks like the DGEBAL error didn't happen with CSOLNP and version 2.13, correct?

In addition, I think I've tried to install what is needed for NPSOL, but have failed.

At minimum, you need to reload the package first. The download instructions actually recommend going a step further and restarting R (that's what I do).

tbates's picture
Offline
Joined: 07/31/2009 - 14:25
CIs

e-mail me the model saved as an Rda and I will have a look at CIs for you.

save(m1, file="m1.rda")

noamm's picture
Offline
Joined: 10/10/2018 - 17:20
Thanks

Thanks, Iv'e emailed you the Rda file. I really appreciate your help.

AdminNeale's picture
Offline
Joined: 03/01/2013 - 14:09
Identification?

There are no MZOS pairs - unsurprising, but this means that only one of Ra or Rc can be estimated. Re cannot be estimated because E doesn’t generate resemblance between twins. So please try fixing Rc and Re to zero.