---
title: "Dose-Response Analysis"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{Dose-Response Analysis}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
resource_files:
- figures/logo.svg
- figures/logo-dark.svg
---
```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>",
fig.width = 6, fig.height = 4, fig.align = "center")
options(width = 84)
library(cellreportR)
```
`cr_dose_response()` fits three standard models to unit-level summaries of one
channel across exposure levels:
* **4PL** — four-parameter logistic in `log10(dose)`, with a free lower
asymptote;
* **3PL** — the same shape with the lower asymptote fixed at zero;
* **linear** — `y ~ x`, for a readout that has no plateau in the range
measured.
## Setting up an exposure axis
Any `cr_experiment` whose `design` carries a numeric `dose` column is
eligible. The demonstration experiment has one already, but its levels are
tied to treatment names rather than spread along an axis, so we widen it into
six exposure levels.
```{r setup-doses}
exp <- cr_example_experiment(seed = 3, n_cells_per_well = 60)
exp$design$dose <- c(Untreated = 1,
CompoundB = 10,
CompoundC = 30,
CompoundA_low = 100,
CompoundA_high = 300,
PosControl = 1000)[exp$design$treatment]
table(exp$design$dose)
```
Doses must be strictly positive when `log_dose = TRUE`; a zero-dose vehicle
level has no place on a log axis, and the function drops such units rather
than taking `log10(0)`.
## Fitting
`marker_2` is the readout that falls as exposure rises, which is the shape a
half-maximal *inhibitory* concentration describes.
```{r fit}
fit <- cr_dose_response(exp,
channel = "marker_2",
model = "4pl",
log_dose = TRUE)
fit$model
fit$params
```
`a` is the upper asymptote, `d` the lower, `e` the inflection point on the
`log10(dose)` scale and `b` the slope there.
## IC50 / EC50
```{r ic50}
cr_ic50(fit)
```
The estimate is back-transformed to the dose scale, with the interval carried
through from the standard error of `e`.
## Plotting
```{r plot-dr}
cr_plot_dose_response(fit)
```
## Comparing models
Fixing the lower asymptote at zero moves the inflection point, and with it the
IC50 — a reminder that the number is a property of the model as much as of the
data:
```{r compare-models}
fit_3pl <- cr_dose_response(exp, channel = "marker_2", model = "3pl")
rbind(
cbind(model = "4pl", cr_ic50(fit)),
cbind(model = "3pl", cr_ic50(fit_3pl))
)
```
## Checking that the fit converged
`cr_dose_response()` calls `stats::nls()` with `warnOnly = TRUE`, so a fit that
runs out of iterations is **returned rather than raised**. Two checks catch it:
```{r converged}
c(model = fit$model, any_na_se = anyNA(fit$params$std_error))
```
* `fit$model` is `"linear"` instead of the model you asked for when `nls()`
failed outright and the linear fallback took over.
* `std_error` is `NA` for every parameter when `nls()` returned without
converging. The point estimates from such a fit should not be quoted, and
neither should the IC50 derived from them.
A readout that *rises* with exposure is the usual cause. The logistic is
fitted from starting values that describe a falling curve, so an increasing
response often does not converge. Fit the channel that falls, or use the
linear model:
```{r increasing}
lin <- cr_dose_response(exp, channel = "marker_1", model = "linear")
lin$params
```
Here the slope is the change in the unit-level median of `marker_1` per
tenfold increase in exposure.
## Where this fits
`vignette("end-to-end")` runs the full screening pipeline, in which the
exposure axis is handled as a set of contrasts against the vehicle of each
unit's own batch rather than as a fitted curve. Curve fitting is the right
tool when the exposure levels are dense enough to define a shape; contrasts
are the right tool when they are not.
