tbates / umx

Making Structural Equation Modeling (SEM) in R quick & powerful
https://tbates.github.io/
44 stars 17 forks source link

How to add CIs in umx #109

Closed qingwending closed 4 years ago

qingwending commented 4 years ago

I made a common pathway model in umx. I added CIs through adding umxCI(model, run="yes"). However, the output are NAs and there is a warning that "3: The nonlinear constraints and bounds could not be satisfied. The problem may have no feasible solution." I wonder what it means and what I can do to deal with it.

Here is my information about version and optimizer. umx version: 3.0.1 OpenMx version: 2.14.11 [GIT v2.14.11] R version: R version 3.5.3 (2019-03-11) Platform: x86_64-apple-darwin15.6.0 MacOS: 10.14.6 Default optimizer: SLSQP

Here is my script. setwd("....") getwd()

Load Data; return a data frame data <- read.spss(".....", to.data.frame = TRUE)

Pick the variables selDVs = c("x","y","z") mzData <- subset(data, zygosity=="1") dzData <- subset(data, zygosity=="2") umx_set_optimizer(opt = "SLSQP")

Common Pathway Models cp <- umxCP(selDVs = selDVs, dzData = dzData, mzData = mzData, sep = "",tryHard = "yes") cp <- umxCI(cp, run = "yes")

Here is the output lbound estimate ubound lbound Code ubound Code a_cp_r1c1 NA 0.806 NA NA NA c_cp_r1c1 NA 0.000 NA NA NA e_cp_r1c1 NA 0.592 NA NA NA top.as_std[1,1] NA 0.000 NA NA NA top.as_std[2,1] NA 0.000 NA NA NA top.as_std[3,1] NA 0.000 NA NA NA top.as_std[1,2] NA 0.000 NA NA NA top.as_std[2,2] NA 0.000 NA NA NA top.as_std[3,2] NA 0.000 NA NA NA top.as_std[1,3] NA 0.000 NA NA NA top.as_std[2,3] NA 0.000 NA NA NA top.as_std[3,3] NA 0.067 NA NA NA top.cs_std[1,1] NA 0.452 NA NA NA top.cs_std[2,1] NA 0.000 NA NA NA top.cs_std[3,1] NA 0.000 NA NA NA top.cs_std[1,2] 0 0.000 NA 3 NA top.cs_std[2,2] NA 0.172 NA NA NA top.cs_std[3,2] NA 0.000 NA NA NA top.cs_std[1,3] NA 0.000 NA NA NA top.cs_std[2,3] NA 0.000 NA NA NA top.cs_std[3,3] NA 0.315 NA NA NA top.es_std[1,1] NA 0.708 NA NA NA top.es_std[2,1] NA 0.000 NA NA NA top.es_std[3,1] NA 0.000 NA NA NA top.es_std[1,2] NA 0.000 NA NA NA top.es_std[2,2] NA 0.472 NA NA NA top.es_std[3,2] NA 0.000 NA NA NA top.es_std[1,3] NA 0.000 NA NA NA top.es_std[2,3] NA 0.000 NA NA NA top.es_std[3,3] NA 0.510 NA NA NA top.cp_loadings_std[1,1] NA 0.543 NA NA NA top.cp_loadings_std[2,1] NA 0.864 NA NA NA top.cp_loadings_std[3,1] NA 0.797 NA NA NA expMean_Rumination1 NA 24.225 NA NA NA expMean_Anxiety1 NA 38.829 NA NA NA expMean_Depression1 NA 41.378 NA NA NA as_r1c1 NA 0.000 NA NA NA as_r2c2 NA 0.000 NA NA NA as_r3c3 NA 0.590 NA NA NA cs_r1c1 NA 2.677 NA NA NA cs_r2c2 NA 1.516 NA NA NA cs_r3c3 NA 2.766 NA NA NA es_r1c1 NA 4.194 NA NA NA es_r2c2 NA 4.152 NA NA NA es_r3c3 NA 4.486 NA NA NA cp_loadings_r1c1 NA 3.213 NA NA NA cp_loadings_r2c1 NA 7.598 NA NA NA cp_loadings_r3c1 NA 7.006 NA NA NA [1] "NA: " [1] "3: The nonlinear constraints and bounds could not be satisfied. The problem may have no feasible solution."

mcneale commented 4 years ago

It looks to me like a data or identification problem. You seem to have 3 observed variables per person. However, you have specified 1 common pathway factor, and residuals that also generate covariance (as_r1c2 for example). I would try without the r1c2 r1c3 r2c3 paths, and check with @tbates about this default (as long as it doesn’t make him ruminate, anxious or depressed :)). I expect that mxCheckIdentification(cp) would complain.

It is often helpful to run the psych package describe() function on the data after reading them in - to let you inspect N, mean, variance, max, min. This can help developers provide support.

tbates commented 4 years ago

Hi, Yes: a 1-common factor model with uncorrelated residuals and three manifests comes back as locally unidentified (working example model at bottom of post):

mxCheckIdentification(m1)
$non_identified_parameters
[1] "cs_r2c2" "cs_r3c3"

Perhaps a data issue in this case?

Regarding apparent correlated as: The default for correlated as (freeLowerA) is FALSE (In a print out from OpenMx, off-diagonal values are included because as_std is an algebra, so OpenMx prints the algebra as a square result).

m1$top$as
LowerMatrix 'as' 
$free
      [,1]  [,2]  [,3]
[1,]  TRUE FALSE FALSE
[2,] FALSE  TRUE FALSE
[3,] FALSE FALSE  TRUE

PS: for CIs, I'd use m1 = umxConfint(m1, parm= "smart") the smart option knows enough about CP models to only spend time estimating meaningful parameters in the model.

Here's a working example with 1 common factor and 3 manifests to play with. It returns fine, but struggles to compute CIs.


require(umx)
  data(GFF)
  mzData = subset(GFF, zyg_2grp == "MZ")
  dzData = subset(GFF, zyg_2grp == "DZ")

  m1 = umxCP("new", selDVs = c("qol", "hap", "sat") , sep = "_T", nFac = 1, dzData = dzData, mzData = mzData, tryHard = "yes", optimizer = "SLSQP")

 Solution found!  Final fit=68774.132 (started at 314346.08)  (11 attempt(s): 11 valid, 0 errors)

new -2 × log(Likelihood) = 68774.132
## Common Factor paths

|                |     A|     C|     E|
|:---------------|-----:|-----:|-----:|
|Common.factor.1 | 0.612| 0.266| 0.745|
## Loading of each trait on the Common Factors

|    |   CP1|
|:---|-----:|
|qol | 0.689|
|hap | 3.667|
|sat | 4.499|
## Specific-factor loadings

|           |   qol|hap   |sat   |
|:----------|-----:|:-----|:-----|
|Specific a | 0.408|0.945 |1.149 |
|Specific c | 0.207|.     |.     |
|Specific e | 0.744|2.367 |2.131 |
mcneale commented 4 years ago

Oh I see re: CIs for fixed elements. It seems to me a bit dumb to report NAs for fixed mxMatrix elements - these should be the values themselves (since they can never take any other value, all the confidence intervals are the same as the value).

Anyway, IMO, 3 variable 1 factor CP model should be identified, as long as there are sufficient data. When any zygosity-specific N approaches 2*NumVariables things will go awry in identification land. Some high correlations between variables may also cause problems. This is especially likely when variables are subsets of each other - everyone who reports depression also reports rumination, for example. This data pattern can be a problem even if only one twin in one zygosity has it.

tbates commented 4 years ago

reporting CIs as the found value instead of NA is logged as an issue, but going a bit stale. https://github.com/OpenMx/OpenMx/issues/78

tbates commented 4 years ago

Great info on identification @mcneale ! I'll close this issue here and leave a link back to the OpenMx umx sub forum where this thread started:

https://openmx.ssri.psu.edu/node/4550