[R-sig-ME] MCMCglmm random slope model
Shae Crain
@h@ecr@|n @end|ng |rom gm@||@com
Mon May 13 18:43:43 CEST 2019
Dear List,
Hi, I'm Shae and I'm a grad student working on the chronobiology and
behavioral profile of a spider (Frontinella communis). I'm looking into
whether there's a possible correlation between antipredator behavior and
time of day. My data is quite zero-inflated and overdispersed...therefore
I'm using MCMCglmm.
Quick synopsis of my project so I'm not just throwing a model at yall:
I assayed wariness levels in indivs over a 24hr period (all at the same
time points), twice. Responses are however long they pause (seconds) and
whether they exhibited huddling behavior. I pooled the measurements into
whether the time point was in light and or dark conditions.
I've been using Dingemanse and Dochtermann (2013) instructional paper,
worked through their random intercept examples and now I'm trying to fit
univariate random slope models to answer my further questions regarding the
plasticity of Frontinella's antipredator strategies.(I wish to do bivariate
RR as well but one step at a time).
There are no Bayesians in my department so I'm sending out a lifeline to be
sure that I am doing this correctly (I'm going against my advisor's wishes
to run a binomial logistic regression). I am very new to R and have spent
the last few weeks teaching myself (so thank you for all of your very
helpful papers!)
xdevh is my centered covariate
xjh are individual level averages
I'm using the prior from Bolker et al (2012)'s owl example
>prior_overdisp <- list(R=list(V=diag(c(1,1)),nu=0.002,fix=2),
G=list(list(V=diag(c(1,1e-6)),nu=0.002,fix=2)))
pmc.rs<-MCMCglmm(p~trait-1+at.level(trait,1):xjh+at.level(trait,1):xdevh,
random = ~idh(trait):id, family="zipoisson",rcov=~idh(trait):units,pr=TRUE,
prior = prior_overdisp,,nitt=1300000,thin=1000,burnin=300000,
data=fullpuff ,verbose=FALSE,start=list(QUASI=FALSE))
# I increased the nitt to avoid autocorrelation of the zi variables...I'm
still at a loss on how to interpret them. The zeros in my responses are
there for a reason; if indivs were very bold and refused to huddle/pause at
all.
My model summary:
Iterations = 300001:1299001
Thinning interval = 1000
Sample size = 1000
DIC: 1555.507
G-structure: ~idh(trait):id
post.mean l-95% CI u-95% CI eff.samp
traitp.id 0.701283 0.361328 1.117891 1000
traitzi_p.id 0.000001 0.000001 0.000001 0
R-structure: ~idh(trait):units
post.mean l-95% CI u-95% CI eff.samp
traitp.units 0.6587 0.4952 0.818 1000
traitzi_p.units 1.0000 1.0000 1.000 0
Location effects: p ~ trait - 1 + at.level(trait, 1):xjh + at.level(trait,
1):xdevh
post.mean l-95% CI u-95% CI eff.samp pMCMC
traitp 3.7466 3.2505 4.2624 1000.00 <0.001 ***
traitzi_p -30.6866 -48.3946 -11.6482 2.15 <0.001 ***
at.level(trait, 1):xjh 0.4811 0.3001 0.6716 1000.00 <0.001 ***
at.level(trait, 1):xdevh 0.2335 0.1286 0.3399 942.41 <0.001 ***
---
I do not feel secure in whether the model fit...and I do not know how to
approach the effect sizes of my zi variables. I have a nagging feeling that
I did the model incorrectly and the traceplots for the zi variables are far
from ideal.
In conclusion, my most pressing questions are:
1. Does my model appear to make sense? Would random=~us(trait):id be more
appropriate? Could it be the prior?
2. Could the effect sizes of my zi variables be a product of the model
fitting poorly? Or are they okay? Should I just run the model longer?
3. Whenever I plot my observed responses over my model, some indivs will
have estimates far larger than should be reasonable (like 7500). Could that
be just an artifact and can it be controlled for? Or is it a diagnostic of
a poor fit?
Thank you all so very much in advance for reading this and I'll appreciate
any feedback at all.
-Respectfully,
Shae
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list