Hi,
The residual variance is not identifiable in the likelihood for
ordinal models - you need to fix it at something (e.g. 1):
prior=list(R=list(V=1, fix=1), ...)
Note that under family="ordinal" this is equivalent to probit
regression but with a variance of 2 (rather than the standard 1) in
the normal cumulative distribution function (i.e. pnorm).
family="threhold" fits the same model but fixing the unit variance at
1 gives the standard probit regression. This implementation samples
the posterior more efficiently too.
Cheers,
Jarrod
Quoting Luka Rubinjoni <rubinjoni at gmail.com> on Wed, 11 Feb 2015
15:13:59 +0100:
Dear list members,
I'm trying to analyze ordinal data with MCMCglmm, but I keep getting
inconsistent results. I'm a biologist with no experience in Bayesian
inference and mixed models, these are my first attempts...
I'm exploring herbivory load on a dioecious plant (separate male and
female sex). I used leaf damage as an estimate of herbivory load, by
assigning each leaf with a 1 (no damage) to 4 (heavy/severe damage)
score. I have sampled about 500 plants from 12 locations, and have about
10 leaves per plant. Plants do not have a constant number of leaves, it
ranges from 6 to 17. All plants were collected from natural habitats -
there was no experimental setup for a full factorial design. Also,
uneven number of plants was sampled from various locations.
I wish to determine if one of the sexes (male or female) is more
susceptible to herbivore damage. I would also like to compare sites
according to the observed damage, and explore sex*site interaction.
Further, I would like to see if some other plant traits (except sex)
have an impact on herbivory: plant height, plant mass, leaf number, leaf
mass, flower/fruit mass & number, reproductive allocation...
I tried running MCMCglmm procedure with the default prior (without
declaring the prior), but I get awkward results (negative DICs, large
differences in DICs between runs...). I'm not sure if I'm making any
other mistakes, too. Do you have any suggestions for the prior
specification?
Best regards,
Luka Rubinjoni
Faculty of Biology, University of Belgrade
P.S. Here's summary of the data, and output from one of the models:
ID LeafID score N
K24M : 17 Min. : 1 1:2248 Min. : 6.00
M9M : 17 1st Qu.:1480 2:2982 1st Qu.:10.00
J10M : 16 Median :2960 3: 564 Median :11.00
K18M : 16 Mean :2960 4: 124 Mean :11.35
M17F : 16 3rd Qu.:4439 3rd Qu.:12.00
M19M : 16 Max. :5918 Max. :17.00
(Other):5820
PlantH FlowerN InternodeN
Min. : 94.0 Min. : 1.0 Min. : 4.000
1st Qu.:204.0 1st Qu.: 5.0 1st Qu.: 5.000
Median :231.0 Median :11.0 Median : 6.000
Mean :233.7 Mean :15.6 Mean : 6.131
3rd Qu.:265.0 3rd Qu.:24.0 3rd Qu.: 7.000
Max. :398.0 Max. :68.0 Max. :10.000
PlantM LeafM Sex Site
Min. :0.1690 Min. :0.1170 F:2891 A2S : 875
1st Qu.:0.4760 1st Qu.:0.2950 M:3027 A : 805
Median :0.6120 Median :0.3910 K : 784
Mean :0.6816 Mean :0.4316 M : 566
3rd Qu.:0.8500 3rd Qu.:0.5310 DJ : 479
Max. :2.3020 Max. :1.5610 P : 454
(Other):1955
RA InfloresenceN FlowerM
Min. :0.001348 Min. :1.000 Min. :0.00100
1st Qu.:0.011674 1st Qu.:4.000 1st Qu.:0.00700
Median :0.024943 Median :4.000 Median :0.01700
Mean :0.065849 Mean :4.492 Mean :0.04613
3rd Qu.:0.097990 3rd Qu.:5.000 3rd Qu.:0.05600
Max. :0.838617 Max. :9.000 Max. :0.56100
NA's :2891
CMC200kIDSite <- MCMCglmm(score~Sex, ~ID+Site, rcov=~units,
family="ordinal", data=skorovi, DIC=TRUE, nitt=200000, thin=100,
burnin=5000)
...
MCMC iteration = 200000
Acceptance ratio for latent scores = 0.897869
Acceptance ratio for cutpoint set 1 = 0.779000
Iterations = 5001:199901
Thinning interval = 100
Sample size = 1950
DIC: -9100.507
G-structure: ~ID
post.mean l-95% CI u-95% CI eff.samp
ID 12.02 0.3165 24.33 1.436
~Site
post.mean l-95% CI u-95% CI eff.samp
Site 17.89 0.1927 52.11 5.023
R-structure: ~units
post.mean l-95% CI u-95% CI eff.samp
units 40.44 0.198 80.42 1.478
Location effects: score ~ Sex
post.mean l-95% CI u-95% CI eff.samp pMCMC
(Intercept) 2.38948 0.01827 5.96931 5.605 0.0359 *
SexM 0.11478 -0.64723 0.90102 1950.000 0.7426
---
Signif. codes: 0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1
Cutpoints:
post.mean l-95% CI u-95% CI eff.samp
cutpoint.traitscore.1 10.74 2.051 16.73 1.347
cutpoint.traitscore.2 16.79 3.203 26.18 1.323
autocorr(MCMC200kIDSite$Sol)
, , (Intercept)
(Intercept) SexM
Lag 0 1.0000000 0.02834661
Lag 100 0.4652665 0.12773192
Lag 500 0.4637719 0.12885704
Lag 1000 0.4608056 0.10518225
Lag 5000 0.4426579 0.10261122
, , SexM
(Intercept) SexM
Lag 0 0.02834661 1.000000000
Lag 100 0.10172481 0.000939134
Lag 500 0.13601706 0.003533253
Lag 1000 0.13280788 0.040197635
Lag 5000 0.11814053 0.070899772
autocorr(MCMC200kIDSite$VCV)
, , ID
ID Site units
Lag 0 1.0000000 0.7125238 0.9873699
Lag 100 0.9770221 0.7118302 0.9862788
Lag 500 0.9692561 0.7054513 0.9799585
Lag 1000 0.9644127 0.6951924 0.9737191
Lag 5000 0.9104725 0.6701422 0.9193135
, , Site
ID Site units
Lag 0 0.7125238 1.0000000 0.7214728
Lag 100 0.7095693 0.4868717 0.7170873
Lag 500 0.7028373 0.5282474 0.7136912
Lag 1000 0.7045606 0.5254454 0.7081390
Lag 5000 0.6686113 0.4824315 0.6739685
, , units
ID Site units
Lag 0 0.9873699 0.7214728 1.0000000
Lag 100 0.9864973 0.7192773 0.9971748
Lag 500 0.9810497 0.7137308 0.9910396
Lag 1000 0.9743296 0.7063965 0.9838989
Lag 5000 0.9213108 0.6722353 0.9303135