Summary of logistic growth MCMC output from stan

logistic_analysis_stan_summary.R

library(cowplot)
library(coda)
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)

load("logistic_mcmc_output.RData")

levels(data$Beach)
### Get phis for each beach/year
phi2_beach4 = smp_list$beta2[,1] + smp_list$re2[,,1]
phi1_beach4 = exp(smp_list$beta1[,1] + smp_list$re1[,,1])
phi3_beach4 = exp(smp_list$beta3[,1] + smp_list$re3[,,1])

phi2_beach7 = smp_list$beta2[,2] + smp_list$re2[,,2]
phi1_beach7 = exp(smp_list$beta1[,2] + smp_list$re1[,,2])
phi3_beach7 = exp(smp_list$beta3[,2] + smp_list$re3[,,2])

phi2_beachN = smp_list$beta2[,3] + smp_list$re2[,,3]
phi1_beachN = exp(smp_list$beta1[,3] + smp_list$re1[,,3])
phi3_beachN = exp(smp_list$beta3[,3] + smp_list$re3[,,3])

phi2_beachS = smp_list$beta2[,4] + smp_list$re2[,,4]
phi1_beachS = exp(smp_list$beta1[,4] + smp_list$re1[,,4])
phi3_beachS = exp(smp_list$beta3[,4] + smp_list$re3[,,4])

### Get median phis over years..
# phi2_beach4 = smp_list$beta2[,1]
# phi1_beach4 = exp(smp_list$beta1[,1])
# phi3_beach4 = exp(smp_list$beta3[,1])



### 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))

### Calc. prop. born by areal surrvey (day=54)
surv_prop_7 = 1/(1 + exp(-(54-phi2_beach7)/phi3_beach7))
summary(mcmc(surv_prop_7))
HPDinterval(mcmc(surv_prop_7))

### Calc. 90th percentile birth date
perc90_7 = phi2_beach7 - phi3_beach7*log(1/0.9 - 1)
summary(mcmc(perc90_7))
HPDinterval(mcmc(perc90_7))