Failure to obtain CI's for a umx model

Posted on
No user picture. noamm Joined: 10/10/2018
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
Picture of user. tbates Joined: 07/31/2009

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 - 10:38
Picture of user. AdminRobK Joined: 01/24/2014

In reply to by tbates

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:33
Picture of user. AdminRobK Joined: 01/24/2014

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
Picture of user. tbates Joined: 07/31/2009

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
No user picture. noamm Joined: 10/10/2018

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
Picture of user. AdminRobK Joined: 01/24/2014

In reply to by noamm

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
No user picture. noamm Joined: 10/10/2018

In reply to by AdminRobK

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.

Replied on Mon, 08/12/2019 - 10:51
Picture of user. AdminRobK Joined: 01/24/2014

In reply to by noamm

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
Picture of user. tbates Joined: 07/31/2009

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
No user picture. noamm Joined: 10/10/2018

In reply to by tbates

Thanks, Iv'e emailed you the Rda file. I really appreciate your help.
Replied on Sun, 08/25/2019 - 01:14
Picture of user. AdminNeale Joined: 03/01/2013

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.