Bayesian phylogenetic regression incorporating phylogenetic uncertainty by sampling from multiple trees (work in progress)

phylostan_Ktrees.R

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

# load data
data(geospiza)
dat <- geospiza$geospiza.data

# create fake sample of trees
tr <- drop.tip(geospiza$geospiza.tree, 'olivacea')
trees <- vector('list', 15)
for(i in 1:15) trees[[i]] <- tr

# get the inverse of the phylovcv, order all to match data
invA <- lapply(trees, function(x) solve(vcv.phylo(x))[row.names(dat), row.names(dat)])

stdat <- list(
	N=dim(dat)[1],
	y=dat$wingL,
	x=dat$tarsusL,
	invA=invA,
	Ntree=length(invA), #number of trees
	p=rep(1/length(invA), length(invA)) #probability of sampling each tree
)

phymodel <- "
data {
	int<lower=0> N; //species
	int<lower=0> Ntree; //number of trees
	vector[N] y; 
	real x[N];
	matrix[N,N] invA[Ntree]; //array of inverse of phylovcv?
}

parameters {
	real alpha; //intercept
	real beta; //slope
	real<lower=0> tau; // scaling factor
	simplex[Ntree] theta; //mixture sampling?
}

transformed parameters{
	real sigma; //regression error
	sigma <- 1/sqrt(tau);
}

model {
	vector[N] mu; //multivariate normal mean

	//priors
	alpha ~ normal(0, 10);
	beta ~ normal(0, 10);
	sigma ~ uniform(0,1000);
	tau ~ gamma(1,1);
	theta ~ uniform(0,1);
	
	for(n in 1:N){
		mu[n] <- alpha+beta*x[n];
	  }
	y ~ multi_normal_prec(mu, tau*invA[Ntree]);	
	
}
"

phystan = stan(model_name="phylogenetic regression", model_code = phymodel, data=stdat , iter = 1000, chains = 2)