Seeking posterior predictive distributions for 727 parameters (lakes) that are destinations for seaplane pilots but where a portion of the flights originates from source lakes harboring aquatic invasive species. I am modeling this like trials using a hierarchical binomial model just like the rats tumor example to keep it simple. Have data for over 25000 flights but some of the lakes only receive a handfull of flights.

Model runs fine but:

- max tree depth reached each time I increase, most recently to 15, but still same warning messages
- cannot open pairs() plot as suggested Error in plot.new(): figure margins too large
- half of the lakes have Rhat larger than 1.1

Some of the Rhat >1.1 lakes have chains that run parallel horizontally, not mixing

Increased iterations from 2000 to 10000 and most recently 20000 (took 4 days to run with multiple processors etc.) same warnings - bulk ESS too low
- tail ESS too low

Used launch_shinystan but found**no**issues with divergence. So I am wondering given the very large number of parameters/predictive distributions whether the diagnostics are to be expected and as good/simple as we can get here. The distributions will inform detection efforts to survey lakes that show high median introduction rates (tumor successes) first the ones with low variance then the ones with higher variance.

Here is the model:

```
J <- nrow(data)
n <- data$flights
y <- data$eflights # flights from invader source lakes
library(rstan)
rstan_options(auto_write = TRUE)
Sys.setenv(LOCAL_CPPFLAGS = '-march=native')
options(mc.cores = parallel::detectCores())
datafit3 <- stan(file="STAN/flights_model/flights_model.stan",data=c("J","y","n"), iter = 20000, chains = 4,control=list(adapt_delta=0.8, max_treedepth =15))
```

Stancode

```
// Following Carpenter et al. 2016 p. 16
data {
int<lower=0> J; // number of lakes
int<lower=0> y[J]; // number of flights from elodea-invaded sources into lake j
int<lower=0> n[J]; // number of total flights into lake j (number of trials)
}
parameters {
real<lower=0,upper=1> theta[J]; // chances of success
real<lower=0,upper=1> lambda; // prior mean chance of success
real<lower=0.1> kappa; // prior count
}
transformed parameters {
real<lower=0> alpha; // prior success count
real<lower=0> beta; // prior failure count
alpha = lambda* kappa; //alpha parameter for the prior
beta = (1 - lambda)* kappa; //beta parameter for the prior
}
model {
lambda ~ uniform(0,1); // hyperprior
kappa ~ pareto(0.1,1.5); // hyperprior was uniform(0.1,5) for results as of 04/30/20 creating datafit1, pareto(0.1,1.5) creating datafit2, uniform(0,1) in Carpenter et al.
theta ~ beta(alpha,beta); // prior for the probability p that lake j is invaded
y ~ binomial(n,theta); // likelihood for flights
}
generated quantities {
real<lower=0,upper=1> avg; // avg success
int<lower=0,upper=1> above_avg[J]; // true if theta[j] > mean(theta), meaning the posterior prob. that a given destination is above-average in terms of the elodea introduction rate
int<lower=1,upper=J> rnk[J]; // rank of j
int<lower=0,upper=1> highest[J]; // true if j is highest rank
avg = mean(theta);
for (j in 1:J)
above_avg[j] = (theta[j] > avg);
for (j in 1:J) {
rnk[j] = rank(theta,j) + 1;
highest[j] = rnk[j] == 1;
}
}
```