Attachment | Size |
---|---|

Script binary depvar | 10.49 KB |

Script ordinal depvar | 11.15 KB |

Hello,

I wrote a script to test a threeway GxE where one moderator is continuous (mod1) and the other is binary (mod2) and the dependent variable (depvar) is ordinal (4 categories). While there is no within-twin-variance of mod1 (family-level-variable education of parents), the values of mod2 can differ within a pair (MZ and DZ) (birth weight of twins). So there are six groups in total: MZ high and high; MZ low and low; MZ high and low and the same for DZ's. The aim of the model is to test whether the moderation of the effect of one of the components on depvar by mod1 are the same in the groups defined by mod2 (high weight and low weight). There are no NAs on mod1 or mod2 in the sample used for the analysis. The "oppsite pairs" are always in the same order so that twin 1 is high and twin 2 is low. The model runs fine with mxTryHardOrdinal() and the solution seems robust across optimizers (NPSOL and SLSQP; CSOLNP seems to be very slow: it does not report any progress like the best fit so far, so I stopped it after 20 minutes). But I get the following warning:

The model does not satisfy the first-order optimality conditions to the required accuracy, and no improved point for the merit function could be found during the final linesearch (Mx status RED)

In several threads in this forum I have read that this warning does not necessarily mean that there is something wrong with the model when using the liability threshold approach. But on the other hand I am not sure whether I made a mistake in the script (the multi-group-approach and moderation via the definition variable in the same script gave me some headache...) which causes the warning. So I would appreciate if somebody could check the script and tell me if there is something wrong with it.

For the SLSQP-runs I am using the OpenMx version 2.17 with R 4.0.2 on a Windows (64 bit). For the NPSOL-runs I downgraded R to 3.6.3 on the same machine since the NPSOL-OpenMx does not work with R>=4 at the moment. Also, following the message

Your ordinal model may converge if you reduce mvnRelEps

try: model <- mxOption(model, 'mvnRelEps', mxOption(model, 'mvnRelEps')/5)

I reduced the mvnRelEps, but that did not make the warning disappear. I collapsed the ordinal dependent variable into a binary variable (constraining Var(Y|M=0)=1 and Mean(Y)=0) but here the same issue arises.

I thought about running the models for the high and low group seperately but I am not sure how to deal with the opposite pairs.

Thank you!

The likelihood of ordinal data under the threshold model is evaluated by numerical integration of the multivariate normal distribution over the region of the response. The limited numerical precision of the integration is likely the culprit for the warnings here. I suggest you try a few other starting values and see if the parameter estimates and likelihood are the same. If they are it increases your confidence that you are at the right solution. One way to automate this process is to use mxTryHardOrdinal() instead of mxRun().

Interesting model!

Hi, ok I understand. Thanks for the clarification. I am already using MxTryHardOrdinal() with 20 exhaustive extra tries, so reading your comment it seems that the solution can be considered as an acceptable one.

Thank you!

Hi, I have one more question. When I use a higher number of tries, then some of the models switch the sign of some of the path coefficients, while the model has more or less the same fit. I have read about the sign indeterminacy here in the forum and sometimes it was proposed to fix one of the parameter to be positive. Since the sign of the moderating path coefficients are crucial for my research question I would like to be sure to chose the model with the "right" sign. Could you tell me how to deal with this phenomenon in the models I built?

When a coefficient changes sign, does its new estimate have approximately the same absolute value as its previous estimate? Because if not, then it doesn't seem like sign indeterminacy. In fact, I don't think sign indeterminacy is even possible in a Purcell-style moderation model. I have to ask why the sign of the moderating path coefficients is of substantive interest. The corresponding variance component is quadratic in the coefficients and moderators. To see how the variance components change as function of the moderators, you'll probably have to plot them to get a sense of what's going on.

Also, since birth weight has both between- and within-family variance, you will probably want to fit the so-called "full bivariate" form of the Purcell model. See this paper for reference.

I was using the phone and I meant "One more question" in the title and the fit is not more or less the same. It is the same.

I re-checked the results and it seems as if I commented out the set.seed() so

and the results mxTryHardOrdinal() with different numbers of extra tries differed in the absolute values of the sign-switching coefficients, the SEs and the model fit differed like in the fifth and/or sixth decimal place. Then I realized that I forgot to set a seed with set.seed(). Having done that, the results with different numbers of extra tries were consistent.

Concerning the between and within-pair variance of birth weight: To model the threeway interaction I am using the dichotomized birth weight variable as a grouping variable with 6 groups (For MZ and DZ: High and Low birthweight; High and High; Low and Low). I want to find out whether the moderating path coefficients differ significantly between twins with high and low birth weight by constraining the parameters to be equal. This approach seemed to me to be more feasible than modeling a three-way interaction using the bivariate Purcell model.

I am interested in the effects because of two reasons: 1. the reviewers are interested in the effects (sociologists are more accustomed to look at effects instead of variance components); 2. my plan was to check the significance of the three-way moderation by testing wether there are significant differences in the moderating paths between the high and low group.

But I see your point about checking the (raw) variances and not the path coefficients and your comment made me think more in detail about how to deal with the three-way interaction.

In the two-way case - following your advise and the suggestions of Purcell - I would check the significance of the moderation by testing the moderating path and check the sign/shape of the moderation by plotting the (raw) variances. In my three-way-interaction model I want to answer the question, if the moderation of the ACE components differs significantly between the high and low weight group. But there are cases where the moderation of e.g. A may differ significantly between the groups but the shape of the raw variance in a sensible value range of mod1 for my case (e.g. -1 to +1) is the same. E.g.: Group 1: 1-0.5mod1; Group 2: -1+0.5mod1. The coefficients (-0.5 and +0.5) may differ significantly, but plotting the raw variance of A as a function of mod1 in both groups shows that the shape is the same (it decreases in both cases). Another case: Group 1: 1-0.5mod1; Group 2: -1-0.5mod1. The moderation coefficients (+0.5 and +0.5) may not differ significantly, but plotting the raw variance of A as a function of mod1 in both groups shows that the shape is the not same same in the relevant value range (-1 to +1: decreasing in group 1 and increasing in group 2).

So answering the questions "Moderation significant?" and "Which shape?" with two different statistics (testing path coefficient and plotting raw variance) works fine in the two-way-interaction case but may lead to contradicting results in a three-way-interaction case. Do you have any suggestions how to deal with that situation? Is a multiple group Purcell model the right choice to answer my research question?

To be clear: I questioned the interpretability of the moderation coefficients' signs, and remarked that the variance components are quadratic in the coefficients and moderators, to point out that a given variance component need not necessarily increase or decrease as a function of a moderator across the moderator's observed range. As a function of a moderator, the variance component will be a curve, not a straight line. So, the question for statistical inference is whether or not the moderation effects are distinguishable from zero. Plots or some other description of how the variance components change with the moderators would then provide interpretability of your results.

Dichotomizing birth weight is a mistake. Birth weight is a continuous variable, with a true zero point, and measured in physically meaningful units. Why would you throw that away by dichotomizing it? On top of that, dichotomizing it will make it impossible to detect anything like the curvilinear trends (which I described above) in how each variance component changes as a function of it.

What I'd suggest instead is a "full bivariate" analysis of birth weight and the ordinal trait. In this model, each moderated path coefficient will have the form of

a + b*x + c*y + d*x*y

where 'x' is birth weight, 'y' is parental education, 'a' is the "unmoderated" coefficient, 'b' is the coefficient for birth weight, 'y' is the coefficient for parental education, and 'd' is the coefficient for the product of birth weight and parental education. Be sure to include parental education as a covariate for both birth weight and the ordinal DV.