diff --git a/NAMESPACE b/NAMESPACE
index 9e0d8b0..916d326 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -8,6 +8,7 @@ S3method(print,mipdeval_results)
 S3method(print,mipdeval_results_bayesian_impact)
 S3method(print,mipdeval_results_shrinkage)
 S3method(print,mipdeval_results_stats_summ)
+export(accuracy)
 export(add_grouping_column)
 export(calculate_bayesian_impact)
 export(calculate_shrinkage)
@@ -17,13 +18,19 @@ export(compare_psn_proseval_results)
 export(group_by_dose)
 export(group_by_time)
 export(install_default_literature_model)
+export(is_accurate)
+export(is_accurate_abs)
+export(is_accurate_rel)
 export(new_ode_model)
 export(parse_psn_proseval_results)
 export(reldiff_psn_execute_results)
 export(reldiff_psn_proseval_results)
 export(run_eval)
+export(stats_summ_options)
 export(vpc_options)
 importFrom(PKPDsim,install_default_literature_model)
 importFrom(PKPDsim,new_ode_model)
 importFrom(rlang,.data)
+importFrom(rlang,caller_arg)
+importFrom(rlang,caller_env)
 importFrom(utils,read.csv)
diff --git a/NEWS.md b/NEWS.md
index 740f29b..459af27 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,7 @@
 # mipdeval (development version)
 
+* Added `accuracy()` statistic for calculating the proportion of predicted drug concentrations that fall within specified absolute and relative error margins. This statistic will only be calculated in `run_eval()` when the absolute and relative error margins have been specified in the `.stats_summ_options` argument (see `stats_summ_options()`). (#30)
+
 # mipdeval 0.1.0
 
 * Initial version.
diff --git a/R/calculate_stats.R b/R/calculate_stats.R
index 7de3a71..fd94c7f 100644
--- a/R/calculate_stats.R
+++ b/R/calculate_stats.R
@@ -3,13 +3,19 @@
 #' @param res output object (`mipdeval_results`) from `run_eval()`, or
 #' `data.frame` with raw results.
 #' @param rounding number of decimals to round to.
+#' @param acc_error_abs,acc_error_rel For calculating [accuracy()]: Positive number
+#'   providing an absolute or relative error margin. The cutoff is exclusive of
+#'   the error margin. When `NULL` (the default), accuracy will not be
+#'   calculated and will return `NA` instead.
 #'
 #' @returns tibble
 #'
 #' @export
 calculate_stats <- function(
     res,
-    rounding = 3
+    rounding = 3,
+    acc_error_abs = NULL,
+    acc_error_rel = NULL
 ) {
   if(inherits(res, "mipdeval_results")) {
     res <- res$results
@@ -23,10 +29,67 @@ calculate_stats <- function(
       rmse = rmse(.data$dv, .data$value),
       nrmse = nrmse(.data$dv, .data$value),
       mpe = mpe(.data$dv, .data$value),
-      mape = mape(.data$dv, .data$value)
+      mape = mape(.data$dv, .data$value),
+      accuracy = dplyr::if_else(
+        !is.null(acc_error_abs) & !is.null(acc_error_rel),
+        accuracy(.data$dv, .data$value, acc_error_abs, acc_error_rel),
+        NA
+      )
     ) |>
     dplyr::mutate(dplyr::across(dplyr::everything(), round, rounding)) |>
     dplyr::as_tibble()
   class(out) <- c("mipdeval_results_stats_summ", class(out))
   out
 }
+
+#' Options for summary statistics
+#'
+#' @param ... These dots are reserved for future extensibility and must be empty.
+#' @inheritParams calculate_stats
+#'
+#' @returns A list.
+#' @export
+stats_summ_options <- function(
+    ...,
+    rounding = 3,
+    acc_error_abs = NULL,
+    acc_error_rel = NULL
+) {
+  rlang::check_dots_empty()
+  check_required_accuracy(acc_error_abs, acc_error_rel)
+  out <- list(
+    rounding = vctrs::vec_assert(
+      rounding, ptype = numeric(), size = 1L, arg = "rounding"
+    ),
+    acc_error_abs = vec_assert_or_null(
+      acc_error_abs, ptype = numeric(), size = 1L, arg = "acc_error_abs"
+    ),
+    acc_error_rel = vec_assert_or_null(
+      acc_error_rel, ptype = numeric(), size = 1L, arg = "acc_error_rel"
+    )
+  )
+  structure(out, class = "mipdeval_stats_summ_options")
+}
+
+check_required_accuracy <- function(
+  acc_error_abs = NULL,
+  acc_error_rel = NULL,
+  call = caller_env()
+) {
+  acc_error_abs_null <- is.null(acc_error_abs) & !is.null(acc_error_rel)
+  acc_error_rel_null <- !is.null(acc_error_abs) & is.null(acc_error_rel)
+  if (!acc_error_abs_null & !acc_error_rel_null) return()
+  cli::cli_abort(
+    message = c(
+      paste0(
+        "{.arg acc_error_abs} and {.arg acc_error_rel} must both be specified ",
+        "to calculate accuracy."
+      ),
+      "x" = dplyr::case_when(
+        acc_error_abs_null ~ "{.arg acc_error_abs} is NULL.",
+        acc_error_rel_null ~ "{.arg acc_error_rel} is NULL."
+      )
+    ),
+    call = call
+  )
+}
diff --git a/R/mipdeval-package.R b/R/mipdeval-package.R
index 3499750..78c0013 100644
--- a/R/mipdeval-package.R
+++ b/R/mipdeval-package.R
@@ -2,7 +2,7 @@
 "_PACKAGE"
 
 ## usethis namespace: start
-#' @importFrom rlang .data
+#' @importFrom rlang .data caller_arg caller_env
 #' @importFrom utils read.csv
 ## usethis namespace: end
 NULL
diff --git a/R/misc.R b/R/misc.R
index d0a1ab0..40b13ec 100644
--- a/R/misc.R
+++ b/R/misc.R
@@ -26,6 +26,25 @@ is_timevarying <- function(.data, .cols) {
   )
 }
 
+#' Assert an argument has known prototype and/or size or is NULL
+#'
+#' @inheritParams vctrs::vec_assert
+#'
+#' @returns Either throws an error or returns `x`, invisibly.
+vec_assert_or_null <- function(
+    x,
+    ptype = NULL,
+    size = NULL,
+    arg = caller_arg(x),
+    call = caller_env()
+) {
+  if (!is.null(x)) {
+    vctrs::vec_assert(x = x, ptype = ptype, size = size, arg = arg, call = call)
+  } else {
+    NULL
+  }
+}
+
 #' Root-mean-squared error
 #'
 #' @param obs observations vector
@@ -68,6 +87,63 @@ mpe <- function (obs, pred) {
   sum((obs - pred)/obs)/length(obs)
 }
 
+#' Accuracy
+#'
+#' Accuracy provides a measure of clinical suitability, defined by whether model
+#' predicted drug concentrations fall within an absolute OR relative error
+#' margin of the measured concentrations.
+#'
+#' @param obs Observations vector.
+#' @param pred Predictions vector.
+#' @param error_abs,error_rel Positive number providing an absolute or relative
+#'   error margin. The cutoff is exclusive of the error margin. Defaults to `0`,
+#'   meaning no predictions fall within the error margin.
+#'
+#' @returns For `is_accurate()`, `is_accurate_abs()`, and `is_accurate_rel()`: A
+#'   logical vector indicating whether or not each predicted drug concentration
+#'   was considered accurate according to the specified absolute or relative
+#'   error margin(s).
+#'
+#'   For `accuracy()`: A single value between 0 and 1 indicating the proportion
+#'   of predicted drug concentrations that fell within the specified absolute
+#'   and relative error margins.
+#'
+#' @examples
+#' # Does the predicted drug concentration fall within 0.5 mg/L error margin?
+#' is_accurate_abs(6, 5, error_abs = 0.5)
+#'
+#' # Does the predicted drug concentration fall within 25% error margin?
+#' is_accurate_rel(6, 5, error_rel = 0.25)
+#'
+#' # Does the predicted drug concentration fall within 0.5 mg/L OR 25%?
+#' is_accurate(6, 5, error_abs = 0.5, error_rel = 0.25)
+#'
+#' # What proportion of predicted drug concentrations fell within 0.5 mg/L OR 25%?
+#' accuracy(rnorm(10, 6), rnorm(10, 5), error_abs = 0.5, error_rel = 0.25)
+#'
+#' @export
+accuracy <- function(obs, pred, error_abs = 0, error_rel = 0) {
+  mean(is_accurate(obs, pred, error_abs, error_rel))
+}
+
+#' @rdname accuracy
+#' @export
+is_accurate <- function(obs, pred, error_abs = 0, error_rel = 0) {
+  is_accurate_abs(obs, pred, error_abs) | is_accurate_rel(obs, pred, error_rel)
+}
+
+#' @rdname accuracy
+#' @export
+is_accurate_abs <- function(obs, pred, error_abs = 0) {
+  abs(pred - obs) < error_abs
+}
+
+#' @rdname accuracy
+#' @export
+is_accurate_rel <- function(obs, pred, error_rel = 0) {
+  (pred/obs > 1 - error_rel) & (pred/obs < 1 + error_rel)
+}
+
 #' Weighted sum-of-squares of residuals
 #'
 #' @inheritParams rmse
@@ -83,3 +159,4 @@ ss <- function(obs, pred, w = NULL) {
   if(sum(w) == 0) return(NA)
   sum(w * (obs - pred)^2)
 }
+
diff --git a/R/run_eval.R b/R/run_eval.R
index 55890e4..23d738f 100644
--- a/R/run_eval.R
+++ b/R/run_eval.R
@@ -28,6 +28,8 @@
 #'   approach has been called "model predictive control (MPC)"
 #'   (www.page-meeting.org/?abstract=9076) and may be more predictive than
 #'   "regular" MAP in some scenarios. Default is `FALSE`.
+#' @param .stats_summ_options Options for summary statistics. This must be the
+#'   result from a call to [stats_summ_options()].
 #' @param .vpc_options Options for VPC simulations. This must be the result from
 #'   a call to [vpc_options()].
 #' @param threads number of threads to divide computations on. Default is 1,
@@ -59,11 +61,15 @@ run_eval <- function(
   weight_prior = 1,
   censor_covariates = TRUE,
   incremental = FALSE,
+  .stats_summ_options = stats_summ_options(),
   .vpc_options = vpc_options(),
   threads = 1,
   progress = TRUE,
   verbose = TRUE
 ) {
+  # Evaluate options at the start so any errors are triggered immediately.
+  .stats_summ_options <- .stats_summ_options
+  .vpc_options <- .vpc_options
 
   if(progress) { # configure progressbars
     progressr::handlers(global = TRUE)
@@ -169,7 +175,12 @@ run_eval <- function(
   # res is NULL when vpc_options(..., vpc_only = TRUE).
   if (!is.null(res)) {
     if(verbose) cli::cli_progress_step("Calculating forecasting statistics")
-    out$stats_summ <- calculate_stats(out)
+    out$stats_summ <- calculate_stats(
+      out,
+      rounding = .stats_summ_options$rounding,
+      acc_error_abs = .stats_summ_options$acc_error_abs,
+      acc_error_rel = .stats_summ_options$acc_error_rel
+    )
     out$shrinkage <- calculate_shrinkage(out)
     out$bayesian_impact <- calculate_bayesian_impact(out)
     cli::cli_progress_done()
diff --git a/man/accuracy.Rd b/man/accuracy.Rd
new file mode 100644
index 0000000..5d4bb38
--- /dev/null
+++ b/man/accuracy.Rd
@@ -0,0 +1,55 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/misc.R
+\name{accuracy}
+\alias{accuracy}
+\alias{is_accurate}
+\alias{is_accurate_abs}
+\alias{is_accurate_rel}
+\title{Accuracy}
+\usage{
+accuracy(obs, pred, error_abs = 0, error_rel = 0)
+
+is_accurate(obs, pred, error_abs = 0, error_rel = 0)
+
+is_accurate_abs(obs, pred, error_abs = 0)
+
+is_accurate_rel(obs, pred, error_rel = 0)
+}
+\arguments{
+\item{obs}{Observations vector.}
+
+\item{pred}{Predictions vector.}
+
+\item{error_abs, error_rel}{Positive number providing an absolute or relative
+error margin. The cutoff is exclusive of the error margin. Defaults to \code{0},
+meaning no predictions fall within the error margin.}
+}
+\value{
+For \code{is_accurate()}, \code{is_accurate_abs()}, and \code{is_accurate_rel()}: A
+logical vector indicating whether or not each predicted drug concentration
+was considered accurate according to the specified absolute or relative
+error margin(s).
+
+For \code{accuracy()}: A single value between 0 and 1 indicating the proportion
+of predicted drug concentrations that fell within the specified absolute
+and relative error margins.
+}
+\description{
+Accuracy provides a measure of clinical suitability, defined by whether model
+predicted drug concentrations fall within an absolute OR relative error
+margin of the measured concentrations.
+}
+\examples{
+# Does the predicted drug concentration fall within 0.5 mg/L error margin?
+is_accurate_abs(6, 5, error_abs = 0.5)
+
+# Does the predicted drug concentration fall within 25\% error margin?
+is_accurate_rel(6, 5, error_rel = 0.25)
+
+# Does the predicted drug concentration fall within 0.5 mg/L OR 25\%?
+is_accurate(6, 5, error_abs = 0.5, error_rel = 0.25)
+
+# What proportion of predicted drug concentrations fell within 0.5 mg/L OR 25\%?
+accuracy(rnorm(10, 6), rnorm(10, 5), error_abs = 0.5, error_rel = 0.25)
+
+}
diff --git a/man/calculate_stats.Rd b/man/calculate_stats.Rd
index 29f44b7..9f15f17 100644
--- a/man/calculate_stats.Rd
+++ b/man/calculate_stats.Rd
@@ -4,13 +4,18 @@
 \alias{calculate_stats}
 \title{Calculate basic statistics, like RMSE, MPE, MAPE for forecasted data}
 \usage{
-calculate_stats(res, rounding = 3)
+calculate_stats(res, rounding = 3, acc_error_abs = NULL, acc_error_rel = NULL)
 }
 \arguments{
 \item{res}{output object (\code{mipdeval_results}) from \code{run_eval()}, or
 \code{data.frame} with raw results.}
 
 \item{rounding}{number of decimals to round to.}
+
+\item{acc_error_abs, acc_error_rel}{For calculating \code{\link[=accuracy]{accuracy()}}: Positive number
+providing an absolute or relative error margin. The cutoff is exclusive of
+the error margin. When \code{NULL} (the default), accuracy will not be
+calculated and will return \code{NA} instead.}
 }
 \value{
 tibble
diff --git a/man/run_eval.Rd b/man/run_eval.Rd
index 49e40fa..1c49a8d 100644
--- a/man/run_eval.Rd
+++ b/man/run_eval.Rd
@@ -19,6 +19,7 @@ run_eval(
   weight_prior = 1,
   censor_covariates = TRUE,
   incremental = FALSE,
+  .stats_summ_options = stats_summ_options(),
   .vpc_options = vpc_options(),
   threads = 1,
   progress = TRUE,
@@ -81,6 +82,9 @@ approach has been called "model predictive control (MPC)"
 (www.page-meeting.org/?abstract=9076) and may be more predictive than
 "regular" MAP in some scenarios. Default is \code{FALSE}.}
 
+\item{.stats_summ_options}{Options for summary statistics. This must be the
+result from a call to \code{\link[=stats_summ_options]{stats_summ_options()}}.}
+
 \item{.vpc_options}{Options for VPC simulations. This must be the result from
 a call to \code{\link[=vpc_options]{vpc_options()}}.}
 
diff --git a/man/stats_summ_options.Rd b/man/stats_summ_options.Rd
new file mode 100644
index 0000000..26cfa64
--- /dev/null
+++ b/man/stats_summ_options.Rd
@@ -0,0 +1,29 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/calculate_stats.R
+\name{stats_summ_options}
+\alias{stats_summ_options}
+\title{Options for summary statistics}
+\usage{
+stats_summ_options(
+  ...,
+  rounding = 3,
+  acc_error_abs = NULL,
+  acc_error_rel = NULL
+)
+}
+\arguments{
+\item{...}{These dots are reserved for future extensibility and must be empty.}
+
+\item{rounding}{number of decimals to round to.}
+
+\item{acc_error_abs, acc_error_rel}{For calculating \code{\link[=accuracy]{accuracy()}}: Positive number
+providing an absolute or relative error margin. The cutoff is exclusive of
+the error margin. When \code{NULL} (the default), accuracy will not be
+calculated and will return \code{NA} instead.}
+}
+\value{
+A list.
+}
+\description{
+Options for summary statistics
+}
diff --git a/man/vec_assert_or_null.Rd b/man/vec_assert_or_null.Rd
new file mode 100644
index 0000000..46e1e3d
--- /dev/null
+++ b/man/vec_assert_or_null.Rd
@@ -0,0 +1,39 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/misc.R
+\name{vec_assert_or_null}
+\alias{vec_assert_or_null}
+\title{Assert an argument has known prototype and/or size or is NULL}
+\usage{
+vec_assert_or_null(
+  x,
+  ptype = NULL,
+  size = NULL,
+  arg = caller_arg(x),
+  call = caller_env()
+)
+}
+\arguments{
+\item{x}{A vector argument to check.}
+
+\item{ptype}{Prototype to compare against. If the prototype has a
+class, its \code{\link[vctrs:vec_ptype]{vec_ptype()}} is compared to that of \code{x} with
+\code{identical()}. Otherwise, its \code{\link[=typeof]{typeof()}} is compared to that of
+\code{x} with \code{==}.}
+
+\item{size}{A single integer size against which to compare.}
+
+\item{arg}{Name of argument being checked. This is used in error
+messages. The label of the expression passed as \code{x} is taken as
+default.}
+
+\item{call}{The execution environment of a currently
+running function, e.g. \code{caller_env()}. The function will be
+mentioned in error messages as the source of the error. See the
+\code{call} argument of \code{\link[rlang:abort]{abort()}} for more information.}
+}
+\value{
+Either throws an error or returns \code{x}, invisibly.
+}
+\description{
+Assert an argument has known prototype and/or size or is NULL
+}
diff --git a/tests/testthat/test-calculate_stats.R b/tests/testthat/test-calculate_stats.R
new file mode 100644
index 0000000..954763e
--- /dev/null
+++ b/tests/testthat/test-calculate_stats.R
@@ -0,0 +1,23 @@
+# stats_summ_options() --------------------------------------------------------
+test_that("stats_summ_options() works", {
+  actual <- stats_summ_options(
+    rounding = 2, acc_error_abs = 3, acc_error_rel = 4
+  )
+  expected <- list(rounding = 2, acc_error_abs = 3, acc_error_rel = 4)
+  expect_s3_class(actual, "mipdeval_stats_summ_options")
+  expect_identical(unclass(actual), expected)
+})
+
+test_that("stats_summ_options() throws error for incomplete accuracy args", {
+  expect_error(
+    stats_summ_options(rounding = 2, acc_error_abs = NULL, acc_error_rel = 4),
+    "`acc_error_abs` is NULL"
+  )
+  expect_error(
+    stats_summ_options(rounding = 2, acc_error_abs = 3, acc_error_rel = NULL),
+    "`acc_error_rel` is NULL"
+  )
+  expect_no_error(
+    stats_summ_options(rounding = 2, acc_error_abs = NULL, acc_error_rel = NULL)
+  )
+})
diff --git a/tests/testthat/test-error_metrics.R b/tests/testthat/test-error_metrics.R
new file mode 100644
index 0000000..341c6eb
--- /dev/null
+++ b/tests/testthat/test-error_metrics.R
@@ -0,0 +1,45 @@
+test_that("Root mean squared error", {
+  pred <- c(7, 10, 12, 10, 10, 8, 7,  8, 11, 13, 10, 8)
+  obs  <- c(6, 10, 14, 16,  7, 5, 5, 13, 12, 13,  8, 5)
+  expect_equal(round(rmse(obs, pred), 3), 2.915)
+})
+
+test_that("Normalized root mean squared error", {
+  pred <- c(7, 10, 12, 10, 10, 8, 7,  8, 11, 13, 10, 8)
+  obs  <- c(6, 10, 14, 16,  7, 5, 5, 13, 12, 13,  8, 5)
+  mean_obs <- mean(obs)
+  expect_equal(nrmse(obs, pred), sqrt(8.5)/ mean_obs)
+})
+
+test_that("Mean absolute percent error", {
+  pred <- c(7, 10, 12, 10, 10, 8, 7,  8, 11, 13, 10, 8)
+  obs  <- c(6, 10, 14, 16,  7, 5, 5, 13, 12, 13,  8, 5)
+  expect_equal(round(mape(obs, pred), 3), 0.286)
+})
+
+test_that("Accuracy", {
+  obs <- 6
+  pred <- 5
+
+  # No observations should be considered accurate when the error margin is 0
+  # (the default):
+  expect_false(is_accurate_abs(5, 5, error_abs = 0))
+  expect_false(is_accurate_rel(5, 5, error_rel = 0))
+  expect_false(is_accurate(obs, pred, error_abs = 0, error_rel = 0))
+  expect_false(is_accurate(obs, pred))
+
+  expect_true(is_accurate(obs, pred, error_abs = 2, error_rel = 0))
+  expect_false(is_accurate(obs, pred, error_abs = 0.5, error_rel = 0))
+  expect_true(is_accurate(obs, pred, error_abs = 0, error_rel = 0.5))
+  expect_false(is_accurate(obs, pred, error_abs = 0, error_rel = 0.05))
+
+  expect_equal(
+    accuracy(
+      c(1, 2, 3, 4, 5),
+      c(1.1, 2.5, 3.1, 4.5, 5.1),
+      error_abs = 0.2,
+      error_rel = 0.05
+    ),
+    0.6
+  )
+})
diff --git a/tests/testthat/test-run_eval.R b/tests/testthat/test-run_eval.R
index 4f50953..69493ed 100644
--- a/tests/testthat/test-run_eval.R
+++ b/tests/testthat/test-run_eval.R
@@ -11,6 +11,7 @@ test_that("Basic run with vanco data + model works", {
     ruv = mod_obj$ruv,
     fixed = mod_obj$fixed,
     censor_covariates = FALSE, # shouldn't matter, since no time-varying covs
+    .stats_summ_options = stats_summ_options(acc_error_abs = 0.5, acc_error_rel = 0.25),
     .vpc_options = vpc_options(skip = TRUE),
     progress = FALSE
   )
@@ -43,6 +44,7 @@ test_that("Basic run with vanco data + model works", {
     model = "pkvancothomson",
     data = nm_vanco,
     censor_covariates = FALSE, # shouldn't matter, since no time-varying covs
+    .stats_summ_options = stats_summ_options(acc_error_abs = 0.5, acc_error_rel = 0.25),
     .vpc_options = vpc_options(skip = TRUE),
     progress = FALSE
   )
@@ -148,12 +150,10 @@ test_that("Incremental Bayes method works", {
 "map_ipred", "pred", "pred"),
       apriori = c(0, 1, 0, 1, 0, 1),
       rmse = c(4.051, 8.299, 1.96, 4.319, 4.789, 8.299),
-      nrmse = c(0.411,
-    0.232, 0.199, 0.121, 0.486, 0.232),
-      mpe = c(-0.339, -0.004,
-    -0.045, -0.002, -0.415, -0.004),
-      mape = c(0.445, 0.246, 0.166,
-    0.118, 0.558, 0.246)
+      nrmse = c(0.411, 0.232, 0.199, 0.121, 0.486, 0.232),
+      mpe = c(-0.339, -0.004, -0.045, -0.002, -0.415, -0.004),
+      mape = c(0.445, 0.246, 0.166, 0.118, 0.558, 0.246),
+      accuracy = c(NA_real_, NA_real_, NA_real_, NA_real_, NA_real_, NA_real_)
     ),
     class = c("mipdeval_results_stats_summ","tbl_df", "tbl", "data.frame"),
     row.names = c(NA, -6L))