massdata = read.csv("garden.csv")
sites = massdata$Site   # This code cleans up site as a vector for the biomass model 
for (i in 1:length(sites)) {
  if (sites[i] >= 10) {
    sites[i] = sites[i] - 1
  } 
}

massdata = data.frame(sites, massdata)
light.dat<-unique(data.frame(massdata$Light,massdata$sites))

library("R2jags")


mass_data = list(mass = c((massdata$totmass), rep(NA, times = 100)), light = c(light.dat$massdata.Light, seq(1,100,1)), site = c(massdata$sites, seq(22,121,1)))
params = list("a", "b", "mu", "p", "alpha1", "alpha2") 

sink("ranmodelgroE.txt") #the sink function drops the code as a text file into your working directory
cat('
    
    model{
    for(i in 1:481) {
    mass[i] ~ dgamma(rr[i], lambda[i])
    rr[i] <- mu[site[i]]/alpha1
    lambda[i] <- 1/alpha1
    p[i] ~ dgamma(rrr[i], lambda[i])
    rrr[i] <- nmu[site[i]]/alpha1

    }  
    
        for(j in 1:121) {
         mu[j] ~ dgamma(nr[j], nlambda[j])
         nr[j] <- ((a*light[j])/(b + light[j]))/alpha2
          nlambda[j] <- 1/alpha2
        nmu[j] ~ dgamma(nr[j], nlambda[j])

        }
    
    
    alpha1 ~ dunif(1, 10)
    alpha2 ~ dunif(1, 10)


# Sampling distribution priors 
# question is: does information from mesocosm improve the model? 

    a ~ dnorm(18.216, (1/2.3^2))
    b ~ dnorm(40.434, (1/10.496^2))
    
    
    }
    
    ',fill=TRUE)
sink()

set.seed(101)

ran.MvMM101 = jags(mass_data ,inits=NULL,unlist(params),model.file="ranmodelgroE.txt",n.chains=5,n.iter=12000,n.burnin=2000,n.thin= 50)

source("generalBAYESIANcode.R") # contains function gipar used below 

pars = gipar(ran.MvMM101)
mass_preds = pars$p[,382:481]   # Note there are 1000 samples from the posterior


## Each column in the matrix represents predictions for a different biomass level 
## each row is a sample from the posterior (n = 1000)

x = rep(0, times = 20*100)
pred_mass = matrix(x, nrow = 20, ncol = 100)
mean_mass = rep(0, times = 100)
median_mass = rep(0, times = 100)
prob = rep(0, times = 100)

for (i in 0:19) {
  prob[i] = 25 + 50*i 
  for (j in 1:100) {
    c = sort(mass_preds[,j])
    pred_mass[i,j] = c[prob[i]]
    mean_mass[j] = mean(c)
    median_mass[j] = median(c)
  }
}

## The above represents data sampled for the 2.5-97.5 percentile prediction interval

# below is code for plotting the posterior predictions against data


library("ggplot2")
library("reshape2")

cen_data = data.frame(mean_mass, median_mass, light = seq(1,100,1))

light = seq(1, 100, 1)

pred_contours = data.frame(pred_mass)
pred_contour1 = t(pred_contours)
argh = melt(pred_contour1)
light_vals = rep(light, times = 20)
pred_contours = data.frame(argh, light_vals) # check name generated for Var2 used below

contour_plot = ggplot(pred_contours, aes(x = light_vals))
for (i in 1:19) {
  pred_contours_1 = subset(pred_contours, Var2 == i)
  contour_plot = contour_plot + stat_smooth(data  = pred_contours_1, aes(y = value), se = FALSE, linetype = "dashed")
}

biomass_fig = contour_plot + geom_point(aes(x = Light, y = totmass), data = massdata) + stat_smooth(aes(x = light, y = mean_mass), data = cen_data, se = FALSE, colour = "Red")
biomass_fig + theme_bw() + labs(x = "% available light", y = "biomass/plant (g)") + scale_y_continuous(breaks = seq(5, 50, 5))