[R-sig-Geo] How to check & interpret an MCMCglmm ordinal model in R?
Mirela Beloiu
m|re|@_be|o|u @end|ng |rom y@hoo@com
Mon Feb 3 13:57:06 CET 2020
Dear users,
I'm trying to explain changes in tree vitality from 1 to 3 (1=green, 2=damage, 3=dry) using precipitation, diameter and tree species in an MCMCglmm model. I am struggling with two questions:1. How do I check if the model is correct? 2. How do I interpret the summary of the MCMCglmm model in relation to my response variable (the three groups 1,2 and 3)?
Data and R script can be found here https://drive.google.com/drive/folders/1LmgEAssR5FfFw1CkYjygsaawg84dDAwk?usp=sharing . Details are posted here.
###MCMCglmm model, family= ordinal
prior1<-list(R=list(V=diag(1),nu=0.002))
m1 <-MCMCglmm(VIT_2018~ pp18 +DBH+Species,
family = "ordinal", data = comsp,prior=prior1,pr = TRUE,
nitt = 60000, burnin =30000, thin = 50)
summary(m1)
Iterations = 30001:59951
Thinning interval = 50
Sample size = 600
DIC: -146008.7
R-structure: ~units
post.mean l-95% CI u-95% CI eff.samp
units 2697 1016 5304 6.075
Location effects: VIT_2018 ~ pp18 + DBH + Species
post.mean l-95% CI u-95% CI eff.samp pMCMC
(Intercept) 60.64229 31.71644 92.78733 11.49 < 0.002 **
pp18 -0.07020 -0.11518 -0.03556 15.47 < 0.002 **
DBH 0.17357 -0.01797 0.37281 58.26 0.06000 .
SpeciesBetula pendula -12.88510 -29.42191 1.87807 58.40 0.08667 .
SpeciesCarpinus betulus 15.57570 0.96439 30.89528 45.04 0.02000 *
SpeciesCorylus avellana -3.81337 -18.31965 13.12771 600.00 0.59000
SpeciesCrataegus spec. -14.90077 -36.86880 4.24729 286.00 0.10333
SpeciesFagus sylvatica -15.03559 -29.04547 -2.07809 60.56 0.00667 **
SpeciesFrangula alnus 20.10817 -0.10598 38.66354 73.88 0.01333 *
SpeciesQuercus spec. -9.09458 -24.52595 6.42502 293.33 0.26000
SpeciesSambucus nigra 24.29894 2.58339 46.29901 49.66 0.02333 *
SpeciesSorbus aucuparia 39.56930 22.97282 63.15175 10.13 < 0.002 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Cutpoints:
post.mean l-95% CI u-95% CI eff.samp
cutpoint.traitVIT_2018.1 80.46 55.38 117.9 4.204
#### Estimating Credible Intervals
HPDinterval(mcmc(randomprior1$Sol[,"(Intercept)"]))
# lower upper
#var1 31.71644 92.78733
#attr(,"Probability")
#[1] 0.95
Any suggestions/ideas would be really helpful. Thank you for your time and help!
Mirela
[[alternative HTML version deleted]]
More information about the R-sig-Geo
mailing list