yoshidk6 / rstanemax

Simple Emax model fit with Stan
https://yoshidk6.github.io/rstanemax/
GNU General Public License v3.0
5 stars 1 forks source link

Add instruction on calculating credible intervals for EC95 etc. #42

Open yoshidk6 opened 3 years ago

yoshidk6 commented 3 years ago

Not sure where is the best place, maybe vignette? Wrap together with instruction on how to extract stanfit objects?


First extract MCMC samples

https://yoshidk6.github.io/rstanemax/articles/emaxmodel.html#set-covariates-1

extract_stanfit(fit) %>% 
  rstan::extract(pars = c("ec50", "gamma"))

Calculate EC95 from MCMC samples

From this (https://yoshidk6.github.io/rstanemax/reference/stan_emax.html), unnamed

, one can solve for 0.95 = EC95^gamma / (ec50^gamma + EC95^gamma) which would be EC95 = EC50 * (19)^(1/gamma)

Calculate credible intervals

Now that you have posterior samples of EC95 (by default 4000 of them), one can calculate percentiles to obtain credible intervals

Code example

library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(rstanemax)
#> Loading required package: Rcpp

data(exposure.response.sample)

fit.emax <- stan_emax(response ~ exposure, data = exposure.response.sample)
#> 
#> SAMPLING FOR MODEL 'emax' NOW (CHAIN 1).
#> Chain 1: 
#> ............

as.data.frame(extract_stanfit(fit.emax), pars = c("ec50", "gamma")) %>% 
  mutate(ec95 = `ec50[1]` * (19)^(1/gamma)) %>% 
  summarize(P05 = quantile(ec95, probs = 0.05),
            P50 = quantile(ec95, probs = 0.5),
            P95 = quantile(ec95, probs = 0.95))
#>        P05      P50      P95
#> 1 899.3797 1368.342 2167.361

Created on 2021-06-14 by the reprex package (v0.3.0)

Note

This is just as simple as multiplying the CI of EC50 by 19.
The method described above is general one and applicable to more complex cases, e.g. to calculate the credible interval of the concentration that is needed to achieve certain response value (instead of the XX% of Emax). In such case the equation 0.95 = EC95^gamma / (ec50^gamma + EC95^gamma) needs to be replaced