logistic_analysis_stan.R

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

logistic_re.stan

data {
int<lower=0> N;
int n_year;
vector[N] y;
vector[N] x;
int beach[N];
int year[N];
}

parameters {
vector[4] beta1;
vector[4] beta2;
vector[4] beta3;
matrix[n_year, 4] beta_s;
matrix[n_year, 4] re1;
matrix[n_year, 4] re2;
matrix[n_year, 4] re3;
vector[4] tau1;
vector[4] tau2;
vector[4] tau3;
}

transformed parameters {
  vector[N] phi1;
  vector[N] phi2;
  vector[N] phi3;
  vector[N] mu;
  matrix[n_year,4] sigma;
  for(i in 1:N){
    phi1[i] <- exp(beta1[beach[i]] + re1[year[i],beach[i]]);
    phi2[i] <- beta2[beach[i]] + re2[year[i],beach[i]];
    phi3[i] <- exp(beta3[beach[i]] + re3[year[i],beach[i]]);
    mu[i] <- phi1[i] / (1+exp(-(x[i]-phi2[i])/phi3[i]));
    }
  sigma <- exp(beta_s);
}

model {
  for(j in 1:4){
    for(i in 1:n_year){
        re1[i,j] ~ normal(0,tau1[j]);
        re2[i,j] ~ normal(0,tau2[j]);
        re3[i,j] ~ normal(0,tau3[j]);
       beta_s[i,j] ~ normal(0,100);
    }
  }
  beta1 ~ normal(0,100);
  beta2 ~ normal(0,100);
  beta3 ~ normal(0,100);
  tau1 ~ uniform(0,1000);
  tau2 ~ uniform(0,1000);
  tau3 ~ uniform(0,1000);
  for(i in 1:N){
    y[i] ~ normal(mu[i], sigma[year[i],beach[i]]);
  }
}