Reposting for readability (thanks to notice from Jarrod Hadfield).
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 (yet hopefully valid)
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) = Beta0_j + Beta1_j * CON_jkl + Beta2_j * HET_jkl + phi_k
The term phi_k represents a random effect for sampling plot in a frequentist
sense, and has a distribution N(0,Sigma_phi_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):
Beta_mj = Gamma_m0 + Gamma_m1 * SPTrt_j + SIGMA_Beta
Species-specific coefficients (Beta_mj) are modeled using a multivariate normal
distribution. SIGMA_Beta is the variance-covariance matrix of the Beta's and it
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 (eta_j) at individual-level
regresion:
logit(p_jkl) = Beta0_j + Beta1_j * CON_jkl + Beta2_j * HET_jkl + phi_k + eta_j
with eta_j distributed as a multivariate normal with mean 0 and variance
covariance matrix delta2*SIGMA_Phylo, with SIGMA_Phylo derived from the
phylogenetic tree [from refs: delta2 can be interpreted as the residual
variance and SIGMA_Phylo as a correlation structure when SIGMA_Phylo 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
(SIGMA_phylo) from function in package 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 SIGMA_phylo
# 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,
eta_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) 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]
} # U is a matrix containing intercept and trait variables
}
...
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
_______________________________________________
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/
