library(rstan)
devtools::source_gist("8f675b9aa74ab9900fb4")
load("Sig_test_DJ.RData")
data = all_model
data$Year = factor(data$Year)
data = droplevels(data)

### Add beach variables
data$Beach4 = ifelse(data$Beach=="4", 1, 0)
data$Beach7 = ifelse(data$Beach=="7", 1, 0)
data$BeachN = ifelse(data$Beach=="N", 1, 0)
data$BeachS = ifelse(data$Beach=="S", 1, 0)
data$IslandM = ifelse(data$Island=="MARMOT", 1, 0)
data$IslandU = ifelse(data$Island=="UGAMAK", 1, 0)

sat_model = list(
  phi1 = ~ (Beach-1),
  phi2 = ~ (Beach-1),
  phi3 = ~ (Beach-1),
  sigma_group =  ~ (Year-1):(Beach-1)
)

fit_inits = fit_logistic(
  formula=m_pups~M1, 
  model_list=sat_model, 
  data=data
)


data_list = list(
  N = nrow(data),
  n_year = length(levels(data$Year)),
  y = data$m_pups,
  x = data$M1,
  beach = as.integer(data$Beach),
  year = as.integer(data$Year)
)

inits = function(){
  list(
    beta1=fit_inits$results$beta$phi_1$beta,
    beta2=fit_inits$results$beta$phi_2$beta,
    beta3=fit_inits$results$beta$phi_3$beta,
    re1 = matrix(0, 10, 4),
    re2 = matrix(0, 10, 4),
    re3 = matrix(0, 10, 4),
    tau1=rep(1,4), tau2=rep(1,4), tau3=rep(1,4)
  )
}

fit = stan(file="logistic_re.stan", data=data_list, init=inits, chains=1)

smp_list = extract(fit)

### Get mean date for each beach/year
phi2_beach4 = smp_list$beta2[,1]
phi1_beach4 = exp(smp_list$beta1[,1])
phi2_beach7 = smp_list$beta2[,2]
phi1_beach7 = exp(smp_list$beta1[,2])

### Posterior summary
summary(mcmc(phi2_beach4))
HPDinterval(mcmc(phi2_beach4))

### Beach effect at Marmot
marmot_beach_eff = phi2_beach7-phi2_beach4
summary(mcmc(marmot_beach_eff))
HPDinterval(mcmc(marmot_beach_eff))
# Beach 7 has a mean date of birth 2-5.5 days later than beach 4!

### Calculate marmot "average"
phi2_marmot = (phi2_beach4*phi1_beach4 + phi2_beach7*phi1_beach7)/(phi1_beach4 + phi1_beach7)
summary(mcmc(phi2_marmot))
HPDinterval(mcmc(phi2_marmot))