Double Observer Distance Sampling in STAN

Double observer distance sampliung coded up in STAN

Setup

Load cmdstanr as well as severla packages (data.table and dplyr for data manipulation). Packges like ggplot2, bayesplot and cowplot are used for plotting model results.

Show code

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  cache = FALSE,
  message = FALSE, 
  warning = FALSE
)

library(cmdstanr)
library(data.table)
library(dplyr)
library(ggplot2)
library(bayesplot)
library(cowplot)
library(VicmapR)
library(kableExtra)
options(mc.cores=4)
check_cmdstan_toolchain(fix = TRUE, quiet = TRUE)
register_knitr_engine(override = FALSE)

# WombatSF <- vicmap_query("datavic:CROWNLAND_PLMGEN_STATE_FOREST") %>%
#   filter(LABELSHORT == "Wombat SF") %>%
#   collect() %>%
#   sf::st_combine()

Stan model

Below we construct a mark-recapture distance sampling model for calculating abundance of individuals across 60 sites.

Show code

// double observer distance sampling model
data {
  int<lower=0> N1;                     // number of observations obs2 detect
  int<lower=0> N2;                     // number of observations obs1 detect
  array[N1] int<lower=0,upper=1> M;    // marked and obs2 detect
  array[N2] int<lower=0,upper=1> C;    // captured and obs1 detect
  real delta;                          // bin width
  vector[N1] D1;                       // distance and obs2 detect, not binned
  vector[N2] D2;                       // distance and obs1 detect, not binned
  int<lower=1> n_site;                 // sites
  int<lower=1> n_distance_bins;        // distance bins
  vector[n_distance_bins] midpts;      // midpt of each interval
  array[n_site, n_distance_bins] int y; // observations matrix
  array[n_site] int n_obs;              //number of observations
  int<lower=1> max_int_dist;            // distances used for indexing and predicting
  real<lower=1> max_distance;
  int<lower=0> poisson_ncb;  // number of covariates for poisson model
  matrix[n_site, poisson_ncb] poisson_model_matrix; // poisson model matrix
}

transformed data {
  vector[n_distance_bins] pi; // availability: uniform for line transect

  for (i in 1:n_distance_bins) {
    pi[i] = delta/max_distance; // availability function (line transect)
    //pi[i] = (2 * delta * midpts[i])/max_distance^2 // point transect
  }
}

parameters {
  // mark-recapture parameters
  real alpha1;
  real alpha2;
  real b_distance;
  // abundance parameters
  vector[poisson_ncb] beta_poisson;
  vector[n_site] eps_pois;
  // distance parameters
  real log_sigma; 
}

transformed parameters {
  real<lower=0, upper=1> y01; // observation at distance zero by observer 1 given observer 2 detection
  real<lower=0, upper=1> y02; // observation at distance zero by observer 2 given observer 1 detection
  real<lower=0, upper=1> y0; // observation at distance zero by observer 1 or 2 given other observer detection
  vector[n_distance_bins] p_raw; // detection probability
  vector[n_distance_bins] p_raw_scale; // detection probability for point independence model and accounting for availability
  vector<upper=0>[n_distance_bins] log_p_raw; // log detection probability for point independence model and accounting for availability
  real log_p; // log probability with summed exponential for the multinomial logit model
  real sigma = exp(log_sigma);

  y01 = inv_logit(alpha1); // model prediction for observer 1 at distance = 0
  y02 = inv_logit(alpha2); // model prediction for observer 2 at distance = 0
  y0 = y01 + y02 - y01*y02; // union prediction for observer 1 or observer 2 detecting an individual at distance = 0

    for (i in 1:n_distance_bins) {
      // assuming a half-normal detection fn from line 
      p_raw[i] = exp(-(midpts[i])^2/(2*sigma^2)); // half-normal function (pg 419 of AHM)
      p_raw_scale[i] = y0*p_raw[i]*pi[i]; //  pr(animal occurs and is detected in bin i)
    }

    log_p_raw = log(p_raw_scale);
    log_p = log_sum_exp(log_p_raw);

}

model {
  alpha1 ~ normal(0, 2); // prior for intercept
  alpha2 ~ normal(0, 2); // prior for intercept
  b_distance ~ normal(0, 2); // prior for distance
  beta_poisson ~ normal(0, 2); // prior for poisson model
  eps_pois ~ student_t(4, 0, 1); // prior for overdispersion
  
  M ~ bernoulli_logit(alpha1 + b_distance * D1); // model, conditional on observer 2 detection 
  C ~ bernoulli_logit(alpha2 + b_distance * D2); // model, conditional on observer 1 detection
  
  log_sigma ~ normal(0, 5); // prior for sigma
  
  n_obs ~ poisson_log(poisson_model_matrix * beta_poisson + log_p + eps_pois); 
  
  for (i in 1:n_site) {
  y[i, ] ~ multinomial_logit(log_p_raw);
  }

}

generated quantities {
  array[n_site] int n_unobserved;
  array[n_site] int n;
  // vector<lower=0>[n_site] phi; // occupancy probability
  array[max_int_dist+1] real MR1;
  array[max_int_dist+1] real MR2;
  array[max_int_dist+1] real MRUnion;

  for (i in 1:n_site) { // get site specific counts
    n_unobserved[i] = poisson_rng(exp(poisson_model_matrix[i] * beta_poisson + eps_pois[i]) * (1 - exp(log_p))); 
    n[i] = n_obs[i] + n_unobserved[i];
  }

  for(i in 0:max_int_dist) { // get mark-recapture predictions for distance 0 to max bin distance
    MR1[i+1] = inv_logit(alpha1 + b_distance*i);
    MR2[i+1] = inv_logit(alpha2 + b_distance*i);
    MRUnion[i+1] = MR1[i+1] + MR2[i+1] - (MR1[i+1] * MR2[i+1]); 
  }
}

Based on the above model we simulate double-observer distance sampling data. Basically this involves generating a set of Bernoulli trials for observation processes of a simulated number of true individuals at each site. The success of these observations is based on a half-normal detection function for a line transect (80m maximum half-strip width).

Show code

O1_Effect <- -0.2 # observer 1 variation (not as successful as observer 2)
O2_Effect <- 0.2 # observer 2 variation 
n_site <- 60 # number of sites
det_at_0 <- 0.6 # detection at 0 for 1 observer
y01 <- det_at_0 * (1 + O1_Effect)
y01
y02 <- det_at_0 * (1 + O2_Effect)
y02
y0 <- y01 + y02 - (y01 * y02)
y0
n_distance_bins <- 8 # bins
max_distance <- 80 # you can bias the estimates if this is not set smartly. 
delta<- max_distance/n_distance_bins
transect_length <- 500 # not currently used but may be used as an offset in further models
bin_breakpoints <- seq(0, max_distance, length.out = n_distance_bins + 1)
midpts <- seq(delta/2, max_distance, delta)

sigma <- 30 # sigma detection shape
den = sqrt(2) * sigma
erf <- function(x) 2*pnorm(sqrt(2)*x) - 1
p <- rep(NA, n_distance_bins)
for (b in 1:n_distance_bins) {
  p[b] <- exp(-(midpts[b])^2/(2*sigma^2))
}
p_c <- p 

coords <- data.frame(x = runif(n_site), y = runif(n_site)) 
# coords <- bind_rows(coords, coords)
D <- as.matrix(dist(coords))

# parameter values
mean_lambda <- 15
beta_pois <- c(log(mean_lambda), 2, 1.5)


# prediction coordinates 
# Raster of data 
raster_data <- raster::raster(x = matrix(rev(1:20), nrow = 20) %*% matrix(1:20, nrow = 1), 
                              xmn=0, xmx=1, ymn=0, ymx=1)

names(raster_data) <- 'Elevation'

raster_data[["HBT"]] <- raster::raster(x = matrix(1:20, nrow = 20) %*% matrix(1:20, nrow = 1), 
                              xmn=0, xmx=1, ymn=0, ymx=1)
# normalize
raster::values(raster_data) <- (raster::values(raster_data))^(1/2)

# scale
raster::values(raster_data[[1]]) <- scales::rescale(raster::values(raster_data[[1]]), to = c(-1,1))
raster::values(raster_data[[2]]) <- scales::rescale(raster::values(raster_data[[2]]), to = c(-1,1))

# extract
SiteData <- cbind(raster::extract(raster_data, coords, df = TRUE), coords)


poisson.formula <- as.formula(~ HBT + Elevation)

# Model fitting data 

poisson_model_matrix <- model.matrix(poisson.formula, SiteData)


# response vector
mu <- exp(c(beta_pois[1] + 
              beta_pois[2]*poisson_model_matrix[,2] + 
              beta_pois[3]*poisson_model_matrix[,3] +
              rnorm(n_site, mean = 0, sd = 0.2)))
#mu[zi_calc] <- 0
n_true <- rpois(n_site, mu)

# SiteData$n_true <- n_true
# SiteData %>% 
#   ggplot() + 
#   geom_point(aes(x = x, y = y, colour = n_true))

# Observation Process for Mark Recapture ----------------------------------------------------------

n_obs <- integer()
SFL <- list()
# site error 
# std_er <- rnorm(n = n_site, mean = 0, sd = 0.15)

for(j in 1:n_site) {
  SFL[[j]] <- list()
  
      if(n_true[j] == 0) {
    sample_v <- rep(1, n_distance_bins)
    pr <- 0
    } else {
  samples <- data.frame(db = 1:n_distance_bins) %>%
    left_join(data.frame(idv = rep(1, n_true[j])) %>%
  rowwise() %>%
  mutate(db = sample(x = 1:n_distance_bins, 
                          size = idv, replace = FALSE, 
                          prob = rep(1/n_distance_bins, n_distance_bins))) %>%
  group_by(db) %>%
  summarise(samples = sum(idv)), by = "db") %>%
    as.data.frame()
  
  samples[is.na(samples)] <- 0
  
  sample_v <- samples$samples

    }
  
for(i in 1:n_distance_bins) {
  
  if(sample_v[i] == 0) {
    sample_v[i] <- 1
    pr <- 0
    } else {
      pr <- p_c[i]
    }

   SFL[[j]][[i+1]] <- data.frame(M = rbinom(sample_v[i], size = 1, prob = pr*(det_at_0+O1_Effect)), # observations
                                C = rbinom(sample_v[i], size = 1, prob = pr*(det_at_0+O2_Effect)), # observations
                                D = (bin_breakpoints[i] + bin_breakpoints[i+1])/2, 
                                DReal = sample(x = bin_breakpoints[i]:(bin_breakpoints[i]-1), 
                                               size = sample_v[i], 
                                               replace = TRUE, 
                                               prob = seq(from = p_c[i], to = min(p_c[i+1], p_c[i], na.rm = T), 
                                                          length = length(bin_breakpoints[i]:(bin_breakpoints[i]-1)))), # uniform real distances (pre-binning; will generate better model fit if not binned for MR model)
                                Site = j)
}
  o_sf <- bind_rows(SFL[[j]])
  n_obs[j] <- nrow(o_sf[o_sf$M == 1 | o_sf$C == 1, ])
}

SF <- bind_rows(SFL)
SF <- SF %>% 
  rowwise() %>%
  mutate(Union = max(M, C)) %>%
  ungroup() %>%
  as.data.frame()

y <- as.data.table(SF) %>%
  dcast(Site ~ D, value.var = "Union", fun.aggregate = sum) %>%
  tibble::column_to_rownames("Site") %>%
  as.matrix()

stan_d_mr <- list(N1 = sum(SF$C), 
                  N2 = sum(SF$M), 
                  M = SF[SF$C == 1, "M"], 
                  C = SF[SF$M == 1, "C"], 
                  delta = delta, 
                  transect_length = transect_length,
                  D1 = SF[SF$C == 1, "DReal"], 
                  D2 = SF[SF$M == 1, "DReal"], 
                  n_site = n_site,
                  n_distance_bins = n_distance_bins,
                  midpts = midpts,
                  y = y,
                  n_obs = n_obs, 
                  max_int_dist = ceiling(max_distance), 
                  poisson_ncb = ncol(poisson_model_matrix), 
                  poisson_model_matrix = poisson_model_matrix,
                  max_distance = max_distance)

# MCMC settings
ni <- 2000
nw <- 1000
nt <- 1
nb <- 1000
nc <- 4

inits = lapply(1:nc, function(i) list(log_sigma=runif(1,1,5),
                                      alpha1=runif(1),
                                      alpha2=runif(1),
                                      b_distance=0, 
                                      gp_sigma = runif(1), 
                                      gp_phi = runif(1, min = 0, max = 5), 
                                      gp_eta = runif(1), 
                                      gp_z = runif(1),
                                      # beta_hurdle = runif(1), 
                                      beta_poisson = runif(1)))

fit <- mod$sample(data = stan_d_mr, 
                  parallel_chains = nc, 
                  show_messages = FALSE, 
                  save_warmup = FALSE,
                  iter_sampling = ni, 
                  iter_warmup = nw, 
                  init = inits)

Model diagnostics

Inspect traceplots and the distribution of two key model parameters (log_lambda and log_sigma).

Show code

# Model Diagnostics
diag1 <- bayesplot::mcmc_trace(fit$draws(c("lp__", "log_sigma", "log_p", "b_distance", "alpha1", "alpha2", "beta_poisson")))
diag2 <- bayesplot::mcmc_pairs(fit$draws(c("log_sigma", "log_p")))
plot_grid(diag1, diag2, ncol = 1)

Figure 1: Model diagnostics for the mark-recapture distance-sampling model.

Estimate population size/s

For each individual site (as well as the whole population) you can estimate the population size. This is done in the generated quantities block within the STAN model. From this we can draw the posterior distribution and compare it to the true values.

Show code

# Model predictions 
true_vals_df <- data.frame(variable = paste0("n[", 1:n_site, "]"), true = n_true)

modelled_n <- fit$draws('n', format = "draws_matrix", inc_warmup = FALSE) %>%
  as.data.frame() %>%
  tibble::rowid_to_column(var = "draw") %>%
  as.data.table() %>%
  melt(id.var = "draw") %>%
  # group_by(variable) %>%
  # arrange(desc(value)) %>%
  # slice_sample(n = 6000, weight_by = c(rep(0, 1000), rep(1, 6000), rep(0, 1000))) %>%
  group_by(draw) %>%
  summarise(PopSize = sum(value)) 

modelled_n_zi <- fit$draws('n', format = "draws_matrix", inc_warmup = FALSE) %>%
  as.data.frame() %>%
  tibble::rowid_to_column(var = "draw") %>%
  as.data.table() %>%
  melt(id.var = "draw") %>%
  group_by(variable) %>%
  summarise(PopSizeSite = mean(value),
            PopSizeMedian = median(value),
            PopSizeLCI = quantile(value, 0.10),
            PopSizeUCI = quantile(value, 0.90)) %>%
  ungroup()

modelled_n_unobserved <- fit$draws('n_unobserved', format = "draws_matrix", inc_warmup = FALSE) %>%
  as.data.frame() %>%
  tibble::rowid_to_column(var = "draw") %>%
  as.data.table() %>%
  melt(id.var = "draw") %>%
  # group_by(variable) %>%
  # arrange(desc(value)) %>%
  # slice_sample(n = 6000, weight_by = c(rep(0, 1000), rep(1, 6000), rep(0, 1000))) %>%
  group_by(draw) %>%
  summarise(PopSize = sum(value))

pop1 <- modelled_n %>%
  ggplot() +
  geom_histogram(aes(x = PopSize), colour = "DarkBlue", fill = "LightBlue", lwd = 0.2, bins = 40) +
  geom_vline(data = true_vals_df, aes(xintercept = sum(true)), lwd = 1, linetype = "dashed", colour = "DarkRed") +
  # geom_vline(data = modelled_n_zi, aes(xintercept = sum(PopSizeMedian)), lwd = 1, linetype = "solid", colour = "DarkGreen") +
  # geom_vline(data = modelled_n_zi, aes(xintercept = sum(PopSizeUCI)), lwd = 1, linetype = "solid", colour = "DarkGreen") +
  scale_y_discrete(expand = c(0, 0)) +
  ylab("Density") +
  xlab("Total Population") +
  theme_classic() 
  # ggtitle("Trimmed mean estimates of total population \nsize using the middle 75 % of estimates")

pop2 <- modelled_n %>%
  ggplot() +
  geom_histogram(aes(x = PopSize), colour = "DarkBlue", fill = "LightBlue", lwd = 0.2, bins = 40) +
  geom_histogram(aes(x = PopSize), colour = "Tomato", fill = "DarkRed", lwd = 0.2, bins = 40, data = modelled_n_unobserved) +
  geom_vline(data = true_vals_df, aes(xintercept = sum(true)), lwd = 1, linetype = "dashed", colour = "DarkRed") +
  # geom_vline(data = modelled_n_zi, aes(xintercept = sum(PopSizeMedian)), lwd = 1, linetype = "solid", colour = "DarkGreen") +
  scale_y_discrete(expand = c(0, 0)) +
  ylab("Density") +
  xlab("Population: Unobserved and Total") +
  theme_classic() 
  # ggtitle("Trimmed mean estimates of total population\n size and the total unobserved population \nusing the middle 75 % of estimates")

pop3 <- bayesplot::mcmc_intervals(fit$draws("n")) +
  geom_point(data = true_vals_df, aes(x = true, y = variable), fill = "Red", shape = 21, size = 2, position =  position_nudge(y = 0)) +
  geom_point(data = modelled_n_zi, aes(x = PopSizeSite, y = variable), fill = "DarkGreen", shape = 25, size = 2, position =  position_nudge(y = 0.5)) +
  scale_y_discrete(expand = c(0.01, 0)) 

plot_grid(plot_grid(pop1, pop2, ncol = 1), pop3, ncol = 2, labels = "AUTO")

Figure 2: (A) Overall population size estimates vs the true value (red line) as well as estimates for each sites (B) vs the true simulated value (Red dot). Plot B shows posterior distribution of individuals at each site, with the true value presented as a red dot. The median estimate for each site is presented as the blue dot with MCMC intervals (blue lines). Finally, the green triangle is the mean estimate for the site.

View the predicted detection function

The model used here assumed a point-independence arrangement. Thus the detection function is multiplied by the probability of detection at distance zero. Given that we are conducting double-observer distance sampling the detection function shown below has been modified to be for the union of detections of observer 1 and 2. That is; at least one observer detects the animal.

Show code

det_prob_pred <- fit$summary(c("MRUnion", "MR1", "MR2"))
det_prob_pred$Distance <- rep(seq(from = 0, to = max_distance, by = 1), 3)
det_prob_pred$Observer <- rep(c("Obs1 U Obs2", "Obs1", "Obs2"), each = max_distance+1)

det_prob_MR <- det_prob_pred %>%
  ggplot(aes(x = Distance, y = mean, fill = Observer, colour = Observer)) +
  geom_hline(yintercept = y0, colour = "#395D9CFF", linetype = "dashed", lwd = 1.5, alpha = 0.75) +
  geom_hline(yintercept = y02, colour = "#60CEACFF", linetype = "dashed", lwd = 1.5, alpha = 0.75) +
  geom_hline(yintercept = y01, colour = "#0B0405FF", linetype = "dashed", lwd = 1.5, alpha = 0.75) +
  geom_ribbon(aes(ymin = q5, ymax = q95), alpha = 0.5) +
  geom_line() +
  viridis::scale_fill_viridis(discrete = T, option = "G", end = 0.8) +
  viridis::scale_colour_viridis(discrete = T, option = "G", end = 0.8) +
  ylim(c(0,1)) +
  scale_x_continuous(expand = c(0, 0)) +
  ggtitle("Mark Recapture Model") +
  theme_classic() +
  theme(legend.position = "top")

y_combined <- colSums(y) %>%
  as.data.frame() %>%
  `colnames<-`("Count") %>%
  tibble::rownames_to_column("Distance") %>%
  mutate(Distance = as.integer(Distance), 
         Prop = Count/(max(Count)))

hnormsim <- rep(NA, 99)
for (b in 1:99) {
  hnormsim[b] <- exp(-(b+1)*b/(2*sigma^2))
}

y0_summary <- fit$summary("y0")


hnormsim_scaled <- data.frame(HNS = hnormsim + hnormsim - hnormsim^2, # union for half-normal
                              MRDS = y0_summary$mean*hnormsim,
                              MRDS_LCI = y0_summary$q5*hnormsim,
                              MRDS_UCI = y0_summary$q95*hnormsim, 
                              Distance = 2:100)


detection_plot_HN <- ggplot(aes(x = Distance), data = hnormsim_scaled) +
  geom_col(aes(y = Prop), fill = "grey70", colour = "grey20", width = 10, data = y_combined, alpha = 0.7) +
  geom_line(aes(y = HNS)) +
  geom_ribbon(aes(ymin = MRDS_LCI, ymax = MRDS_UCI), alpha = 0.5, fill = "Red") +
  geom_line(aes(y = MRDS), colour = "Red") +
  ggtitle("Mark Recapture Distance Sampling Model", subtitle = "Red = g(0) x p(dist)") +
  theme_classic()

plot_grid(det_prob_MR, detection_plot_HN, ncol = 2)

Figure 3: Comparison of logistic regression for the mark-recapture model and the distance-sampling model. The histogram shows the proportional detection rates from the simulated data. The red line shows the model predictions for the detection function; accounting for g(0), or the detection probability being less than 1 at distance 0

Model Results

Comparison of the model results against the true values used during simulation.

Show code

ParamsofInterest <- list()
ParamsofInterest[["MRParams"]] <- c("y0", "y01", "y02", "b_distance")
ParamsofInterest[["DSParams "]]<- c("log_sigma", "p_raw")
ParamsofInterest[["MRDSParams"]] <- c("log_p", "log_p_raw")
# ParamsofInterest[["hurdleParams"]] <- c("beta_hurdle")
ParamsofInterest[["poissonParams"]] <- c("beta_poisson")
# ParamsofInterest[["GPParams"]] <- c("gp_w")
Summary_Table <- fit$summary(unlist(ParamsofInterest))
Summary_Table$Group <- c(rep("Mark Recapture", 4), 
                         rep("Distance Sampling", 9), 
                         rep("Mark Recapture Distance Sampling", 9), 
                         # rep("Hurdle Model", ncol(hurdle_model_matrix)), 
                         rep("Poisson Model", ncol(poisson_model_matrix))) 
                         # rep("Gaussian Process", n_site))

Summary_Table$TrueSim <- c(y0, y01, y02, NA, 
                           log(sigma), p, 
                           log(sum(exp(log(p*y0*delta/max_distance)))), 
                           log(p*y0*delta/max_distance), 
                           # beta_bern, 
                           beta_pois)

Summary_Table <- Summary_Table %>%
  mutate(InRange = case_when(TrueSim > q5 & TrueSim < q95 ~ TRUE, 
                             TRUE ~ FALSE))
                        
Summary_Table %>%
  dplyr::select(Group, variable, TrueSim, mean, sd, q5, q95, rhat) %>%
  kbl("html", digits = 3, caption = "Comparison of estimated parameters against the true value used during simulation") %>%
  kable_styling(bootstrap_options = "striped") %>%
  collapse_rows(1:2) %>%
    row_spec(which(Summary_Table$InRange), bold = TRUE, background = "DarkGreen", color = "white") 

Table 1: Comparison of estimated parameters against the true value used during simulation

Group	variable	TrueSim	mean	sd	q5	q95	rhat
Mark Recapture	y0	0.854	0.856	0.024	0.815	0.893	1.001
	y01	0.480	0.452	0.034	0.396	0.509	1.000
	y02	0.720	0.737	0.035	0.678	0.794	1.001
	b_distance	NA	-0.049	0.007	-0.059	-0.038	1.001
Distance Sampling	log_sigma	3.401	3.462	0.039	3.400	3.529	1.000
	p_raw[1]	0.986	0.988	0.001	0.986	0.989	1.000
	p_raw[2]	0.882	0.895	0.008	0.882	0.908	1.000
	p_raw[3]	0.707	0.735	0.018	0.706	0.764	1.000
	p_raw[4]	0.506	0.547	0.026	0.505	0.590	1.000
	p_raw[5]	0.325	0.369	0.029	0.323	0.418	1.000
	p_raw[6]	0.186	0.226	0.026	0.185	0.272	1.000
	p_raw[7]	0.096	0.126	0.021	0.095	0.162	1.000
	p_raw[8]	0.044	0.064	0.014	0.043	0.089	1.000
Mark Recapture Distance Sampling	log_p	-0.920	-0.863	0.045	-0.936	-0.787	1.000
	log_p_raw[1]	-2.251	-2.248	0.028	-2.297	-2.205	1.001
	log_p_raw[2]	-2.362	-2.347	0.029	-2.397	-2.301	1.000
	log_p_raw[3]	-2.584	-2.544	0.037	-2.606	-2.485	1.000
	log_p_raw[4]	-2.917	-2.840	0.055	-2.930	-2.750	1.000
	log_p_raw[5]	-3.362	-3.235	0.083	-3.372	-3.098	1.000
	log_p_raw[6]	-3.917	-3.728	0.120	-3.925	-3.531	1.000
	log_p_raw[7]	-4.584	-4.320	0.165	-4.593	-4.048	1.000
	log_p_raw[8]	-5.362	-5.010	0.219	-5.372	-4.654	1.000
Poisson Model	beta_poisson[1]	2.708	2.694	0.162	2.428	2.963	1.005
	beta_poisson[2]	2.000	2.197	0.359	1.615	2.786	1.000
	beta_poisson[3]	1.500	1.565	0.351	0.992	2.151	1.002