## gamma model for microstegium biomass with negative binomial fecundity predictions 

massdata = read.csv("garden.csv")
seed_frame = read.csv("light_exp_seed.csv")
seed_data = seed_frame


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)), seeds = seed_data$seeds[1:24], biomass = seed_data$biomass[1:24])
params = list("a", "b", "alpha2", "alpha1", "alpha", "pseeds")


sink("ranmodelgroE.txt") #the sink function drops the code as a text file into your working directory
cat('
        
    # Biomass Likelihood and Biomass/Fecundity Prediction Models 
    
  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

# Now we predict seeds with negative binomial 
    pseeds[i] ~ dnbinom(ppp[i], rll[i]) 
    rll[i] <- (b1*(p[i]^b2) + 0.00001)/(alpha - 1)  # Inherits b1, b2, and alpha from seed production model 
    ppp[i] <- 1/alpha
    }  
    
    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])
    
    }
    
    # Likelihood model for seed production as a function of biomass (mesocosm dataset)
    for(k in 1:length(seeds)) {
    seeds[k] ~ dnbinom(pp[k], ll[k])
   # nmu[k] <- (ll[k]*(1-pp[k]))/pp[k]
  #  s[k] <- alpha*nmu[k]
    
    pp[k] <- 1/alpha
    ll[k] <- (b1*(biomass[k]^b2))/(alpha - 1)
    
    }
    
    # Priors for gamma models of biomass 

    alpha1 ~ dunif(1, 10)
    alpha2 ~ dunif(1, 10)
    
     # Priors for Michaelis-Menten 
     # based on sampling distribution from mesocosm 

    a ~ dnorm(18.216, (1/2.3^2))
    b ~ dnorm(40.434, (1/10.496^2))    

    # Priors for negative binomial seed model 
    
    b1 ~ dunif(100,300)
    b2 ~ dunif(1,2)
    alpha ~ dunif(1, 1000)
    
    
    }
    
    ',fill=TRUE)
sink()

set.seed(101)

ran.MvMM1001 = 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 the gipar function used below 
pars = gipar(ran.MvMM1001) # gipar function extracts parameters from the JAGS object 

seed_preds = pars$pseeds[,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_seeds = matrix(x, nrow = 20, ncol = 100)
mean_seeds = rep(0, times = 100)
median_seeds = 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(seed_preds[,j])
    pred_seeds[i,j] = c[prob[i]]
    mean_seeds[j] = mean(c)
    median_seeds[j] = median(c)
  }
}


central_dat = data.frame(mean_seeds, median_seeds, light = mass_data$light[22:121])

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

# below is code for plotting the prediction intervals with ggplot 

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

light = seq(1, 100, 1)
pred_contours = data.frame(pred_seeds)
pred_contour1 = t(pred_contours)
argh = melt(pred_contour1)
light_vals = rep(light, times = 20)
pred_contours = data.frame(argh, light_vals) #note: be sure and check the name 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")
}

contour_preds = contour_plot + theme_bw() + labs(x = "% available light", y = "seeds/plant") 

contour_preds + stat_smooth(data = central_dat, aes(x = light, y = mean_seeds), se = FALSE, colour = "Black")


## Posterior inference on seed probabilities
## note that due to monte carlo error, probability estimates
## vary slightly between model runs 


pars <- gipar(ran.MvMM1001)

# 2 % 
post_2cent = data.frame(preds = sort(pars$pseeds[,383]))
str(post_2cent)
zeros = subset(post_2cent, preds == 0)
sufficient = subset(post_2cent, preds >= 40)

hist_2cent = ggplot() + geom_histogram(aes(x = preds), data = post_2cent, binwidth = 200)
hist_2cent + theme_bw() + labs(x = "seeds/plant", y = "count")

#result: 673 zeros. So probability of no seeds at 2% light is 0.673
# 260 observations are greater than or equal to forty. so probability is 0.260 


# 5% 
post_5cent = data.frame(preds = sort(pars$pseeds[,386]))
zeros_5 = subset(post_5cent, preds == 0)
sufficient_5 = subset(post_5cent, preds >= 40)
str(sufficient_5)

hist_5cent = ggplot() + geom_histogram(aes(x = preds), data = post_5cent)
hist_5cent + theme_bw() + labs(x = "seeds/plant", y = "count")

## Results. probability of 0 seeds is 0.366 
# probability of exceeding 40 seeds is 0.501 

# 10% 
post_10cent = data.frame(preds = sort(pars$pseeds[,391]))
zeros_10 = subset(post_10cent, preds == 0)
sufficient_10 = subset(post_10cent, preds >= 40)
str(sufficient_10)


hist_10cent = ggplot() + geom_histogram(aes(x = preds), data = post_10cent)
hist_10cent + theme_bw() + labs(x = "seeds/plant", y = "count")

# probability of 0 seeds is 0.188
# probability of 40+ seeds is 0.706 

# 80% 

post_80cent = data.frame(preds = sort(pars$pseeds[,461]))
zeros_80 = subset(post_80cent, preds == 0)

sufficient_80 = subset(post_80cent, preds >= 40)
str(sufficient_80)

hist_80cent = ggplot() + geom_histogram(aes(x = preds), data = post_80cent)
hist_80cent + theme_bw() + labs(x = "seeds/plant", y = "count")


# probability of zero seeds is 0.005
# probability of 40+ seeds is 0.987