##load packages needed:
library(lme4)#
library(brms)
library(HDInterval)
##function for assessing the degree of overdispersion; taken from https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html:
overdisp_fun <- function(model) {
    rdf <- df.residual(model)
    rp <- residuals(model,type="pearson")
    Pearson.chisq <- sum(rp^2)
    prat <- Pearson.chisq/rdf
    pval <- pchisq(Pearson.chisq, df=rdf, lower.tail=FALSE)
    c(chisq=Pearson.chisq,ratio=prat,rdf=rdf,p=pval)
}
##set working directory:
setwd("path")
all.data=read.table(file="data.txt", header=T, sep="\t", stringsAsFactors=T)
##preparatory steps:
all.data$total.nr.in.area=round(78.3*all.data$nr.individuals/2)
all.data=subset(all.data, nr.individuals>0)
##checks:
range(all.data$nr.hosts/all.data$total.nr.in.area)
##exclude unidentified genera:
nrow(all.data)
all.data=droplevels(subset(all.data, Host.genus!="indet. sp. 1" & Host.genus!="indet . sp. 2"))
nrow(all.data)
##some cleaning:
all.data$Host.genus=as.factor(gsub(x=as.character(all.data$Host.genus), pattern="Parkia ", replacement="Parkia", fixed=T))
##add response variable:
all.data$rv=cbind(all.data$nr.hosts, all.data$total.nr.in.area)
##fit miced model (to assess dispersion parameter):
m.f=glmer(rv~1+(1|Host.genus), family=binomial, data=all.data)
overdisp_fun(m.f)

##fit Bayesian model:
m.b=brm(nr.hosts|trials(total.nr.in.area)~1+(1|Host.genus), family=binomial, data=all.data, 
  prior=set_prior("student_t(4, 0, 2.5)", class = "sd", group="Host.genus", coef = "Intercept"),
  chains=10, iter=5000, warmup=2500, control=list(adapt_delta=0.95, max_treedepth=20))
##extract posterior samples:
ests=posterior_samples(m.b)
##extract results:
##link space standard deviation of infestation probabilities per genus:
paste(round(as.vector(hdi(ests[, "sd_Host.genus__Intercept"], credMass=0.99)), 3), collapse=" to ")
##prepare for plot:
pred.mat=as.vector(ests[, "b_Intercept"])+as.matrix(ests[, grepl(x=colnames(ests), pattern="r_Host.genus", fixed=T)])
xnames=colnames(ests)[grepl(x=colnames(ests), pattern="r_Host.genus", fixed=T)]
xnames=gsub(xnames, pattern="r_Host.genus[", replacement="", fixed=T)
xnames=gsub(xnames, pattern=",Intercept]", replacement="", fixed=T)
colnames(pred.mat)=xnames
all.data$prop.with=all.data$nr.hosts/all.data$total.nr.in.area
all.data=all.data[rev(order(all.data$prop.with)), ]
spec.order=rev(sort(apply(pred.mat, 2, median)))
xorder=match(names(spec.order), as.character(all.data$Host.genus))
##plot:
dev.off()
X11(width=7, height=10)
par(mar=c(0.2, 5, 3, 0.2), mgp=c(1.7, 0.3, 0), tcl=-0.15, las=1)
plot(1, 1, ylim=c(nrow(all.data)+1, 0), yaxs="i", xlim=c(0, 0.07), yaxt="n", xlab="", ylab="", xaxt="n")#xlim=range(c(all.data$prop.with, plogis(pred.mat)))
mtext(text="probability [%]", side=3, line=1.2)
axis(side=3, at=seq(0, 0.12, by=0.02), labels=seq(0, 0.12, by=0.02)*100)
mtext(text=paste(all.data$Host.genus[xorder], " (", all.data$nr.hosts[xorder], "/", all.data$total.nr.in.area[xorder], ")", sep=""), side=2, at=1:nrow(all.data), line=0.2, adj=1, cex=0.5)
for(i in 1:nrow(all.data)){
  xspec=gsub(names(spec.order[i]), pattern=" ", replacement=".", fixed=T)
  #xspec=gsub(as.character(all.data$Host.genus[i]), pattern=" ", replacement=".", fixed=T)
  hpdi=as.vector(hdi(pred.mat[, xspec], credMass=0.99))
  xmed=plogis(median(pred.mat[, xspec]))
  to.add=density(pred.mat[, xspec])
  to.add=cbind(x=to.add[["x"]], y=to.add[["y"]])
  to.add=to.add[to.add[, "x"]>=min(pred.mat[, xspec]) & to.add[, "x"]<=max(pred.mat[, xspec]), ]
  to.add=to.add[to.add[, "x"]>=min(hpdi) & to.add[, "x"]<=max(hpdi), ]
  to.add[, "y"]=0.6*to.add[, "y"]/max(to.add[, "y"])
  to.add[, "y"]=i-to.add[, "y"]
  to.add.99=to.add
  to.add.99[, "x"]=plogis(to.add.99[, "x"])
  to.add.99=rbind(to.add.99[1, ], to.add.99, to.add.99[nrow(to.add), ])
  to.add.99[c(1, nrow(to.add.99)), "y"]=i
  #browser()
  hpdi=as.vector(hdi(pred.mat[, xspec], credMass=0.5))
  to.add=to.add[to.add[, "x"]>=min(hpdi) & to.add[, "x"]<=max(hpdi), ]
  to.add[, "x"]=plogis(to.add[, "x"])
  to.add=rbind(to.add[1, ], to.add, to.add[nrow(to.add), ])
  to.add[c(1, nrow(to.add)), "y"]=i
  polygon(to.add, border=NA, col=grey(level=0.8))
  polygon(to.add.99)
  to.add=which(all.data$Host.genus==xspec)
  points(x=all.data$nr.hosts[to.add]/all.data$total.nr.in.area[to.add], y=i, pch=1, cex=0.5)
  points(x=xmed, y=i, pch=4, cex=0.5)
}
abline(v=sum(all.data$nr.hosts)/sum(all.data$total.nr.in.area), lty=3)
dev.copy2pdf(file="FIG.pdf")
savePlot(file="FIG.png", type="png")
##prepare Table:
xx.99=apply(pred.mat, 2, hdi, credMass=0.99)
xx.50=apply(pred.mat, 2, hdi, credMass=0.5)
xx.med=apply(pred.mat, 2, median)
xx=cbind(xx.99[1, ], xx.50[1, ], xx.med, xx.50[2, ], xx.99[2, ])
colnames(xx)=c("lowerCl99", "lowerCl50", "median", "upperCl50", "upperCl99")
xx=plogis(xx)

save.image(file="workspace.RData")