fit_MuSyC_score_by_dose: Fit the MuSyC synergy model by dose

In momeara/MPStats: Tools for analyzing Morphological Profiling data

Fit the MuSyC synergy model by dose

Description

Fit the MuSyC synergy model by dose

Usage

fit_MuSyC_score_by_dose(
  well_scores,
  group_vars = vars(compound),
  C1_prior = brms::prior(normal(0.5, 0.5), nlpar = "C1", lb = 0),
  C2_prior = brms::prior(normal(0.5, 0.5), nlpar = "C2", lb = 0),
  s1_prior = brms::prior(normal(1, 3), nlpar = "s1", lb = -0.1),
  s2_prior = brms::prior(normal(1, 3), nlpar = "s2", lb = -0.1),
  log10alpha_prior = brms::prior(normal(0, 0.5), nlpar = "log10alpha", lb = 0),
  E0_prior = brms::prior(beta(1, 1), nlpar = "E0", lb = 0, ub = 1),
  E1_prior = brms::prior(beta(1, 1), nlpar = "E1", lb = 0, ub = 1),
  E2_prior = brms::prior(beta(1, 1), nlpar = "E2", lb = 0, ub = 1),
  E3_alpha_prior = brms::prior(normal(0, 2), nlpar = "E3alpha", lb = 0),
  C1_init = function() {     as.array(runif(1, 0, 2)) },
  C2_init = function() {     as.array(runif(1, 0, 2)) },
  s1_init = function() {     as.array(runif(1, -0.1, 2)) },
  s2_init = function() {     as.array(runif(1, -0.1, 2)) },
  log10alpha_init = function() {     as.array(0) },
  E0_init = function() {     as.array(rbeta(1, 1, 1)) },
  E1_init = function() {     as.array(rbeta(1, 1, 1)) },
  E2_init = function() {     as.array(rbeta(1, 1, 1)) },
  E3_alpha_init = function() {     as.array(rnorm(1, 0, 2)) },
  combine = FALSE,
  verbose = FALSE,
  iter = 8000,
  cores = 4,
  stan_model_args = list(verbose = FALSE),
  control = list(adapt_delta = 0.99, max_treedepth = 12),
  model_evaluation_criteria = c("loo", "bayes_R2"),
  ...
)

Arguments

well_scores	data.frame with columns: dose1, dose2, n_positive, count, [<group_vars>]
group_vars	quosures list dplyr::vars(...) columns to over when fitting synergy model
C1_prior	prior distribution for Ed when d1=d1_IC50, d2=0
C2_prior	prior distribution for Ed when d1=0, d2=d2_IC50
s1_prior	prior distribution for d(Ed)/d(d1) when d1=d1_IC50, d2=0
s2_prior	prior distribution for d(Ed)/d(d2) when d1=0, d2=d2_IC50
log10alpha_prior	prior distribution for alpha synergy parameter
E0_prior	prior distribution for Ed when d1=0, d2=0
E1_prior	prior distribution for Ed when d1=Inf, d2=0
E2_prior	prior distribution for Ed when d1=0, d2=Inf
E3_alpha_prior	prior distribution for Ed scaled by alpha when d1=Inf, d2=Inf
C1_init	initial sampling distribution for the C1 parameter
C2_init	initial sampling distribution for the C2 parameter
s1_init	initial sampling distribution for the s1 parameter
s2_init	initial sampling distribution for the s2 parameter
log10alpha_init	intial sampling distribution for the alpha parameter
E0_init	initial sampling distribution for the E0 parameter
E1_init	initial sampling distribution for the E1 parameter
E2_init	initial sampling distribution for the E2 parameter
E3_alpha_init	initial sampling distribution for the E3 parameter
combine	combine the grouped models into a single brms model
verbose	verbose output
iter	number of stan NUTS sampling steps
cores	number of cores used for sampling
stan_model_args	stan model arguments
control	stan control arguments The bernoulli_inf(n_positive / count) = Ed ~ MuSyC(d1, d2, C_params, E_params, s_params, alpha) To improve numeric stability, the d1 and d2 and C1 and C2 variables are scaled to improve numeric stability: d1 = dose1/max(dose1) d2 = dose2/max(dose2) drug1_IC50 = C1 * max(dose1) drug2_IC50 = C2 * max(dose2) Functional form: Ed ~ ( C1^h1 * C2^h2 * E0 + d1^h1 * C2^h2 * E1 + C1^h1 * d2^h2 * E2 + d1^h1 * d2^h2 * E3 * alpha ) / ( C1^h1 * C2^h2 + d1^h1 * C2^h2 + C1^h1 * d2^h2 + d1^h1 * d2^h2 * alpha ) ############################################## # Proof of the definitions of the parameters # ############################################## Claim: When d1=0 and d2=0 then Ed = E0 Ed = (C1^h1 * C2^h2 * E0) / (C1^h1 * C2^h2) = E0 Claim: When d1=0 and d2 -> Inf then Ed = E2 Ed = (C2^h2 * E0 + d2^h2 * E2) / (C2^h2 + d2^h2) = (d2^h2 * E2) / (d2^h2) = E2 When d1>0 and d2 -> Inf then Ed Ed = (C1^h1 * C2^h2 * E0 + C1^h1 * C2^h2 * E2) / (C1^h1 * C2^h2 + C1^h1 * C2^h2) = (E0 + E2) / 2

momeara/MPStats documentation built on July 19, 2022, 3:34 p.m.