--- title: "Simple Emax model fit with Stan" author: "Kenta Yoshida" date: "`r Sys.Date()`" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Simple Emax model fit with Stan} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` ```{r message=FALSE} library(rstanemax) library(dplyr) library(ggplot2) set.seed(12345) ``` This vignette provide an overview of the workflow of Emax model analysis using this package. # Typical workflow ## Model run with `stan_emax` function `stan_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`. ```{r, results="hide"} 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) ``` ```{r} fit.emax ``` `plot()` function shows the estimated Emax model curve with 95% credible intervals of parameters. ```{r plot_example, fig.show='hold'} 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). ```{r plot_example_log, fig.show='hold'} 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. ```{r} class(extract_stanfit(fit.emax)) ``` ## Prediction of response with new exposure data `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`. ```{r} 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. ```{r} 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. ```{r} obs.formatted <- extract_obs_mod_frame(fit.emax) ``` These are particularly useful when you want to plot the estimated Emax curve. ```{r plot_with_pp, fig.show='hold'} 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. ```{r} posterior.fit.emax <- extract_param(fit.emax) posterior.fit.emax ``` # Fix parameter values in Emax model 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. ```{r, results="hide"} 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) ``` ```{r} fit.emax.sigmoidal ``` You can compare the difference of posterior predictions between two models (in this case they are very close to each other): ```{r plot_with_gamma_fix, fig.width = 6, fig.height = 4, fig.show='hold'} 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") ``` # Set covariates You can specify categorical covariates for Emax, EC50, and E0. See help of `stan_emax()` for the details. ```{r, results="hide"} 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) ``` ```{r plot_with_cov, fig.width = 6, fig.height = 4, fig.show='hold'} fit.cov plot(fit.cov) ``` You can extract MCMC samples from raw stanfit and evaluate differences between groups. ```{r compare_emax, fig.show='hold'} 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) ```