Bayesian phylogenetic regression

phylostan_1tree.R

require(rstan)
require(geiger)
require(MCMCglmm)

# load data
data(geospiza)

# get the inverse of the phylogenetic variance-covariance matrix
tr <- drop.tip(geospiza$geospiza.tree, 'olivacea')
invA <- solve(vcv.phylo(tr))
# Alternatively, using MCMCglmm - not sure if different, results seem the same
# invA <- as.matrix(inverseA(tr, 'TIPS')$Ainv)

# get the data, order to match invA
dat <- geospiza$geospiza.data
dat <- dat[row.names(invA),]

# prepare the data in STAN format
stdat <- list(
  N=dim(dat)[1],
	y=dat$wingL,
	x=dat$tarsusL,
	invA=invA
)

# STAN model
phymodel <- "
data {
	int<lower=0> N; //Number of species
	vector[N] y; //y variable (vector because of multivariate error)
	real x[N]; //x variable
	cov_matrix[N] invA; //inverse of the phylovcv matrix
}

parameters {
	real alpha; //intercept
	real beta; //slope
	real<lower=0> sigma; //regression error
}

transformed parameters{
	real tau;
	tau <- pow(1/sigma, 2); //scaling factor for the precision matrix
}

model {
	vector[N] mu; //multivariate normal mean
	//priors
	alpha ~ normal(0, 1000);
	beta ~ normal(0, 1000);
	sigma ~ uniform(0,1000);
	tau ~ gamma(1,1);

	for(n in 1:N){
	  mu[n] <- alpha+beta*x[n];
	  }
	y ~ multi_normal_prec(mu, tau*invA); //precision matrix error
}
"

#simple run, compiles OK, results seem fine
#lots of errors of non-positive-definite invA matrix (as also happens in MCMCglmm)
#ignore?
phystan <- stan(model_name="phylogenetic regression", model_code = phymodel, 
  data=stdat , iter = 1000, chains = 2)