[R-sig-ME] Estimates for groups within fixed effects & also: enough info for estimating variance? MCMCglmm
Jarrod Hadfield
j.hadfield at ed.ac.uk
Tue Mar 6 10:49:19 CET 2018
Hi,
1/ It is not a good idea to remove the intercept in this instance. It is
forcing the regressions to go through the origin so that when Zdbh and
Znnd are zero the probability of being in either class is 50:50 (0 on
the logit scale). The random effects are not deviations from the first
species but deviations from the average species. Only with fixed effects
are the contrasts (usually) with the first level. If you want random
slopes you need to replace ~species with us(Znnd):species but with only
23 species and a categorical response you will probably not get very
precise estimates.
2/ i.i.d species effects are often hard to separate from
phylogenetically correlated species effects, and the phylogenetic
heritability (aka Pagel's lambda) is often poorly estimated. This is
particularly the case with categorical response variables and you have
the added difficulty that numerical issues are often encountered if
there is support for heritabilities close to one. Presumably there is
variation in the outcome within species? If you do encounter numerical
problems then you can use trunc=TRUE in the call to MCMCglmm which keeps
the latent variable from going into extreme-probability territory or you
can assume the heritability=1 by using MCMCgmmRAM:
http://jarrod.bio.ed.ac.uk/MCMCglmmRAM_2.24.tar.gz
If you have a lot of replication within species (?) my feeling is that
you can probably get away with fitting both terms.
Cheers,
Jarrod
On 05/03/2018 17:08, Drager, Andrea Pilar wrote:
>
> Hi all,
>
> I've fitted a model with a binary response at the individual level,
> two continuous fixed effects (measured at the individual level but
> containing info on 23 species), and a single random effect of
> "species". My model is mixing well and the results fit the data, but
> what I am really interested in is getting estimates for the fixed
> effect by species, rather than globally.
>
> It is my understanding (per a 2016 reply to a post to this list:
> parameter estimates for all factor levels MCMCglmm) that by
> suppressing the intercept here, I am estimating the first "level" (one
> species?) for each effect, and that the remaining effects are
> deviations from these levels. How do I code random slopes and
> intercepts for species within each of the fixed effects, rather than
> deviations from a common distribution?
>
> priorS = list(R = list(V = 1, nu = 0, fix = 1), G = list(G1 = list(V
> = 1, nu = 0.002)))
>
> smodel <-MCMCglmm(fl~ -1 + Zdbh + Znnd,
> random = ~ species,
> family = "categorical", verbose=F, pr = TRUE,
> start=list(QUASI=FALSE),
> data=IHF,prior=priorS,
> nitt=500000,burnin=5000,thin=100)
>
>
> A second part of the question is that I would like to include a random
> effect to measure the influence of phylogenetic autocorrelation. I
> know how to code this using a distance matrix (see below), however, I
> am concerned that statistically, "species" and "phylo" may not have
> enough information to partition the variance in a meaningful way.
> Perhaps running three models (random=species, random=phylo,
> random=species + phylo) and comparing goodness-of-fit could be useful?
> I would greatly appreciate any insight. Many thanks!
>
> prior = list(R = list(V = 1, nu = 0, fix = 1), G = list(G1 = list(V =
> 1,nu = 0.002),
> G2=list(V=1,nu=0.002)))
>
> model <-MCMCglmm(fl ~ Zdbh_mm + Znnd,
> random = ~species + phylo,
> family = "categorical", verbose=F, pr = TRUE,
> start=list(QUASI=FALSE),
> ginverse=list(phylo=inv.phylo$Ainv),data=IHF,prior=prior,
> nitt=500000,burnin=5000,thin=100)
> Andrea Pilar Drager
> PhD. student
> Ecology and Evolutionary Biology, Rice University
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.
More information about the R-sig-mixed-models
mailing list