Dear all,
I am plant ecologist using a 2-level hierarchical Bayesian model (inplemented
using JAGS and R2Jags) to study how plant functional traits affect the
sensibility of species to the
effects of neighboring trees on survival. I would like to extend my model to
account
for phylogenetic dependence and evaluate whether results I have found are due
to trait
functionality or merely to evolutionary relatedness (assuming a Brownian model
of evolution). After some reading (e.g. Freckleton et al 2002, Blomberg et al
2012; de Villemereuille et al. 2012, Abadi et al. 2014), and taking the most
simple approach, I think I can do so by adding to the model a species random
term with a variance covariance matrix derived from the phylogenetic tree. I am
doubting though whether the term should be added to the individual (first
level) or the species (second level)
regression and, due to my novelty to phylogenetic and specially Bayesian
methods, whether I am approaching the problem (and implementing a solution)
correctly.
A more specific doubt is how to scale to height 1 the variance-covariance
matrix derived from the phylogenetic tree. I would greatly appreciate any input
from the specialized community.
Below I describe the general model followed by the addition of the
phylogenetic component (I am copying only parts of my code as the mail is
already long but I can provide more details if needed or I have not expressed
myself clearly).
# SURVIVAL model WITHOUT phylogenetic dependence
In the individual (first-level) regression of the model, survival
(S) of an individual i, of species j, in plot k is modeled as a Bernoulli
process and a function of the densities
of conspecific (CON) and heterospecific neighbors (HET):
S_ijkl = Bernoulli (pjkl)
logit(pjkl) = â0j + â1j * CON_jkl + â2j * CON_jkl + ök
The term ök represents a random effect for sampling plot in a
frequentist sense, and has a distribution N(0,óö2).
In the species (second-level) regression of the model, the
first-level regression coefficients for each species j are modeled as a linear
function of species traits (SPTrt):
âmj = ãm0 + ãm1 * SPTrtj + Óâ
Species-specific coefficients (âmj) are modeled using a
multivariate normal distribution. The variance-covariance matrix of the â's,
ÓB, is modeled using a scaled inverse-Wishart distribution following Gelman and
Hill (2007).
# SURVIVAL model WITH phylogenetic relatedness (with details on
implementation in R)
# Option 1: Model with species random term (çj) at
individual-level regresion:
logit(pjkl) = â0j + â1j * CON_jkl + â2j * CON_jkl + ök + çj
with çj distributed as a multivariate normal with mean 0 and
variance covariance matrix ä2Ó, with Ó derived from the phylogenetic
tree (ä2 can be interpreted as the residual variance and Ó as a
correlation structure when Ó is scaled to a height of 1 - de Vemeuirille et al
2012, Abadi et al. 2014)
# Implementation
# First estimate the vcv matrix:
pvarcov <- vcv.phylo(SMTree) # phylogenetic variance covariance
matrix (Ó) from function in pckge ape
pvarcov175<-pvarcov[sort.list(as.numeric(rownames(pvarcov))),
sort.list(as.numeric(colnames(pvarcov)))] #order species
Omega <- solve(pvarcov175[1:J,1:J])# inverse of Ó
# The regression model (partial R/ JAGS code)
Full.poly.phy<-function (){
for (i in
1:n){
y[i]
~ dbern (p.bound[i])
p.bound[i]
<- max(0, min(1, p[i]))
logit(p[i])
<- inprod(B[SPECIES[i],],X[i,]) + a.plot[PLOT[i]] + tau.phy*epspe[SPECIES[i]]
}
...
# Definition and priors of the phylogenetic dependence species
random term, çj (i.e. tau.phy*epspe[SPECIES[i]]):
epspe[1:J] ~ dmnorm(mophy[], Omega[,]) # Omega is the inverse of
the covariance matrix obtained from the phylogeny data
for(j in 1:J) { mophy[j] <- 0 }
tau.phy <-
pow(sigma.phy,-2)
sigma.phy ~
dunif(0,100)
# Option 2, with the phylogenetic
term at the species level regresion:
Full.poly.phy<-function (){
for (i in
1:n){
y[i]
~ dbern (p.bound[i])
p.bound[i]
<- max(0, min(1, p[i]))
logit(p[i])
<- inprod(B[SPECIES[i],],X[i,]) + a.plot[PLOT[i]]
}
...
# Species specific intercepts and slopes
for (j in 1:J){
for
(k in 1:K){
B[j,k]<-xi[k]*B.raw[j,k]
}
B.raw[j,1:K]
~ dmnorm (B.raw.hat[j,], Tau.B.raw[,])#
}
# Species level regression (intercept and slopes as a function of
traits (U)) with phylogenetic random term
for (j in
1:J){
for
(k in 1:K){
B.raw.hat[j,k]
<- inprod(G.raw[k,], U[j,]) + tau.phy*epspe[j]
}
}
...
epspe[1:J] ~ dmnorm(mophy[], Omega[,]) # inverse of the covariance
matrix
for(j in 1:J) { mophy[j] <- 0 }
tau.phy <-
pow(sigma.phy,-2)
sigma.phy ~
dunif(0,100)
Many thanks in advance for any input,
Edwin
[[alternative HTML version deleted]]
_______________________________________________
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
