Failure to obtain CI's for a umx model

Posted on Thu, 07/25/2019 - 07:45

noamm Joined: 10/10/2018

Forums

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.

Replied on Thu, 07/25/2019 - 08:45

tbates Joined: 07/31/2009

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.

Replied on Thu, 07/25/2019 - 09:47

jpritikin Joined: 05/23/2012

DGEBAL

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

Replied on Thu, 07/25/2019 - 10:38

AdminRobK Joined: 01/24/2014

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.

Replied on Thu, 07/25/2019 - 10:23

AdminRobK Joined: 01/24/2014

mxVersion?

What do you get from `mxVersion()`?

Replied on Thu, 07/25/2019 - 10:33

AdminRobK Joined: 01/24/2014

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.

Replied on Fri, 07/26/2019 - 11:44

tbates Joined: 07/31/2009

speleology :-)

FYI `umxCI_boot` function is from here [https://openmx.ssri.psu.edu/thread/2598](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](https://github.com/tbates/umx/issues/95) is closed for the version 3.0 `umx` release this Fall.

Replied on Wed, 08/07/2019 - 08:56

noamm Joined: 10/10/2018

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.

Replied on Wed, 08/07/2019 - 10:10

AdminRobK Joined: 01/24/2014

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.

Replied on Sun, 08/11/2019 - 07:25

noamm Joined: 10/10/2018

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 1         am_r1c1 2         am_r2c2 3         cm_r1c1 4         cm_r2c2 5         em_r1c1 6         em_r2c2 7         af_r1c1 8         af_r2c2 9         cf_r1c1 10        cf_r2c2 11        ef_r1c1 12        ef_r2c2 13        Ra_r2c1 14        Rc_r2c1 15        Re_r2c1 16 x_mean_m top.expMeanGm 17   y_mean_m top.expMeanGm 18 x_mean_f top.expMeanGf 19   y_mean_f top.expMeanGf ubound
 1
 2
 3
 4
 5
 6
 7
 8
 9
 10
 11
 12
 13      1
 14      1
 15      1
 16
 17
 18
 19

matrix row       col      Estimate  Std.Error A lbound
 top.am   1         1  0.7895632830 0.16331118    1e-04
 top.am   2         2  0.0556255100 0.26090790    1e-04
 top.cm   1         1  0.0001000004 0.15521944 ! 1e-04!
 top.cm   2         2  0.5058710146 0.21668825    1e-04
 top.em   1         1  0.6265102953 0.08777096    1e-04
 top.em   2         2  0.7983627259 0.11031195    1e-04
 top.af   1         1  1.2427580390 0.20415855    1e-04
 top.af   2         2  0.0001000354 0.28206959 ! 1e-04!
 top.cf   1         1  0.0001000000 0.18284594 ! 1e-04!
 top.cf   2         2  0.1961368709 0.15562563    1e-04
 top.ef   1         1  0.3057593395 0.05631552    1e-04
 top.ef   2         2  0.9729683133 0.13836033    1e-04
 top.Ram   1         2 -0.2442409485 0.11593499       -1
 top.Rcm   1         2  0.0629173480 0.24414338       -1
 top.Rem   1         2  0.0883929800 0.03553627       -1
 1 x_1 -0.0529413841 0.05435899
 1   y_1  0.1286657693 0.05590730
 1 x_1  0.0457339934 0.05631440
 1   y_1 -0.0864138755 0.04856195

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.

Replied on Mon, 08/12/2019 - 10:51

AdminRobK Joined: 01/24/2014

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).

Replied on Mon, 08/12/2019 - 12:25

tbates Joined: 07/31/2009

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")

Replied on Wed, 08/14/2019 - 11:03

noamm Joined: 10/10/2018

Thanks

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

Replied on Sun, 08/25/2019 - 01:14

AdminNeale Joined: 03/01/2013

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.

News

Recent Posts

Failure to obtain CI's for a umx model

thoughts to help with CIs

DGEBAL

two remarks

mxVersion?

umx bug?

speleology :-)

Some responses

OK, so you don't have an

Thanks again,

some tips

CIs

Thanks

Identification?

News

Recent Posts