| Attachment | Size |
|---|---|
 joint_ord_cont_model.R | 13.22 KB |
Hi everyone,
I am trying to run a joint ordinal-continuous model.
In my model I have a continuous variable (varA) and an ordinal variable with 1 threshold (varB).
There are 5 zygosity groups (MZM, DZM, MZF, DZF, DOS) and a covariate age which
I modeled separately for men and women.
The model is actually running perfectly fine (no errors) and the estimates for the thresholds for varB
look good and are what I would have expected.
However, the estimates for the means are all way too low when I compare them with the raw data.
Yet the pattern of the means between the groups is similar to the raw data (higher means in males than females).
Does anyone have an idea as to why my means are lower than they should be?
Attached is my script.
Thank you very much in advance!
Jorien
Printer-friendly version
I don't see anything obviously wrong with your script. It would help if you provided the direct and maximum-likelihood estimates of the means, to show the discrepancy between them.
Also, how much missing data are you dealing with, and how big is the sample size?
thanks for the quick reply!
These are the direct means:
MZM DZM MZF DZF DOS-M DOS-F
220.4 225.7 136.7 120.9 222.6 114.6
And the maximum-likelihood estimates of the means:
MZM DZM MZF DZF DOS-M DOS-F
184.0 197.2 91.4 84.4 170.5 80.6
The sample size is 10,368 twins, of which 6,866 from complete pairs and 3,502 from incomplete pairs
(1,425 MZM, 907 DZM, 3,541 MZF, 1,948 DZF and 2,547 DOS).
The strange thing is that the thresholds for the ordinal variable seem fine.
Would be great if you can figure it out.
How are you computing the means from your MxModel? I suspect that you're pulling the values of the intercepts for each kinship group, from the MxMatrices named 'MMZM', 'MDZM', etc. Since you are including the regression onto age in your model for the phenotypic mean, the model-expected mean ("yhat," so to speak) likely differs from person to person within a group, and therefore, it doesn't really make sense to talk about just one "mean" for a group. If the regression slope w/r/t age is positive, then I would anticipate that the intercepts would be smaller than the "direct" grand means (from
colMeans()?) of the dataset.FYI, to find the model-expected mean for a given row of your dataset, you could do something like
You would change the value of
defvar.rowas needed.ohhh wait I think you are right.. when I run the model without age,
then the maximum likelihood estimates are:
MZM DZM MZF DZF DOS-M DOS-F
223.2 234.3 141.3 124.5 221.8 113.4
Which is very similar to the direct means.
Makes sense now that I think about it (should have thought about it before :-) )
Thank you so much!