# Bivariate twin model (output)

7 posts / 0 new
Offline
Joined: 03/07/2017 - 09:59
Bivariate twin model (output)

Hello,

I've done some univariate twin model analyses in OpenMx recently, but now I'm trying my hand at bivariate models.
I'm using a script by Hermine (http://ibg.colorado.edu/cdrom2016/maes/MultivariateAnalysis/mulACEc2.R).

So far I succeeded in running the model and obtaining indicators of fit and path loadings. But I've got a couple of questions:
1) However, the output I get does not produce viable 95%CI's (the model summary shows lbounds converging on zero and no ubounds). How can I obtain proper intervals around my path loadings?

2) And to answer my research question, about the influence of A, C, and E on the correlation between variables, I need the cross-trait within twin correlation, and the cross-trait cross-twin correlation. So far, I haven't been able to find these (but maybe they are in this script and I just read over them).

3) So far I've compared the satured model with models where A, C, and E are included for both variables, but I would also like to try models that include e.g. ACE for variable 1 and CE for variable 2, or CE for 1 and E for 2. Where can I find a script for such models?

Cheers,
Jizzo

Offline
Joined: 01/24/2014 - 12:15
answers
1) However, the output I get does not produce viable 95%CI's (the model summary shows lbounds converging on zero and no ubounds). How can I obtain proper intervals around my path loadings?

Could you provide the actual output? Try running summary() on your fitted MxModel, with argument verbose=TRUE. That will print the details of the confidence-limit solutions. There might not be anything wrong with the lower limits, but upper limits of NA could indicate a problem. I've learned from experience that novice users are sometimes confused by CI results that are in fact correct (e.g., threads 4253, 3999, 3446).

In any event, there are three reasons why a confidence limit would be displayed as NA. First, it could be equal to the point estimate, in which case it may actually be correct if the point estimate is is on a boundary of the feasible space (or is actually fixed, like a diagonal element of a correlation matrix). Second, it could be on the wrong side of the point estimate, such as an upper limit that is less than the point estimate, which is obviously wrong. Third, the corresponding change in -2logL from the MLE could differ too much from its target value (which for a 95% CI is about 3.841); in this case, the limit is wrong in that the resulting confidence interval is too narrow or too wide.

2) And to answer my research question, about the influence of A, C, and E on the correlation between variables, I need the cross-trait within twin correlation, and the cross-trait cross-twin correlation. So far, I haven't been able to find these (but maybe they are in this script and I just read over them).

You'll want to convert the model-expected MZ and DZ covariance matrices into correlation matrices. See here and here for suggestions.

3) So far I've compared the satured model with models where A, C, and E are included for both variables, but I would also like to try models that include e.g. ACE for variable 1 and CE for variable 2, or CE for 1 and E for 2. Where can I find a script for such models?

You can modify the script you're using to selectively fix the appropriate path coefficients to zero, with omxSetParameters(). For instance, to fit a model that has ACE for the first trait and CE for the second, fix to zero the [1,2] and [2,2] elements of the "A" path-coefficient matrix. In your case, the command is probably something like omxSetParameters(fitACE, labels=c("a12","a22"), free=FALSE, values=0).

Offline
Joined: 03/07/2017 - 09:59
First of all, thanks for your

First of all, thanks for your comments.

1) The output for the ACE model is as follows:

Summary of mulACEc

free parameters:
name matrix row col     Estimate   Std.Error A lbound ubound
1   mean_SSRTK  meanG   1   1 -4.218330819 0.072103637
2  mean_MARSHK  meanG   1   2  1.999382994 0.042575069
3        a_1_1      a   1   1  0.404005238 0.328936236    1e-04
4        a_2_1      a   2   1  0.224590135 0.126932020      -10
5        a_2_2      a   2   2  0.000100000 0.336833820 ! 1e-04!
6        c_1_1      c   1   1  0.554463862 0.188160774    1e-04
7        c_2_1      c   2   1  0.056627754 0.147461734      -10
8        c_2_2      c   2   2  0.000100000 0.990393931 ! 1e-04!
9        e_1_1      e   1   1  0.912123905 0.064540022    1e-04
10       e_2_1      e   2   1 -0.032189192 0.060461897      -10
11       e_2_2      e   2   2  0.713286712 0.041982278    1e-04

Model Statistics:
|  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
Model:             11                    657             1765.1873
Saturated:             NA                     NA                    NA
Independence:             NA                     NA                    NA
Number of observations/statistics: 167/668

Information Criteria:
|  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
AIC:      451.18729              1787.1873                       NA
BIC:    -1597.33464              1821.4852                1786.6576
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: 2017-06-12 14:34:56
Wall clock time (HH:MM:SS.hh): 00:00:00.72
optimizer:  SLSQP
OpenMx version number: 2.7.11
Need help?  See help(mxSummary) 

2) I got the correlation matrix alright.
3) Thanks, no further questions.

4) New issue:
The best fitting model we used (so far) was a CE model, with the following path loadings:
C1 --> Var1 0.5674
C1 --> Var2 0.2017
C2 --> Var2 0.1838
E1 --> Var1 0.8234
E1 --> Var2 -0.0098
E2 --> Var2 0.9620
This negative path loading troubles me..
More so, because I need this to calculate to what extent C and E explain the variance in the correlation.
With a cross-trait, within-twin correlation of 0.1083, this leads to impossible proportions of explained variance.
C [ 0.5674 * 0.2017 ] / 0.1083 = 1.06
E [ 0.8234 * -0.0098] / 0.1083 = -0.07

Offline
Joined: 01/24/2014 - 12:15
1) Concerning the output-

1) Concerning the output--confidence-interval results appear in a separate table from the "free parameters" table. The "lbound" and "ubound" columns of the "free parameters" table contain the lower and upper boundaries on the allowed range of the free parameters. Confidence intervals aren't automatically calculated. Are you including an MxInterval object (such objects are created with mxCI()) in your MxModel, and passing argument intervals=TRUE to mxRun()?

4)
Basically, the interpretation of those quantities as proportions comes apart if any of the components that sum to the phenotypic correlation are negative. You are not the first person to be confused when this happens, which I think reflects a shortcoming in how twin modeling is typically taught.

Offline
Joined: 03/07/2017 - 09:59
Intervals

1) Thanks for the comment. With input from other forum posts, I have updated my model to include

 ciACE     <- mxCI(c("VC[1,1]", "a_1_1", "a_2_1", "a_2_2", "c_1_1", "c_2_1", "c_2_2", "e_1_1", "e_2_1", "e_2_2"))

This does provide me with intervals in the output as intended. However, the output I am mostly interested in are the ACE  Standardized Path Coefficients (pre-multiplied by inverse of standard deviations) and the squared standardized path coefficients. Now the next step would be to obtain the 95%CI's around the standardized and squared standardized coefficients, but I can't quite figure this out.
Offline
Joined: 03/07/2017 - 09:59
Question

Bit of a formatting error there, my question is:

This does provide me with intervals in the output as intended. However, the output I am mostly interested in are the ACE Standardized Path Coefficients (pre-multiplied by inverse of standard deviations) and the squared standardized path coefficients. Now the next step would be to obtain the 95%CI's around the standardized and squared standardized coefficients, but I can't quite figure this out.

Offline
Joined: 01/24/2014 - 12:15
Include MxAlgebra names in mxCI()'s 'reference'

It's hard to make specific suggestions without seeing the script you're currently working from, but the simple answer is to include in mxCI() argument reference the names of the MxAlgebra objects that contain the standardized and squared standardized path coefficients, e.g.

ciACE     <- mxCI(c("VC[1,1]", "a_1_1", "a_2_1", "a_2_2", "c_1_1", "c_2_1", "c_2_2", "e_1_1", "e_2_1", "e_2_2",
"stdPathA","stdPathC","stdPathE","stdPathA2","stdPathC2","stdPathE2"))

or whatever the case may be. Note that mxCI() requires the MxAlgebras' names (i.e., the character strings provided to argument name of mxAlgebra()), and not their R symbols (i.e., the thing to the left of assignment operator  <-  in your syntax).

I don't know whether you've put 'ciACE' into the "container" MxModel or into one of its submodels ("groups"), so if the algebras for the standardized coefficients are in a submodel but 'ciACE' is in the container model, you might have to prefix the MxAlgebra names with the name of the appropriate submodel--say, "MZ." (for instance, "MZ.stdPathA"), or whatever the case may be.

If you do an advanced Content search of the OpenMx website, for forum topics containing the phrase "confidence intervals", you might find additional hints or illustrative examples.

Log in or register to post comments