In maomlab/BayesPharma: Tools for Bayesian Analysis of Non-Linear Pharmacology Models
Hill Equation {#sec-hill}
In this case study, we are going to reanalyze the dose responses of 4
Kappa Opioid receptor (KOR) antagonists from a study performed by Margolis et
al. (-@Margolis2020-bm) using the BayesPharma
package. Whole cell electrophysiology in acute rat
midbrain slices was used to evaluate the pharmacological properties of
four novel KOR antagonists: BTRX-335140, BTRX-395750, PF-04455242, and
JNJ-67953964.
Originally, the dose-response analysis was performed using
the drc package in R, which implements the minimization of negative
log likelihood function and reduces to least square estimation for a
continuous response. The data were normalized to % baseline agonist response,
measuring the extent of blockade of agonist response in the presence of
different concentrations of antagonist. These data were then fit to a 4-parameter
log-logistic dose response model, setting the top (max agonist response) to
100% and estimating the IC50, its variance, and the bottom (min agonist
response).
Fitting the sigmoid model
Using the BayesPharma package, we can re-fit the sigmoid model with a
negative slope, and fix the top parameter to 100 as the response is
normalized to a no-antagonist baseline.
For the prior, we are going to use a normal distribution because the
response values are continuous. First, we will run the analysis with the
top (max response) parameter prior set to a constant value of 100
because top is normalized to 100 and is the default broad prior for the ic50,
hill, and bottom parameters. Broad priors represent unbiased uncertainty and
provide an opportunity for extreme responses.
The level of informativeness of the prior will affect how much
influence the prior has on the model. Here is more information on
prior choice
recommendations.
\scriptsize
kor_prior <- BayesPharma::sigmoid_antagonist_prior(top = 100)
kor_prior
## prior class coef group resp dpar nlpar lb ub source
## normal(-6, 2.5) b ic50 <NA> <NA> user
## normal(-1, 1) b hill <NA> 0.01 user
## constant(100) b top <NA> <NA> user
## normal(0, 0.5) b bottom <NA> <NA> user
\normalsize
Prior predictive checks
Following the Bayesian workflow, before fitting the model, it is good to check
the prior predictive distributions to see if they are compatible with the domain
expertise. So, before running the model, we will verify that the prior
distributions cover a plausible range of values for each parameter. To do this,
we want to sample only from the prior distributions by adding
sample_prior = "only" as an argument to the sigmoid_model
function. We will use the default response distribution of the model
(family = gaussian()).
\scriptsize
kor_sample_prior <- BayesPharma::sigmoid_model(
data = kor_antag |> dplyr::select(substance_id, log_dose, response),
formula = BayesPharma::sigmoid_antagonist_formula(),
prior = kor_prior,
init = BayesPharma::sigmoid_antagonist_init(),
sample_prior = "only")
\normalsize
And then plot of the prior predictive distributions:
\scriptsize
kor_sample_prior |>
BayesPharma::plot_density_distribution()
KOR antagonists prior distribution
\normalsize
To sample from the model, we will use the Stan NUTs Hamiltonian Monte
Carlo, and initialize the parameters to the prior means to help with
model convergence, using the default values of ec50 = -9, hill = -1,
top = 100, bottom = 0.
\scriptsize
kor_model <- BayesPharma::sigmoid_model(
data = kor_antag |> dplyr::select(substance_id, log_dose, response),
formula = BayesPharma::sigmoid_antagonist_formula(
predictors = 0 + substance_id),
prior = kor_prior,
init = BayesPharma::sigmoid_antagonist_init())
\normalsize
Analyzing model fit
The brms generated model summary shows the formula that the expected
response a is sigmoid function of the log_dose with four parameters, and
a shared Gaussian distribution. Each parameter is dependent on the
substance_id. Since we want to fit a separate model for each substance, we
include a 0 + to indicate that there is no common intercept. The data
consists of 73 data points and the posterior sampling was done in 4
chains each with 8000 steps with 4000 steps of warm-up. The population
effects for each parameter summarize the marginal posterior distributions,
as well as the effective sample size in the bulk and tail. This gives
an indication of the sampling quality, with an ESS of > 500 samples
being good for this type of model.
\scriptsize
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: response ~ sigmoid(ic50, hill, top, bottom, log_dose)
## ic50 ~ 0 + substance_id
## hill ~ 0 + substance_id
## top ~ 0 + substance_id
## bottom ~ 0 + substance_id
## Data: data (Number of observations: 73)
## Draws: 4 chains, each with iter = 8000; warmup = 4000; thin = 1;
## total post-warmup draws = 16000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## ic50_substance_idBTRX_335140 -8.84 0.20 -9.20 -8.42 1.00 15027 8248
## ic50_substance_idBTRX_395750 -8.24 0.41 -8.93 -7.36 1.00 11072 5925
## ic50_substance_idJNJ -9.15 0.31 -9.77 -8.50 1.00 15946 10770
## ic50_substance_idPF -6.14 1.08 -7.69 -3.32 1.00 7646 5332
## hill_substance_idBTRX_335140 -1.47 0.59 -2.85 -0.57 1.00 14309 9824
## hill_substance_idBTRX_395750 -0.91 0.52 -2.28 -0.26 1.00 10745 7238
## hill_substance_idJNJ -1.01 0.51 -2.38 -0.41 1.00 14567 11520
## hill_substance_idPF -0.31 0.24 -0.92 -0.03 1.00 7908 5627
## bottom_substance_idBTRX_335140 -0.00 0.49 -0.96 0.95 1.00 17746 11293
## bottom_substance_idBTRX_395750 0.01 0.50 -0.97 0.99 1.00 18472 10902
## bottom_substance_idJNJ -0.01 0.49 -0.99 0.96 1.00 17703 11821
## bottom_substance_idPF 0.01 0.50 -0.97 0.99 1.00 16850 11426
## top_substance_idBTRX_335140 100.00 0.00 100.00 100.00 NA NA NA
## top_substance_idBTRX_395750 100.00 0.00 100.00 100.00 NA NA NA
## top_substance_idJNJ 100.00 0.00 100.00 100.00 NA NA NA
## top_substance_idPF 100.00 0.00 100.00 100.00 NA NA NA
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 32.17 2.90 27.13 38.38 1.00 14788 11389
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
\normalsize
Traceplot: The model ran without warning messages, meaning there were no parameter
value problems or MCMC conflicts. The bulk and tail ESS indicate high
resolution and stability. The R-hat for each parameter equals 1.00 and
the traceplot shows the chains mixed well, indicating the chains
converged.
\scriptsize
kor_model |>
bayesplot::mcmc_trace()
plot of chunk kor-model-traceplot
\normalsize
Compare prior and posterior marginal distributions: Displayed below is a
plot for the prior and posterior distributions of the parameters (prior is pink
and posterior is teal). This can be useful for comparing the density
distribution of the prior and posterior produced by the model:
\scriptsize
BayesPharma::plot_prior_posterior_densities(
model = kor_model,
title_label = "")
KOR antagonists model, compare prior and posterior distributions for each substance
\normalsize
Displayed below is a plot of the posterior distributions for each
parameter with the confidence intervals and mean. This is a useful
visual of the model results and can highlight the mode and high-density
intervals:
\scriptsize
BayesPharma::plot_posterior_density(
kor_model,
title_label = "")
## Error in `dplyr::filter()`:
## ℹ In argument: `!...`.
## Caused by error in `.data[[".variable"]]`:
## ! Column `.variable` not found in `.data`.
\normalsize
Displayed below is a plot of a sample of 100 sigmoid dose-response
curves from the posterior distribution (purple) and the median
quantile intervals:
\scriptsize
BayesPharma::plot_posterior_draws(
model = kor_model,
title = "")
KOR antagonists, posterior draws
\normalsize
Comparing alternative models
To test the sensitivity of the analysis to the prior, we can re-fit the model
with a more informative prior:
\scriptsize
## prior class coef group resp dpar nlpar lb ub source
## normal(-8.5, 0.5) b ic50 <NA> <NA> user
## normal(-1, 0.5) b hill <NA> 0.01 user
## constant(100) b top <NA> <NA> user
## normal(10, 15) b bottom <NA> <NA> user
\normalsize
Re-fitting the model with the kor_prior2 gives:
\scriptsize
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: response ~ sigmoid(ic50, hill, top, bottom, log_dose)
## ic50 ~ 0 + substance_id
## hill ~ 0 + substance_id
## top ~ 0 + substance_id
## bottom ~ 0 + substance_id
## Data: data (Number of observations: 73)
## Draws: 4 chains, each with iter = 8000; warmup = 4000; thin = 1;
## total post-warmup draws = 16000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat
## ic50_substance_idBTRX_335140 -8.82 0.21 -9.23 -8.37 1.00
## ic50_substance_idBTRX_395750 -8.58 0.31 -9.16 -7.94 1.00
## ic50_substance_idJNJ -8.95 0.31 -9.53 -8.32 1.00
## ic50_substance_idPF -8.18 0.45 -9.09 -7.30 1.00
## hill_substance_idBTRX_335140 -1.19 0.37 -1.99 -0.57 1.00
## hill_substance_idBTRX_395750 -1.06 0.39 -1.90 -0.42 1.00
## hill_substance_idJNJ -0.85 0.33 -1.65 -0.39 1.00
## hill_substance_idPF -0.72 0.45 -1.72 -0.07 1.00
## bottom_substance_idBTRX_335140 1.71 10.47 -19.21 21.94 1.00
## bottom_substance_idBTRX_395750 15.44 11.10 -7.30 36.02 1.00
## bottom_substance_idJNJ -2.77 9.82 -23.22 15.35 1.00
## bottom_substance_idPF 31.22 11.59 7.03 52.03 1.00
## top_substance_idBTRX_335140 100.00 0.00 100.00 100.00 NA
## top_substance_idBTRX_395750 100.00 0.00 100.00 100.00 NA
## top_substance_idJNJ 100.00 0.00 100.00 100.00 NA
## top_substance_idPF 100.00 0.00 100.00 100.00 NA
## Bulk_ESS Tail_ESS
## ic50_substance_idBTRX_335140 13087 10643
## ic50_substance_idBTRX_395750 12995 10960
## ic50_substance_idJNJ 13299 11776
## ic50_substance_idPF 13156 10947
## hill_substance_idBTRX_335140 14798 9835
## hill_substance_idBTRX_395750 13437 9765
## hill_substance_idJNJ 13283 11670
## hill_substance_idPF 8776 5539
## bottom_substance_idBTRX_335140 14068 12061
## bottom_substance_idBTRX_395750 12671 10835
## bottom_substance_idJNJ 12510 10615
## bottom_substance_idPF 10495 10156
## top_substance_idBTRX_335140 NA NA
## top_substance_idBTRX_395750 NA NA
## top_substance_idJNJ NA NA
## top_substance_idPF NA NA
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 31.86 2.83 26.91 37.99 1.00 14470 11441
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
\normalsize
Comparing the Two Models Using LOO-Comparison: One way to evaluate the
quality of a model is for each data-point, re-fit the model with remaining
points and evaluate the log probability of the point in the posterior
distribution. Taking the expectation across all points gives the Expected Log
Pointwise predictive Density (ELPD). Since this is computationally challenging
to re-fit the model for each point, if the model fits the data reasonably well,
then the ELPD can be approximated using the Pareto smoothed importance sampling
(PSIS). Using the LOO package, Pareto k value for each data point is computed
that indicate how well that data point is fit by the model, where k < 0.5 is
good, 0.5 <= k < 0.7 is OK, and 0.7 <= k is bad. Evaluating the model for
the KOR antagonists shows that the model fits the data well.
\scriptsize
##
## Computed from 16000 by 73 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -360.5 6.5
## p_loo 6.9 1.1
## looic 721.0 13.1
## ------
## MCSE of elpd_loo is 0.0.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 1.4]).
##
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
\normalsize
Since ELPD is a global measure of model fit, it can be used to compare models.
Using loo_compare from the LOO package returns the elpd_diff and se_diff
for each model relative the model with the lowest ELPD. The kor_model2, the
model with more informative prior, is the preferred model, but not
significantly.
\scriptsize
## No problematic observations found. Returning the original 'loo' object.
## elpd_diff se_diff
## kor_model2 0.0 0.0
## kor_model -0.9 1.2
\normalsize
Comparing against models fit with the drc Package:
Here we will analyze the KOR antagonist data using the drc package and
compare it to the results from the BayesPharma analysis.
We will fix the top to 100 and fit the ic50, hill, and
bottom.
\scriptsize
drc_models <- kor_antag |>
dplyr::group_by(substance_id) |>
dplyr::group_nest() |>
dplyr::mutate(
model = data |>
purrr::map(~drc::drm(
response ~ log_dose,
data = .x,
fct = drc::L.4(fixed = c(NA, NA, 100, NA),
names = c("hill", "bottom", "top", "ic50")))))
drc_models |>
dplyr::mutate(summary = purrr::map(model, broom::tidy, conf.int = TRUE)) |>
tidyr::unnest(summary) |>
dplyr::arrange(term, substance_id) |>
dplyr::select(-data, -model, -curve)
## # A tibble: 12 × 8
## substance_id term estimate std.error statistic p.value conf.low conf.high
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 BTRX_335140 bottom 1.31 19.4 0.0675 9.47e- 1 -40.0 42.6
## 2 BTRX_395750 bottom 29.5 9.40 3.14 7.85e- 3 9.20 49.8
## 3 JNJ bottom -18.1 26.7 -0.681 5.04e- 1 -73.7 37.4
## 4 PF bottom 39.4 30.8 1.28 2.22e- 1 -27.0 106.
## 5 BTRX_335140 hill 4.06 9.20 0.441 6.65e- 1 -15.5 23.7
## 6 BTRX_395750 hill 9.82 164. 0.0600 9.53e- 1 -344. 364.
## 7 JNJ hill 1.17 0.580 2.02 5.69e- 2 -0.0378 2.38
## 8 PF hill 1.13 1.33 0.855 4.08e- 1 -1.73 4.00
## 9 BTRX_335140 ic50 -8.91 0.308 -28.9 1.42e-14 -9.57 -8.26
## 10 BTRX_395750 ic50 -8.97 0.505 -17.8 1.70e-10 -10.1 -7.88
## 11 JNJ ic50 -8.77 0.670 -13.1 2.89e-11 -10.2 -7.37
## 12 PF ic50 -7.96 1.27 -6.29 2.78e- 5 -10.7 -5.23
\normalsize
Displayed below is the comparison of results from drc and
BayesPharma for each parameter of the dose-response curve. Here we
see that the Bayesian method provides a distribution curve as evidence
and has smaller confidence intervals than most of the standard errors
provided by the drc method.
KOR antagonists conditional effects. The blue lines are samples from the `BayesPharma` kor_model posterior distribution, the orange line is the conditional mean, and the purple line is the conditional mean for the `drc` model fit.
maomlab/BayesPharma documentation built on Aug. 24, 2024, 8:45 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Hill Equation {#sec-hill}
In this case study, we are going to reanalyze the dose responses of 4
Kappa Opioid receptor (KOR) antagonists from a study performed by Margolis et
al. (-@Margolis2020-bm) using the BayesPharma
package. Whole cell electrophysiology in acute rat
midbrain slices was used to evaluate the pharmacological properties of
four novel KOR antagonists: BTRX-335140, BTRX-395750, PF-04455242, and
JNJ-67953964.
Originally, the dose-response analysis was performed using
the drc package in R, which implements the minimization of negative
log likelihood function and reduces to least square estimation for a
continuous response. The data were normalized to % baseline agonist response,
measuring the extent of blockade of agonist response in the presence of
different concentrations of antagonist. These data were then fit to a 4-parameter
log-logistic dose response model, setting the top (max agonist response) to
100% and estimating the IC50, its variance, and the bottom (min agonist
response).
Fitting the sigmoid model
Using the BayesPharma package, we can re-fit the sigmoid model with a
negative slope, and fix the top parameter to 100 as the response is
normalized to a no-antagonist baseline.
For the prior, we are going to use a normal distribution because the
response values are continuous. First, we will run the analysis with the
top (max response) parameter prior set to a constant value of 100
because top is normalized to 100 and is the default broad prior for the ic50,
hill, and bottom parameters. Broad priors represent unbiased uncertainty and
provide an opportunity for extreme responses.
The level of informativeness of the prior will affect how much influence the prior has on the model. Here is more information on prior choice recommendations.
\scriptsize
kor_prior <- BayesPharma::sigmoid_antagonist_prior(top = 100) kor_prior ## prior class coef group resp dpar nlpar lb ub source ## normal(-6, 2.5) b ic50 <NA> <NA> user ## normal(-1, 1) b hill <NA> 0.01 user ## constant(100) b top <NA> <NA> user ## normal(0, 0.5) b bottom <NA> <NA> user
\normalsize
Prior predictive checks
Following the Bayesian workflow, before fitting the model, it is good to check
the prior predictive distributions to see if they are compatible with the domain
expertise. So, before running the model, we will verify that the prior
distributions cover a plausible range of values for each parameter. To do this,
we want to sample only from the prior distributions by adding
sample_prior = "only" as an argument to the sigmoid_model
function. We will use the default response distribution of the model
(family = gaussian()).
\scriptsize
kor_sample_prior <- BayesPharma::sigmoid_model( data = kor_antag |> dplyr::select(substance_id, log_dose, response), formula = BayesPharma::sigmoid_antagonist_formula(), prior = kor_prior, init = BayesPharma::sigmoid_antagonist_init(), sample_prior = "only")
\normalsize
And then plot of the prior predictive distributions:
\scriptsize
kor_sample_prior |> BayesPharma::plot_density_distribution()
KOR antagonists prior distribution
\normalsize
To sample from the model, we will use the Stan NUTs Hamiltonian Monte
Carlo, and initialize the parameters to the prior means to help with
model convergence, using the default values of ec50 = -9, hill = -1,
top = 100, bottom = 0.
\scriptsize
kor_model <- BayesPharma::sigmoid_model( data = kor_antag |> dplyr::select(substance_id, log_dose, response), formula = BayesPharma::sigmoid_antagonist_formula( predictors = 0 + substance_id), prior = kor_prior, init = BayesPharma::sigmoid_antagonist_init())
\normalsize
Analyzing model fit
The brms generated model summary shows the formula that the expected
response a is sigmoid function of the log_dose with four parameters, and
a shared Gaussian distribution. Each parameter is dependent on the
substance_id. Since we want to fit a separate model for each substance, we
include a 0 + to indicate that there is no common intercept. The data
consists of 73 data points and the posterior sampling was done in 4
chains each with 8000 steps with 4000 steps of warm-up. The population
effects for each parameter summarize the marginal posterior distributions,
as well as the effective sample size in the bulk and tail. This gives
an indication of the sampling quality, with an ESS of > 500 samples
being good for this type of model.
\scriptsize
## Family: gaussian ## Links: mu = identity; sigma = identity ## Formula: response ~ sigmoid(ic50, hill, top, bottom, log_dose) ## ic50 ~ 0 + substance_id ## hill ~ 0 + substance_id ## top ~ 0 + substance_id ## bottom ~ 0 + substance_id ## Data: data (Number of observations: 73) ## Draws: 4 chains, each with iter = 8000; warmup = 4000; thin = 1; ## total post-warmup draws = 16000 ## ## Regression Coefficients: ## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS ## ic50_substance_idBTRX_335140 -8.84 0.20 -9.20 -8.42 1.00 15027 8248 ## ic50_substance_idBTRX_395750 -8.24 0.41 -8.93 -7.36 1.00 11072 5925 ## ic50_substance_idJNJ -9.15 0.31 -9.77 -8.50 1.00 15946 10770 ## ic50_substance_idPF -6.14 1.08 -7.69 -3.32 1.00 7646 5332 ## hill_substance_idBTRX_335140 -1.47 0.59 -2.85 -0.57 1.00 14309 9824 ## hill_substance_idBTRX_395750 -0.91 0.52 -2.28 -0.26 1.00 10745 7238 ## hill_substance_idJNJ -1.01 0.51 -2.38 -0.41 1.00 14567 11520 ## hill_substance_idPF -0.31 0.24 -0.92 -0.03 1.00 7908 5627 ## bottom_substance_idBTRX_335140 -0.00 0.49 -0.96 0.95 1.00 17746 11293 ## bottom_substance_idBTRX_395750 0.01 0.50 -0.97 0.99 1.00 18472 10902 ## bottom_substance_idJNJ -0.01 0.49 -0.99 0.96 1.00 17703 11821 ## bottom_substance_idPF 0.01 0.50 -0.97 0.99 1.00 16850 11426 ## top_substance_idBTRX_335140 100.00 0.00 100.00 100.00 NA NA NA ## top_substance_idBTRX_395750 100.00 0.00 100.00 100.00 NA NA NA ## top_substance_idJNJ 100.00 0.00 100.00 100.00 NA NA NA ## top_substance_idPF 100.00 0.00 100.00 100.00 NA NA NA ## ## Further Distributional Parameters: ## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS ## sigma 32.17 2.90 27.13 38.38 1.00 14788 11389 ## ## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS ## and Tail_ESS are effective sample size measures, and Rhat is the potential ## scale reduction factor on split chains (at convergence, Rhat = 1).
\normalsize
Traceplot: The model ran without warning messages, meaning there were no parameter
value problems or MCMC conflicts. The bulk and tail ESS indicate high
resolution and stability. The R-hat for each parameter equals 1.00 and
the traceplot shows the chains mixed well, indicating the chains
converged.
\scriptsize
kor_model |> bayesplot::mcmc_trace()
plot of chunk kor-model-traceplot
\normalsize
Compare prior and posterior marginal distributions: Displayed below is a plot for the prior and posterior distributions of the parameters (prior is pink and posterior is teal). This can be useful for comparing the density distribution of the prior and posterior produced by the model:
\scriptsize
BayesPharma::plot_prior_posterior_densities( model = kor_model, title_label = "")
KOR antagonists model, compare prior and posterior distributions for each substance
\normalsize
Displayed below is a plot of the posterior distributions for each parameter with the confidence intervals and mean. This is a useful visual of the model results and can highlight the mode and high-density intervals:
\scriptsize
BayesPharma::plot_posterior_density( kor_model, title_label = "") ## Error in `dplyr::filter()`: ## ℹ In argument: `!...`. ## Caused by error in `.data[[".variable"]]`: ## ! Column `.variable` not found in `.data`.
\normalsize
Displayed below is a plot of a sample of 100 sigmoid dose-response curves from the posterior distribution (purple) and the median quantile intervals:
\scriptsize
BayesPharma::plot_posterior_draws( model = kor_model, title = "")
KOR antagonists, posterior draws
\normalsize
Comparing alternative models
To test the sensitivity of the analysis to the prior, we can re-fit the model with a more informative prior:
\scriptsize
## prior class coef group resp dpar nlpar lb ub source ## normal(-8.5, 0.5) b ic50 <NA> <NA> user ## normal(-1, 0.5) b hill <NA> 0.01 user ## constant(100) b top <NA> <NA> user ## normal(10, 15) b bottom <NA> <NA> user
\normalsize
Re-fitting the model with the kor_prior2 gives:
\scriptsize
## Family: gaussian ## Links: mu = identity; sigma = identity ## Formula: response ~ sigmoid(ic50, hill, top, bottom, log_dose) ## ic50 ~ 0 + substance_id ## hill ~ 0 + substance_id ## top ~ 0 + substance_id ## bottom ~ 0 + substance_id ## Data: data (Number of observations: 73) ## Draws: 4 chains, each with iter = 8000; warmup = 4000; thin = 1; ## total post-warmup draws = 16000 ## ## Regression Coefficients: ## Estimate Est.Error l-95% CI u-95% CI Rhat ## ic50_substance_idBTRX_335140 -8.82 0.21 -9.23 -8.37 1.00 ## ic50_substance_idBTRX_395750 -8.58 0.31 -9.16 -7.94 1.00 ## ic50_substance_idJNJ -8.95 0.31 -9.53 -8.32 1.00 ## ic50_substance_idPF -8.18 0.45 -9.09 -7.30 1.00 ## hill_substance_idBTRX_335140 -1.19 0.37 -1.99 -0.57 1.00 ## hill_substance_idBTRX_395750 -1.06 0.39 -1.90 -0.42 1.00 ## hill_substance_idJNJ -0.85 0.33 -1.65 -0.39 1.00 ## hill_substance_idPF -0.72 0.45 -1.72 -0.07 1.00 ## bottom_substance_idBTRX_335140 1.71 10.47 -19.21 21.94 1.00 ## bottom_substance_idBTRX_395750 15.44 11.10 -7.30 36.02 1.00 ## bottom_substance_idJNJ -2.77 9.82 -23.22 15.35 1.00 ## bottom_substance_idPF 31.22 11.59 7.03 52.03 1.00 ## top_substance_idBTRX_335140 100.00 0.00 100.00 100.00 NA ## top_substance_idBTRX_395750 100.00 0.00 100.00 100.00 NA ## top_substance_idJNJ 100.00 0.00 100.00 100.00 NA ## top_substance_idPF 100.00 0.00 100.00 100.00 NA ## Bulk_ESS Tail_ESS ## ic50_substance_idBTRX_335140 13087 10643 ## ic50_substance_idBTRX_395750 12995 10960 ## ic50_substance_idJNJ 13299 11776 ## ic50_substance_idPF 13156 10947 ## hill_substance_idBTRX_335140 14798 9835 ## hill_substance_idBTRX_395750 13437 9765 ## hill_substance_idJNJ 13283 11670 ## hill_substance_idPF 8776 5539 ## bottom_substance_idBTRX_335140 14068 12061 ## bottom_substance_idBTRX_395750 12671 10835 ## bottom_substance_idJNJ 12510 10615 ## bottom_substance_idPF 10495 10156 ## top_substance_idBTRX_335140 NA NA ## top_substance_idBTRX_395750 NA NA ## top_substance_idJNJ NA NA ## top_substance_idPF NA NA ## ## Further Distributional Parameters: ## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS ## sigma 31.86 2.83 26.91 37.99 1.00 14470 11441 ## ## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS ## and Tail_ESS are effective sample size measures, and Rhat is the potential ## scale reduction factor on split chains (at convergence, Rhat = 1).
\normalsize
Comparing the Two Models Using LOO-Comparison: One way to evaluate the
quality of a model is for each data-point, re-fit the model with remaining
points and evaluate the log probability of the point in the posterior
distribution. Taking the expectation across all points gives the Expected Log
Pointwise predictive Density (ELPD). Since this is computationally challenging
to re-fit the model for each point, if the model fits the data reasonably well,
then the ELPD can be approximated using the Pareto smoothed importance sampling
(PSIS). Using the LOO package, Pareto k value for each data point is computed
that indicate how well that data point is fit by the model, where k < 0.5 is
good, 0.5 <= k < 0.7 is OK, and 0.7 <= k is bad. Evaluating the model for
the KOR antagonists shows that the model fits the data well.
\scriptsize
##
## Computed from 16000 by 73 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -360.5 6.5
## p_loo 6.9 1.1
## looic 721.0 13.1
## ------
## MCSE of elpd_loo is 0.0.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 1.4]).
##
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
\normalsize
Since ELPD is a global measure of model fit, it can be used to compare models.
Using loo_compare from the LOO package returns the elpd_diff and se_diff
for each model relative the model with the lowest ELPD. The kor_model2, the
model with more informative prior, is the preferred model, but not
significantly.
\scriptsize
## No problematic observations found. Returning the original 'loo' object. ## elpd_diff se_diff ## kor_model2 0.0 0.0 ## kor_model -0.9 1.2
\normalsize
Comparing against models fit with the drc Package:
Here we will analyze the KOR antagonist data using the drc package and
compare it to the results from the BayesPharma analysis.
We will fix the top to 100 and fit the ic50, hill, and
bottom.
\scriptsize
drc_models <- kor_antag |> dplyr::group_by(substance_id) |> dplyr::group_nest() |> dplyr::mutate( model = data |> purrr::map(~drc::drm( response ~ log_dose, data = .x, fct = drc::L.4(fixed = c(NA, NA, 100, NA), names = c("hill", "bottom", "top", "ic50"))))) drc_models |> dplyr::mutate(summary = purrr::map(model, broom::tidy, conf.int = TRUE)) |> tidyr::unnest(summary) |> dplyr::arrange(term, substance_id) |> dplyr::select(-data, -model, -curve) ## # A tibble: 12 × 8 ## substance_id term estimate std.error statistic p.value conf.low conf.high ## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> ## 1 BTRX_335140 bottom 1.31 19.4 0.0675 9.47e- 1 -40.0 42.6 ## 2 BTRX_395750 bottom 29.5 9.40 3.14 7.85e- 3 9.20 49.8 ## 3 JNJ bottom -18.1 26.7 -0.681 5.04e- 1 -73.7 37.4 ## 4 PF bottom 39.4 30.8 1.28 2.22e- 1 -27.0 106. ## 5 BTRX_335140 hill 4.06 9.20 0.441 6.65e- 1 -15.5 23.7 ## 6 BTRX_395750 hill 9.82 164. 0.0600 9.53e- 1 -344. 364. ## 7 JNJ hill 1.17 0.580 2.02 5.69e- 2 -0.0378 2.38 ## 8 PF hill 1.13 1.33 0.855 4.08e- 1 -1.73 4.00 ## 9 BTRX_335140 ic50 -8.91 0.308 -28.9 1.42e-14 -9.57 -8.26 ## 10 BTRX_395750 ic50 -8.97 0.505 -17.8 1.70e-10 -10.1 -7.88 ## 11 JNJ ic50 -8.77 0.670 -13.1 2.89e-11 -10.2 -7.37 ## 12 PF ic50 -7.96 1.27 -6.29 2.78e- 5 -10.7 -5.23
\normalsize
Displayed below is the comparison of results from drc and
BayesPharma for each parameter of the dose-response curve. Here we
see that the Bayesian method provides a distribution curve as evidence
and has smaller confidence intervals than most of the standard errors
provided by the drc method.
KOR antagonists conditional effects. The blue lines are samples from the `BayesPharma` kor_model posterior distribution, the orange line is the conditional mean, and the purple line is the conditional mean for the `drc` model fit.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.