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.
The BLAS error code looks like a bug for OpenMx or its CSOLNP component? I've not seen it before
But for profile confidence intervals, which can be tricky to compute, switching optimizer when a CI isn't returned can often help.
For what you are calling "any significance in general", I would be testing if parameters can be dropped or not, with
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
That's strange. We don't even call DGEBAL directly. We'd need a stack trace to diagnose further.
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.
What do you get from
mxVersion()
?Incidentally, this,
, looks like a umx bug. In a multigroup model, the "container" MxModel generally doesn't have an MxData object.
FYI
umxCI_boot
function is from here https://openmx.ssri.psu.edu/thread/2598Now 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.0umx
release this Fall.Hi all, and thanks for all the responses.
Some answers:
Regarding
Iv'e tried to use both optimizers but the first one have failed:
However, when using the SLSQP optimizer, the CI's computation have failed:
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.
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.
I don't know the inner workings of
umxConfint()
. Could you run your model withintervals=TRUE
, and then runsummary()
on it withverbose=TRUE
? That will show diagnostics for why none of the confidence limits was acceptable.Thanks for posting your
mxVersion()
output. Thesummary()
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.Thanks again,
I've tried again with CSOLNP and running the model with
intervals=TRUE
andverbose=TRUE
: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.
.
I would really appreciate any other suggestions.
Thank you,
Noam.
summary()
withverbose=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 withmxRun()
instead? Also, it looks like the DGEBAL error didn't happen with CSOLNP and version 2.13, correct?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).
e-mail me the model saved as an Rda and I will have a look at CIs for you.
save(m1, file="m1.rda")
Thanks, Iv'e emailed you the Rda file. I really appreciate your help.
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.