# The setwd command needs to be changed
setwd ("C:/Users/mstear/Documents/Nematodirus")

#
#install.packages("quadprog")
#install.packages("gtools")
#install.packages("MASS")
#install.packages("mvtnorm")
#install.packages("combinat)")
#install.packages("coda")
#install.packages("genetics")
#install.packages("Matrix")
#install.packages("ape")
#install.packages("kinship2")
#install.packages("MasterBayes")
#install.packages("ggplot2")
#install.packages("MCMCglmm")
#install.packages("gridExtra")



library(gtools)
library(MASS)
library(mvtnorm)
library(combinat)
library(coda)
library(genetics)
library(Matrix)
library(ape)
library(kinship2)
library(MasterBayes)
library(ggplot2)
library(MCMCglmm)
library(gridExtra)  #For putting plots side by side




################################################################################
#Read and process pedigree data.
pedigree <- read.csv("ped95.csv",header=TRUE)
pedigree <- orderPed(pedigree)


#Read and process phenotypic data
data1 <- read.csv("RevTandlemuir 1995 Nematodirus Epg-1.csv",header=TRUE)


data1$animal<-as.factor(data1$animal)
data1$dam<-as.factor(data1$dam)
data1$year<-as.factor(data1$year) 
data1$sex<-as.factor(data1$sex)
data1$field <- as.factor(data1$field)
data1$birthtyp <- as.factor(data1$birthtyp)
data1$age <-as.numeric(data1$age)
data1$nembjun <-as.numeric(paste(data1$nembjun))/12.5
data1$nembjul <-as.numeric(paste(data1$nembjul))/12.5 

###############################################################################
###############################################################################
###############################################################################

###############################################################################
#####Analysis for nembjun

table (data1$nembjun == 0)

nitt_default <- 1.5*10^5
p.var<-var(data1$nembjun,na.rm=TRUE)

#prior1.2<-list(G = list(G1 = list(V = 1, nu = 2), G2 = list(V = 1, nu = 2)), R = list(V = diag(2), n=2, fix =2)) 
 
prior1.3<-list(G = list(G1 = list(V = 1, nu = 1, alpha.mu=0, alpha.V=1), 
                        G2 = list(V = 1, nu = 1, alpha.mu=0, alpha.V=1)), 
                         R = list(V = diag(2), n=1, fix = 2))


startTime=proc.time()

model_zinb <- MCMCglmm(nembjun ~ trait -1 +
  at.level(trait,1):age + at.level(trait,1):field, 
  random = ~idh(at.level(trait,1)):animal + idh(at.level(trait,1)):dam, 
  rcov=~idh(trait):units, 
  family =  "zipoisson",
  data = data1, 
  prior = prior1.3, 
  pedigree = pedigree, 
  singular.ok=TRUE,  
  pr = TRUE, 
  saveX = TRUE, 
  saveZ = TRUE, 
  nitt = nitt_default,
  burnin = 2000, 
  thin = 100,
  verbose = FALSE)
stopTime=proc.time()
show(startTime)
show(stopTime)
elapsedTime=stopTime - startTime
show(elapsedTime)

autocorr(model_zinb$VCV[,"traitnembjun.units"], lags = c(1, 10, 50, 500), relative=TRUE)
heidel.diag(model_zinb$VCV) 

posterior.mode(model_zinb$VCV)
HPDinterval(model_zinb$VCV)
posterior.heritability1.2<-model_zinb$VCV[,"at.level(trait, 1).animal"]/(model_zinb$VCV[,"at.level(trait, 1).animal"] + model_zinb$VCV[,"at.level(trait, 1).dam"]+ model_zinb$VCV[,"traitnembjun.units"])  
HPDinterval(posterior.heritability1.2,0.95)
posterior.mode(posterior.heritability1.2)
median(posterior.heritability1.2)
mean (posterior.heritability1.2)
sd (posterior.heritability1.2)
effectiveSize(posterior.heritability1.2)
summary(model_zinb)

### added on revision following timothee botton comments
cov(model_zinb$VCV)

plot(posterior.heritability1.2) ## plot of estimating posterior heritability in each iteration
## histogram for posterior heritability
histjunallDA95 <-hist(round(posterior.heritability1.2,3), main = "Heritability", 
   xlab="Posterior Heritability", ylab = "Frequency", xlim=c(0,1.0), ylim = c(0,200), col="black", breaks=20)
 
dev.print(png, "histjunallDA95.png", res=600, units="in", width=9, height=12)  
dev.off() # use this code to close it

par("mar")
par(mar=c(2,1,1,2))
diagjunanimalDA95 <- plot(model_zinb$VCV[,"at.level(trait, 1).animal"] , main = "Animal", col="black")
 dev.print(png, "diagjunanimal2DA95.png", res=600, units="in", width=9, height=12)  
 dev.off() # use this code to close it
diagjundamDA95 <- plot(model_zinb$VCV[,"at.level(trait, 1).dam"] , main = "Dam", col="black")
 dev.print(png, "diagjundam2A95.png", res=600, units="in", width = 9, height = 12)  
 dev.off() 
diagjunresidDA95 <- plot(model_zinb$VCV[,"traitnembjun.units"] , main = "Residual", col="black")
 dev.print(png, "diagjunresid2DA95.png", res=600, units="in", width=9, height=12)  
 dev.off()