knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
library(rstanemax) library(dplyr) library(ggplot2) set.seed(12345)
This vignette provide an overview of the workflow of Emax model analysis using this package.
stan_emax
functionstan_emax()
is the main function of this package to perform Emax model analysis on the data.
This function requires minimum two input arguments - formula
and data
.
In the formula
argument, you will specify which columns of data
will be used as exposure and response data, in a format similar to stats::lm()
function, e.g. response ~ exposure
.
data(exposure.response.sample) fit.emax <- stan_emax(response ~ exposure, data = exposure.response.sample, # the next line is only to make the example go fast enough chains = 2, iter = 1000, seed = 12345)
fit.emax
plot()
function shows the estimated Emax model curve with 95% credible intervals of parameters.
plot(fit.emax)
Output of plot()
function (for stanemax
object) is a ggplot
object, so you can apply additional settings as you would do in ggplot2
.
Here is an example of using log scale for x axis (note that exposure == 0 is hanging at the very left, making the curve a bit weird).
plot(fit.emax) + scale_x_log10() + expand_limits(x = 1)
Raw output from rstan
is stored in the output variable, and you can access it with extract_stanfit()
function.
class(extract_stanfit(fit.emax))
posterior_predict()
function allows users to predict the response using new exposure data.
If newdata
is not provided, the function returns the prediction on the exposures in original data.
The default output is a matrix of posterior predictions, but you can also specify "dataframe" or "tibble" that contain posterior predictions in a long format.
See help of rstanemax::posterior_predict()
for the description of two predictions, respHat
and response
.
response.pred <- posterior_predict(fit.emax, newdata = c(0, 100, 1000), returnType = "tibble") response.pred %>% select(mcmcid, exposure, respHat, response)
You can also get quantiles of predictions with posterior_predict_quantile()
function.
resp.pred.quantile <- posterior_predict_quantile(fit.emax, newdata = seq(0, 5000, by = 100)) resp.pred.quantile
Input data can be obtained in a same format with extract_obs_mod_frame()
function.
obs.formatted <- extract_obs_mod_frame(fit.emax)
These are particularly useful when you want to plot the estimated Emax curve.
ggplot(resp.pred.quantile, aes(exposure, ci_med)) + geom_line() + geom_ribbon(aes(ymin=ci_low, ymax=ci_high), alpha = .5) + geom_ribbon(aes(ymin=pi_low, ymax=pi_high), alpha = .2) + geom_point(data = obs.formatted, aes(y = response)) + labs(y = "response")
Posterior draws of Emax model parameters can be extracted with extract_param()
function.
posterior.fit.emax <- extract_param(fit.emax) posterior.fit.emax
You can fix parameter values in Emax model for Emax, E0 and/or gamma (Hill coefficient).
See help of stan_emax()
for the details.
The default is to fix gamma at 1 and to estimate Emax and E0 from data.
Below is the example of estimating gamma from data.
data(exposure.response.sample) fit.emax.sigmoidal <- stan_emax(response ~ exposure, data = exposure.response.sample, gamma.fix = NULL, # the next line is only to make the example go fast enough chains = 2, iter = 1000, seed = 12345)
fit.emax.sigmoidal
You can compare the difference of posterior predictions between two models (in this case they are very close to each other):
exposure_pred <- seq(min(exposure.response.sample$exposure), max(exposure.response.sample$exposure), length.out = 100) pred1 <- posterior_predict_quantile(fit.emax, exposure_pred) %>% mutate(model = "Emax") pred2 <- posterior_predict_quantile(fit.emax.sigmoidal, exposure_pred) %>% mutate(model = "Sigmoidal Emax") pred <- bind_rows(pred1, pred2) ggplot(pred, aes(exposure, ci_med, color = model, fill = model)) + geom_line() + geom_ribbon(aes(ymin=ci_low, ymax=ci_high), alpha = .3) + geom_ribbon(aes(ymin=pi_low, ymax=pi_high), alpha = .1, color = NA) + geom_point(data=exposure.response.sample, aes(exposure, response), color = "black", fill = NA, size=2) + labs(y = "response")
You can specify categorical covariates for Emax, EC50, and E0.
See help of stan_emax()
for the details.
data(exposure.response.sample.with.cov) test.data <- mutate(exposure.response.sample.with.cov, SEX = ifelse(cov2 == "B0", "MALE", "FEMALE")) fit.cov <- stan_emax(formula = resp ~ conc, data = test.data, param.cov = list(emax = "SEX"), # the next line is only to make the example go fast enough chains = 2, iter = 1000, seed = 12345)
fit.cov
plot(fit.cov)
You can extract MCMC samples from raw stanfit and evaluate differences between groups.
fit.cov.posterior <- extract_param(fit.cov) emax.posterior <- fit.cov.posterior %>% select(mcmcid, SEX, emax) %>% tidyr::pivot_wider(names_from = SEX, values_from = emax) %>% mutate(delta = FEMALE - MALE) ggplot2::qplot(delta, data = emax.posterior, bins = 30) + ggplot2::labs(x = "emax[FEMALE] - emax[MALE]") # Credible interval of delta quantile(emax.posterior$delta, probs = c(0.025, 0.05, 0.5, 0.95, 0.975)) # Posterior probability of emax[FEMALE] < emax[MALE] sum(emax.posterior$delta < 0) / nrow(emax.posterior)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.