# R coding for the Bayes A approach of genomic selection

nmarkers <- 10 #number of markers
nrecords <- 50 #number of records
numit <- 1000 #number of iterations

x<-read.table('/Users/guilherme/Dropbox/2016/teaching/Quantitative Genetics/R programs/genomic selection/reference.csv',sep=',',header=TRUE)

y<-as.matrix(x[,1])
x<-as.matrix(x[,2:11])

# Setting up storage vectors and matrices for MCMC output
gStore <- array(0,c(numit,nmarkers))
gvarStore <- array(0,c(numit,nmarkers))
vareStore <- array(0,c(numit))
muStore <- array(0,c(numit))
ittstore <- array(0,c(numit))

# Gibbs sampling

# Starting values and some additional arrays
g <- array(0.01,c(nmarkers))
mu <- 50
gvar <- array(0.1,c(nmarkers))
ones <- array(1,c(nrecords))
e <- array(0,c(nrecords))

for (l in 1:numit) {

# Sampling vare
e <- y - x%*%g - mu
vare <- (t(e)%*%e)/rchisq(1,nrecords-2)

# Sampling the mean
mu <- rnorm(1,(t(ones)%*%y -t(ones)%*%x%*%g)/nrecords,sqrt(vare/nrecords))

# Sampling gvar
for (j in 1:nmarkers) {
# gvar[j] <- (0.002+g[j]*g[j])/rchisq(1,4.012+1) # Meuwissen et al. (2001) prior
# gvar[j] <- (0.002+g[j]*g[j])/rchisq(1,1) # Xu (2003) prior
gvar[j] <- (0.002+g[j]*g[j])/rchisq(1,0.998) # Te Braak et al. (2006) prior
}

# Sampling the marker effects
z <- array(0,c(nrecords))
for (j in 1:nmarkers) {
gtemp <- g
gtemp[j] <- 0
for (i in 1:nrecords) {
z[i] <- x[i,j]
}
mean <- ( t(z)%*%y-t(z)%*%x%*%gtemp-t(z)%*%ones*mu ) /
(t(z)%*%z+vare/gvar[j])
g[j] <- rnorm(1,mean,sqrt(vare/(t(z)%*%z+vare/gvar[j])))
}

# store the sampled values
for (j in 1:nmarkers) {
gStore[l,j] <- g[j]
gvarStore[l,j] <- gvar[j]
}
vareStore[l] <- vare
muStore[l] <- mu
ittstore[l] <- l
}

# For visual inspection of convergence of MCMC for each parameter
# plot(1:1000,gStore[1:1000,1],type='l')

plot(1:1000,vareStore,type='l')
plot(vareStore[1:100],type='l')

# To look at the posterior density of each parameter
# plot(density(gStore[101:1000,1]))
plot(density(vareStore[101:1000]))

# To get the posterior means of each parameter (e.g. gene effects)
# mean(gStore[101:1000,1])