# ===Ancestral behaviors of melanogaster fly===
# =============================================
#Best place to understand code are Hadfield course Notes on GlMM chapter 5.
#https://cran.r-project.org/web/packages/MCMCglmm/vignettes/CourseNotes.pdf

# ---Load the necessary libraries

#library(MCMCglmm)
library(coda,lib.loc="Rpackages")
library(ape,lib.loc="Rpackages")
library(Matrix,lib.loc="Rpackages")
library(MCMCglmm)
library(parallel)
##library(parallel,lib.loc="/Rpackages/")


# ---Global parameters
Nb = 134 #number of behaviors

# ---Set the output file
# sink("ancestral0.out")

# ---Defining phylogeny -> pedigree variable for MCMCglmm
tt="((yakuba:0.09,santomea:0.09)oldyasa:0.44,((sechelia:0.15,simulans:0.15,mauritiana:0.15)oldsim:0.14,melanogaster:0.29)oldmela:0.24);"


flytree<-read.tree(text=tt)
# ---Plot of tree
#par(ask=TRUE)
#plot(flytree,edge.width=2,label.offset=0.1)
#nodelabels()
#tiplabels()
#edgelabels(flytree$edge.length, bg="black", col="white", font=2)
#dev.copy2eps()
#par(ask=FALSE)


# ---Load the dataset: l_i= log p(behavior_i) - log p(behavior_0) 
# 134 behaviors + 0 behavior, 593 flies from 6 species (column-> animal)
LFlydat <- read.table("../../data/LogBehaviorData.txt", header=TRUE, sep="\t", row.names="id")


# ---Defining prior
# non-informative prior
IJ <- (1/(Nb+1))*(diag(Nb)+matrix(1,Nb,Nb))
prior.1<-list(G=list(G1=list(V=IJ/2,n=Nb)), R=list(V=IJ/2,n=Nb))

# ---Creating fixed effects
sfix="cbind("
for (i in 1:(Nb-1)){sfix <-paste(sfix,colnames(LFlydat)[i+1],",",sep="")}
sfix <-paste(sfix,colnames(LFlydat)[Nb+1],") ~ trait-1",sep="")
fixed <- as.formula(sfix)

#---Run MCMCglmm
model0 <- MCMCglmm(fixed=fixed,
                   random = ~us(trait):animal,
                   rcov = ~us(trait):units,
                   data = LFlydat,
                   family = rep("gaussian", Nb),
                   prior=prior.1,
                   pedigree=flytree,
                   pr=TRUE,
                   verbose = FALSE,
                   thin=20)
save(model0, file = "Anc134BhModel0MinusMauNoYak00")
#write.table(model0$Sol, file="ModelSolAllMinSimNoYak00.txt", row.names=FALSE, col.names=FALSE)
#write.table(model0$VCV, file="ShortScriptVCVMCMC.txt", row.names=FALSE, col.names=FALSE)