32 Chapter 19 Exercises - Reproducing the the ‘serial dilution assay’

32.2.1 Model specification

The model uses a normal likelihood to describe the posterior distribution \(p(y|x)\). The mean of the likelihood is defined for a given concentration \(x\) using the standard equation used in the field:

\[ \text{E}[y | x, \beta] = g(x, \beta) = \beta_1 + \frac{\beta_2}{1 + (x/\beta_3)^{-\beta_4}} \\ \]

The model is a scaled and shifted logistic curve. This structure results in the following interpretations for \(\beta\), all of which are restricted to positive values:

1. \(\beta_1\): color intensity when the concentration is 0
2. \(\beta_2\): increase to saturation
3. \(\beta_3\): the inflection point of the curve
4. \(\beta_4\): rate of saturation

Below are the prior distributions for \(\beta\). Note that they are are drastically different scales - this is critical to help the model fit the data.

\[ \beta_1 \sim N(10, 2.5) \\ \beta_2 \sim N(100, 5) \\ \beta_3 \sim N(0, 1) \\ \beta_4 \sim N(0, 2.5) \]

The measurement error of the model, representing the variance in the model’s likelihood is defined as follows:

\[ \tau(\alpha, \sigma_y, g(x, \beta), A) = \lgroup \frac{g(x,\beta)}{A} \rgroup^{2\alpha} \sigma^2_y \]

Here, \(\alpha\), restricted to lie between 0 and 1, allows the variance to be higher for larger measurement values. \(A\) is a constant (set to 30 by the authors) that allows \(\sigma_y\) to be more easily interpreted as the variance from “typical” measurements. Below are the priors for the new variables in the model.

\[ \alpha \sim \text{Beta}(1, 1) \qquad \sigma \sim |N(0, 2.5)| \]

32.2.2 In Stan

Below is the Stan code for the model. It looks very similar to the mathematical description of the model, a nice feature of the Stan probabilistic programming language.

The centrality and variance of the likelihood are calculated separately as g and tau so they can be used in the model and generated quantities block without duplicating the code. The log_lik is calculated so that PSIS-LOO cross validation can be estimated. I also included the ability to provide new data to make predictions over as xnew.

dilution_model_file <- here::here("models", "serial-dilution.stan")
writeLines(readLines(dilution_model_file))

#> data {
#>   int<lower=0> N;           // number of data points
#>   int<lower=0> A;           // constant used in model of measurement error
#>   vector<lower=0>[N] x;     // concentration values
#>   vector<lower=0>[N] y;     // observed color intensity
#>   int<lower=0> M;           // number of new x values
#>   vector<lower=0>[M] xnew;  // new x values
#> }
#> 
#> parameters {
#>   vector<lower=0>[4] beta;
#>   real<lower=0,upper=1> alpha;
#>   real<lower=0> sigma;
#> }
#> 
#> transformed parameters {
#>   vector<lower=0>[N] g;
#>   vector<lower=0>[N] tau;
#> 
#>   for (i in 1:N) {
#>     g[i] = beta[1] + beta[2] / (1 + (x[i] / beta[3]) ^ (-beta[4]));
#>     tau[i] = ((g[i] / A) ^ (2.0 * alpha)) * (sigma ^ 2.0);
#>   }
#> }
#> 
#> model {
#>   // Priors
#>   alpha ~ beta(1, 1);
#>   beta[1] ~ normal(10, 2.5);
#>   beta[2] ~ normal(100, 5);
#>   beta[3] ~ normal(0, 1);
#>   beta[4] ~ normal(0, 2.5);
#>   sigma ~ normal(0, 2.5);
#> 
#>   // Likelihood
#>   for (i in 1:N) {
#>     y[i] ~ normal(g[i], tau[i]);
#>   }
#> }
#> 
#> generated quantities {
#>   vector[N] ypred;
#>   vector[N] log_lik;
#> 
#>   vector[M] g_hat;
#>   vector[M] tau_hat;
#>   vector[M] ynew;
#> 
#>   for (i in 1:N) {
#>     ypred[i] = normal_rng(g[i], tau[i]);
#>     log_lik[i] = normal_lpdf(y[i] | g[i], tau[i]);
#>   }
#> 
#>   for (i in 1:M) {
#>     g_hat[i] = beta[1] + beta[2] / (1 + (xnew[i] / beta[3]) ^ (-beta[4]));
#>     tau_hat[i] = ((g_hat[i] / A) ^ (2.0 * alpha)) * (sigma ^ 2.0);
#>     ynew[i] = normal_rng(g_hat[i], tau_hat[i]);
#>   }
#> }

32.2.3 Sampling

As mentioned above, specifically defining the prior distributions for each \(\beta\) is necessary for MCMC to accurately sample from the posterior. With those helping restrict the range of their values, the model fit very well.

xnew <- seq(0, max(dilution_standards_data$conc), 0.001)
model_data <- list(
  N = nrow(dilution_standards_data),
  A = 30,
  x = dilution_standards_data$conc,
  y = dilution_standards_data$y,
  M = length(xnew),
  xnew = xnew
)

dilution_model <- stan(
  dilution_model_file,
  model_name = "serial-dilution",
  data = model_data,
  refresh = 1000
)

32.2.4 Posterior distributions

The next step is to analyze the posterior draws of the model. We can check the success of MCMC by visualizing the traces of the chains, looking for good mixing (“fuzzy caterpillars”) and checking diagnostic values such as \(\widehat{R}\) and \(n_\text{eff}\). The trace plots are shown below followed by a table of the posteriors with the diagnostic values. Everything looks good suggesting MCMC was successful.

model_pars <- c("beta", "alpha", "sigma")
rstan::stan_trace(dilution_model, pars = model_pars, ncol = 2, alpha = 0.7) +
  scale_x_continuous(expand = expansion(c(0, 0))) +
  scale_y_continuous(expand = expansion(c(0.02, 0.02))) +
  theme(legend.position = "bottom")

print(dilution_model, pars = model_pars)

#> Inference for Stan model: serial-dilution.
#> 4 chains, each with iter=2000; warmup=1000; thin=1; 
#> post-warmup draws per chain=1000, total post-warmup draws=4000.
#> 
#>           mean se_mean   sd  2.5%   25%    50%    75%  97.5% n_eff Rhat
#> beta[1]  14.33    0.02 0.55 13.12 14.08  14.36  14.61  15.33  1074    1
#> beta[2] 102.70    0.10 4.30 94.48 99.67 102.61 105.60 111.32  1755    1
#> beta[3]   0.06    0.00 0.01  0.05  0.06   0.06   0.07   0.08  1633    1
#> beta[4]   1.13    0.00 0.08  0.97  1.07   1.12   1.18   1.29  1870    1
#> alpha     0.72    0.01 0.17  0.33  0.61   0.74   0.85   0.98  1074    1
#> sigma     1.33    0.01 0.24  0.98  1.16   1.30   1.45   1.92   893    1
#> 
#> Samples were drawn using NUTS(diag_e) at Wed Feb  9 06:44:24 2022.
#> For each parameter, n_eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor on split chains (at 
#> convergence, Rhat=1).

The following density plots show the posterior distributions of the model parameters \(\beta\), \(\alpha\), and \(\sigma\).

rstan::stan_dens(
  dilution_model,
  pars = model_pars,
  separate_chains = FALSE,
  alpha = 0.6
) +
  scale_x_continuous(expand = expansion(c(0, 0))) +
  scale_y_continuous(expand = expansion(c(0, 0.02)))

32.2.5 Posterior predictive check

Below is a plot of the posterior predictive distributions of the model on the original data. 1,000 individual simulations are plotted in blue and the mean in black. The simulated curves visually appear to correspond well with the observed data indicating the model has good fit.

dilution_post_pred <- rstan::extract(dilution_model, "ypred")$ypred %>%
  as.data.frame() %>%
  as_tibble() %>%
  set_names(seq(1, ncol(.))) %>%
  mutate(draw = 1:n()) %>%
  pivot_longer(-c(draw), names_to = "idx") %>%
  left_join(
    dilution_standards_data %>% mutate(idx = as.character(1:n())),
    by = "idx"
  )

plt_n_draws <- 1000
plt_draws <- sample(1:max(dilution_post_pred$draw), plt_n_draws)

ppc_mean <- dilution_post_pred %>%
  group_by(conc) %>%
  summarize(value = mean(value)) %>%
  ungroup()

dilution_post_pred %>%
  filter(draw %in% !!plt_draws) %>%
  mutate(grp = glue::glue("{draw}-{rep}")) %>%
  ggplot(aes(x = conc, y = value)) +
  geom_line(aes(group = grp), alpha = 0.05, color = SNS_BLUE) +
  geom_line(group = "a", data = ppc_mean, color = "black") +
  geom_point(data = ppc_mean, color = "black") +
  geom_line(
    aes(y = y, group = rep),
    data = dilution_standards_data,
    color = STAN_RED
  ) +
  geom_point(aes(y = y), data = dilution_standards_data, color = STAN_RED) +
  scale_x_continuous(expand = expansion(c(0, 0))) +
  scale_y_continuous(expand = expansion(c(0.02, 0.02))) +
  labs(
    x = "concentration",
    y = "y",
    title = "Posterior predictive distribution"
  )

I also had the model make posterior predictions on concentrations across the observed range at smaller step-sizes. The mean and 89% HDI are shown in blue below along with the observed data in red. The inset plot is a zoomed-in view of the posterior predictive distribution at the lower concentrations.

ynew_mean <- apply(rstan::extract(dilution_model, pars = "ynew")$ynew, 2, mean)
ynew_hdi <- apply(
  rstan::extract(dilution_model, pars = "ynew")$ynew,
  2,
  bayestestR::hdi,
  ci = 0.89
)
ynew_ppc <- tibble(
  conc = xnew,
  ynew_mean = ynew_mean,
  ynew_hdi_low = purrr::map_dbl(ynew_hdi, ~ unlist(.x)[[2]]),
  ynew_hdi_hi = purrr::map_dbl(ynew_hdi, ~ unlist(.x)[[3]])
)

plot_posterior_pred <- function(ppc_df, obs_df, pt_size = 1.5) {
  ppc_df %>%
    ggplot(aes(x = conc, y = ynew_mean)) +
    geom_ribbon(
      aes(ymin = ynew_hdi_low, ymax = ynew_hdi_hi),
      fill = SNS_BLUE,
      alpha = 0.5
    ) +
    geom_line(group = "a") +
    geom_line(
      aes(y = y, group = rep),
      data = obs_df,
      color = STAN_RED
    ) +
    geom_point(aes(y = y), data = obs_df, size = pt_size, color = STAN_RED) +
    scale_x_continuous(expand = expansion(c(0, 0))) +
    scale_y_continuous(expand = expansion(c(0.02, 0.02)))
}

ppc_plot <- plot_posterior_pred(ynew_ppc, dilution_standards_data) +
  labs(
    x = "concentration",
    y = "y",
    title = "Posterior predictive distribution"
  )

sub_max <- 0.04
sub_ppc_plot <- plot_posterior_pred(
  ynew_ppc %>% filter(conc <= sub_max),
  dilution_standards_data %>% filter(conc <= sub_max),
  pt_size = 0.6
) +
  theme(axis.title = element_blank())

ppc_plot +
  inset_element(sub_ppc_plot, left = 0.5, bottom = 0.05, right = 0.9, top = 0.5)