Hi,
Regarding the blog and the feasibility of MCMCglmm for threshold models:
If y1 is binary and y2 is normal, then the univariate analysis would be:
Ainv<-inverseA(tree)$Ainv
m1<-MCMCglmm(y1~y2, random=~species,ginverse=list(species=Ainv),
data=my.data, prior=my.prior, family="ordinal")
family="ordinal" uses probit link (the original threshold model),
family="categorical" uses logit link.
The bivariate analysis would be
m2<-MCMCglmm(cbind(y2,y1)~trait,
random=~us(trait):species,rcov=us(trait):units,
ginverse=list(species=Ainv), data=my.data, prior=prior,
family=c("ordinal", "gaussian"))
here 2x2 covariance matrices are set up at the phylogenetic and
residual level. If we have these matrices as Va and Ve, then the
assumption of the univariate model is that Va[2,1]/Va[1,1] =
Ve[2,1]/Ve[1,1]. Alternative covariance structures could be used, for
example ~idh(trait):species sets Va[2,1]=0 and ~idh(trait):units sets
Ve[2,1]=0. Other less useful covariance functions for this type of
model might be ~idv(trait):species (Va[2,1]=0, Va[1,1]=Va[2,2]) and
~species (Va[2,1]=Va[1,1]=Va[2,2]).
The prior is important in models with discrete characters, because the
residual variation is not identifiable. I fix it a priori at one, for
example in the univariate model:
prior$R=list(V=1, fix=1)
and the bivariate model
prior$R=list(V=diag(2), fix=2, nu=1.002)
although you need to be more careful here, because the residual
variation in y2 is identifiable and may be sensitive to the prior
(inverse gamma in this instance). Note the ordering of the response
(y2,y1) so that fix=2 fixes Ve[2,2]=1 (the residual variance of the
binary character).
The choice of fixing it at one is entirely arbitrary (most software
chooses zero). Different values do change the estimates, but they can
always be rescaled to what would have been estimated given another
choice (see CourseNotes).
I see there was some discussion of what to do when y1 has more than 2
states. family="ordinal" takes any number of states which are taken to
be ordered. family="categorical" takes any number of states which are
not ordered. "categorical" with probit link could be implemented
(literally change plogis to pnorm at the 2/3 places it appears in the
C code) but my guess is the answers will be indistinguishable in most
cases.
There are some drawbacks to modelling dependencies through correlated
overdispersion (residuals) for non-gaussian data. However, I think
these drawbacks do not apply to discrete data where overdispersion
does not exist.
With regards to the earlier posts on asymmetric transitions, if a
probit link is used pnorm(0,intercept, sqrt(1+Va+Ve)) is
approximately q01/(q10+q01) from the MK models. I have not been
successful at finding an equivalence for q10+q01.
Cheers,
Jarrod
Quoting Joe Felsenstein <j...@gs.washington.edu> on Wed, 29 Aug 2012
12:41:48 -0700:
Theodore Garland Jr wrote:
Check this:
Ives, A. R., and T. Garland, Jr. 2010. Phylogenetic logistic
regression for binary dependent variables. Systematic Biology
59:9-26.
In addition, check this:
Felsenstein, J. 2012. A comparative method for both discrete and
continuous characters using the threshold model. American Naturalist
179: 145-156.
Joe
----
Joe Felsenstein j...@gs.washington.edu
Dept of Genome Sciences and Dept of Biology, Univ. of Washington,
Box 5065, Seattle Wa 98195-5065
_______________________________________________
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.
_______________________________________________
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
