Fitting the Emax model with brms
https://www.kristianbrock.com/post/emax-intro/
library(tidyverse)
library(brms)
df <- tibble(
Exposure = c(10, 25, 50, 75, 100, 150, 300, 400),
Response = c(0.03, 0.04, 0.15, 0.12, 0.25, 0.43, 0.54, 0.47)
)
df %>%
ggplot(aes(x = Exposure, y = Response)) +
geom_point() +
geom_smooth(
method = "glm", se = F,
method.args = list(family = "binomial")
)
## `geom_smooth()` using formula 'y ~ x'
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
model <- bf(
Response ~ E0 + Emax * Exposure^N / (ED50^N + Exposure^N),
E0 + N + Emax + ED50 ~ 1,
nl = TRUE
)
get_prior(model, data = df)
p <- prior(normal(0, .5), nlpar = "E0") +
prior(normal(0, .5), nlpar = "Emax") +
prior(normal(150, 50), nlpar = "ED50", lb = 0) +
prior(normal(0, 1), nlpar = "N", lb = 0)
fit <- brm(
model,
prior = p,
data = df,
cores = 4,
control = list(adapt_delta = .99)
)
## Compiling the C++ model
## Trying to compile a simple C file
## Start sampling
summary(fit)
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: Response ~ E0 + Emax * Exposure^N/(ED50^N + Exposure^N)
## E0 ~ 1
## N ~ 1
## Emax ~ 1
## ED50 ~ 1
## Data: df (Number of observations: 8)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## E0_Intercept -0.01 0.09 -0.22 0.13 1.00 1018 939
## N_Intercept 1.39 0.55 0.53 2.64 1.00 932 1100
## Emax_Intercept 0.67 0.19 0.37 1.11 1.00 980 1126
## ED50_Intercept 142.34 39.34 75.32 228.19 1.00 1332 1436
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.08 0.04 0.04 0.19 1.01 1005 912
##
## Samples were drawn 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).
plot(conditional_effects(fit), points = TRUE)