data {

  int n;
  int Parasite[n];
  int Endo_ecto[n];
  int Par_group[n];
  int Year[n];
  vector[n] D;
  real Mean_abund[n];
  int n_Parasite;
  int n_Year;
  int n_Endo_ecto;
  int n_Par_group;

}

parameters {
  vector[n_Parasite] theta;
  vector[n_Endo_ecto] beta;
  vector[n_Year] lambda;
  real bar_mu; //Mean of theta random
  real<lower=2> phi;
  real<lower=0> sigma_y; //sd of year random
  real<lower=0> bar_sigma; //sd of theta random
}

transformed parameters{
  
}

model {
  // model calculations
  vector[n] mu;                        // transformed linear predictor
  vector[n] A;                         // parameter 1 for beta distn
  vector[n] B;                         // parameter 2 for beta distn
  
  // priors
  phi ~ cauchy(2, 2.5);   
  beta ~ normal(0, 2.2);
  bar_mu ~ normal(0, 1.5);
  sigma_y ~ exponential(1);
  bar_sigma ~ exponential(1);
  
  for (j in 1:n_Year) {
    lambda[j] ~ normal(0, sigma_y);
  }
  
  for (j in 1:n_Parasite) {
    theta[j] ~ normal(bar_mu, bar_sigma);
  }
  
  for (i in 1:n) { 
    mu[i] = theta[Parasite[i]] + beta[Endo_ecto[i]] * Mean_abund[i] + lambda[Year[i]];
	mu[i] = inv_logit(mu[i]);   
  }

  for (i in 1:n) { 					//conversion of mu and phi to A and B beta distribution parameterization
    A[i] = mu[i] * phi;
    B[i] = (1.0 - mu[i]) * phi;
  }

  // likelihood
  D ~ beta(A, B);
  
}

generated quantities {
  
  vector[n] log_lik;
  vector[n] end_mu;
  vector[n] end_A;
  vector[n] end_B;
  
  for (i in 1:n) {
    
    end_mu[i] = theta[Parasite[i]] + beta[Endo_ecto[i]] * Mean_abund[i] + lambda[Year[i]];
	end_mu[i] = inv_logit(end_mu[i]);
    
    end_A[i] = end_mu[i] * phi;
    end_B[i] = (1.0 - end_mu[i]) * phi;
    
    log_lik[i] = beta_lpdf(D[i] | end_A[i], end_B[i]);
  
  }
  
}