| Attachment | Size |
|---|---|
 model.png | 330.63 KB |
Hi,
I would like to fit the twin model in the enclosed figure, but I am not sure if this is feasible.
Can someone help me understanding if the model is identified?
In case, are there any exemples of R code for this model?
Thanks a lot
M
Printer-friendly version
I don't see any reason why it wouldn't work or wouldn't be identified. I don't know of any example script for such a model, though.
Thank you Rob for your reply.
I tried to write the OpenMx syntax for this model. Can I ask you if this seems ok for you?
I also have 2 problems:
1) I found it hard to specify that my distal outcome is binary, so I wrote the script as it was continuous. How can I specify this?
2) My code doesn't work, and I obtain the following Error:
Error: The following error occurred while evaluating the subexpression 'MZ.gamma %&% MZ.phiA + MZ.psiA' during the evaluation of 'MZ.A' in model 'my_model' : Matrices must have same dimensions in as(e1, "dMatrix") + as(e2, "dMatrix")
Here my script:
Thank you very much!
M
That's very strange, because I don't see from your syntax why that error would be occurring. What is your output from
mxVersion()? Maybe you need to encasegamma %&% phiAin parentheses?I'll assume you want to treat it as a threshold variable. There are several changes you'll need to make to your script. First, run the distal outcome through
mxFactor(). Then, you'll need to make an MxMatrix containing the thresholds for the binary variable. That matrix of thresholds needs to go into the MZ and DZ models, and you'll need to provide its name as the value of argumentthresholdstomxExpectationNormal(). Either the intercept or the threshold for the binary variable must be fixed, to identify the location of the latent trait presumed to underlie the observed dichotomy. The scale of the latent trait also needs to be identified. Numerically speaking, the easiest way to do that is to fix its nonshared environmental variance to 1.0. You could instead use an MxConstraint to fix its total variance to 1.0, which provides results that are easier to interpret, but makes optimization more difficult.One comment I have about your script is that you should model the phenotypic mean for the five timepoints in terms of the intercept and slope. That is, there should be a mean intercept and a mean slope, and the model-expected phenotypic mean at a given timepoint would be the mean intercept plus the product of the mean slope and the timepoint's loading on the latent slope.
Thank you again for your reply.
For the Error, I solved the problem updating openMx.
That's exact, I want to treat the binary distal outcome as a threshold variable. I'll try to follow your indications.
For your comment (modeling the phenotypic mean for the five timepoints) I see what you are suggesting, but unfortunately I could not implement that successfully...can you please help me understand how can I modify my code to model these phenotypic means ?
Relatedly, I used mxCheckIdentification() and I obtained that the model is not locally identified, and that the non-identified parameters are all the means of the manifest variables (ManMean in my script). Is it because I didn't model the phenotypic mean as you suggested?
I also would like to show you what I obtained when I ran the model (using mxTryHard, as not possible otherwise):
Not all eigenvalues of Hessian are greater than 0
What does this message mean? Is there something I can do about it ?
When I tried to estimate the following AE model instead of ACE to simplify (AE is what I expect from prior analysis), I obtained the same warning and bizarre results
OK, good!
OK. You'll probably find it helpful to look up an example script that involves a threshold variable in a twin analysis.
I don't know. Perhaps.
That all sounds like it could be symptomatic of unidentification. "Not all eigenvalues of Hessian are greater than 0" suggests that the optimizer didn't find a local minimum of the objective function, and is fairly common when the model is unidentified.
Sure...and I'll also try to sort out some other potential issues for you as well. First, redefine
gammaas:This way, the loadings on the latent intercept will be orthogonal to the loadings onto the latent slope. The orthogonality will improve the numerical stability of the model and make optimization easier. Naturally, it will also change the interpretation of the mean intercept, which will now be the model-expected mean phenotype score at the third, not first, timepoint. With that in mind, define these new objects:
Note that this syntax fixes the mean of the distal outcome to zero, meaning that its threshold (when you get around to that) will need to be a free parameter. You will quite likely also want to adjust the phenotypic mean for age and sex. Now suppose the age of twin #1, sex of twin #1, age of twin #2, and sex of twin #2 have column names in your dataset "age1", "sex1", "age2", "sex2" (and that the sex variables are dummy codes). Then, you would define these additional objects:
This syntax models the effects of age and sex on the time-varying trait as constant across time, which I think is reasonable. You would then redefine the phenotypic mean vector as:
Note that this syntax adjust the distal outcome for the effects of age and sex, but fixes its regression intercept at zero.
I really appreciate your help, thank you so much!
Following your indication I changed my model, and now it has a binary distal outcome specified as threshold (hope I did that in the right way...)
It also has an adjustment for sex (not for age since the twins have the same age).
This is the model:
However, it doesn't work because of the following error:
Error: Unknown reference 'muIntrcpt' detected in the entity 'ManMean' in model 'my_model'
But it seems to me the the muIntrcpt objet is correctly specified and referenced in the model, isn't it?
An additional question: I think that in this model I am not adjusting my distal outcome (t1_6 and t2_6) for the effect of sex...if I wanted to adjust my outcome for sex, I should replace this:
by this:
right?
Thanks again
M
Objects
iLoadings,sLoadings,meanIntercept, andmeanSlopehave to go into themxModel()statement somehow as well (perhaps as part of listpars). The objects have to be part of the MxModel's namespace for all this to work properly.Once you fix that, though, your model still isn't going to run. The syntax I gave you for the phenotypic mean assumed that
selVarswould be defined as it was in your first attempt at a script. But in your newest script, you dropped the two distal outcomes fromselVars. Since you providedselVarsas thedimnamesargument tomxExpectationNormal(), OpenMx is only going to model the ten data columns containing the time-varying trait. But the model-expected, twin-pair phenotypic mean vector and covariance matrix are too big, i.e. 12x12 and 1x12, respectively. That's going to cause an error.Your specification of the threshold looks OK to me (although you could just directly use the 'threshold' MxMatrix, that is, the 'expThresh' MxAlgebra isn't necessary).
The syntax I gave you adjusts the "liability-scale" mean of the distal outcome for sex. You could adjust its threshold for sex instead, similarly to what you propose, except you'd need to define a new 1x1 MxMatrix to hold the free parameter for the effect of sex on the threshold. If you do that, then fix parameter 'bSexd' to zero in
betaSex.Thank you Rob - and sorry for my very basic mistakes!
Now the code works.
I obtained the following Warning:
Warning message:
In model 'my model' Optimizer returned a non-zero status code 5. The Hessian at the solution does not appear to be convex. See ?mxCheckIdentification for possible diagnosis (Mx status RED).
You said that it may me symptomatic of non-identification. Any suggestion on how to solve the problem?
Good to hear.
Status RED can occur even when the model's identified. You could again try
mxCheckIdentification()and see what it says. If it say the model isn't locally identified, then lingering unidentification is likely the cause.Otherwise, see here for a thread about what to do about status RED in general. Status RED is also quite common in analyses of threshold traits, sometimes unavoidably so, even when the optimizer has indeed found a local minimum. In your particular case, once you're pretty sure the model is identified, then be sure you're using the CSOLNP optimizer, and replace
mxRun()in your script withmxTryHardOrdinal().I tried with mxCheckIdentification() but it gave me this error:
Error in dimnames(x) <- dn :
length of 'dimnames' [2] not equal to array extent
Of course, I don't understand what it is about, a part that 2 objects are not of the same size while they are expected to be of the same size.
Is there a way to fix it or to check identification in another way? To understand if my model is really not identified of it is just an issue with the threshold variable
The model with CSOLNP optimizer (or SLSQP, I tried this too) and mxTryHardOrdinal() gives pretty much the same results as the previous one, ie always the same problem of Hessian. For instance, at each iteration of mxTryHard() the result is:
Not all eigenvalues of Hessian are greater than 0: ....
That might be a bug.
I think it's pretty likely to be due to the threshold variable.
If you're getting pretty much the same results--as in, similar parameter estimates and similar -2logL values--then that's a pretty good sign that the model is indeed identified, and that the optimizer is finding a local minimum. I think you are probably OK.
Thank you very much for your help.
Indeed, I obtained the same -2LogL and parameter estimates with all the optimizers.
However, I am not sure that everything is right when I look at the standardized variance components...
Here the results (AE model):
The are 2 values >1, which should not be possible....
Any idea of why it is the case? Or potential solutions?
You are mistaken. Those three "standardized" matrices are respectively the common-factor A, C, and E covariance matrices, elementwise-divided by the total common-factor covariance matrix. Their off-diagonals can certainly be greater than one, or less than zero. More or less this same topic has been discussed before on these forums here, here, here, here, and here.
Perhaps you want to use, e.g.,
cov2cor(phiA)instead?ok, I see. Thank you for your explanations, I really appreciated