library(ape)
library(phylobase)
library(MCMCglmm)
library(parallel)
library(coda)
library(lme4)

# -----------------------------------------------
# - Read maximum credibility phylogenetic tree --
# -----------------------------------------------
fk <- read.nexus("phy.nex")

# -----------------------------------------------
# -Log10-transformation of continuous variables -
# -----------------------------------------------
# read data table
table(fid$Farm)
fid$Start.dist <- log10(fid$Start.dist)
fid$Flock <- log10(fid$Flock)
fid$FID <- log10(fid$FID)

# -----------------------------------------------
# --------------- MCMCglmm model ----------------
# -----------------------------------------------
inv.f <- inverseA(fk, nodes="ALL", scale=TRUE) # Inverse matrix of phylogenetic relatedness
prior1 <- list(R=list(V=1e-8, nu=-2), G=list(G1=list(V=diag(1), nu=0.002))) # Prior settings
mcmc1<-MCMCglmm(FID ~ Farm + Flock + Start.dist, random=~Species, family="gaussian", ginverse=list(Species = inv.f$Ainv), 
prior=prior1, verbose=F, nitt=1000000, burnin=20000, thin=100, data=fid)
summary(mcmc1)
plot(mcmc1)
lambda <- mcmc1$VCV[,"Species"] / (mcmc1$VCV[,"Species"] + mcmc1$VCV[,"units"])
lambda <- posterior.mode(lambda)
summary(lambda)

# -----------------------------------------------
# ----------- Gelman-Rubin statistic ------------
# -----------------------------------------------
mcmcG <- mclapply(1:4, function(i) {
  MCMCglmm(FID ~ Farm + Flock + Start.dist, random=~Species, family="gaussian", ginverse=list(Species = inv.f$Ainv), 
prior=prior1, verbose=F, nitt=1000000, burnin=20000, thin=100, data=fid)
}, mc.cores=4)
mcmcG <- lapply(mcmcG, function(m) m$Sol)
mcmcG <- do.call(mcmc.list, mcmcG)
gelman.plot(mcmcG, auto.layout=F)
plot(mcmcG, ask=F, auto.layout=F)
gelman.diag(mcmcG)

# -----------------------------------------------
# ------------- Linear mixed model --------------
# -----------------------------------------------
library (lme4)
library (lmerTest)
lmer1 <- lmer(FID ~ Farm + Flock + Start.dist + (1|Species), data=fid, REML=F)
summary(lmer1)