diff --git a/DESCRIPTION b/DESCRIPTION
index 9775d40..46e9a4d 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: bmstate
 Type: Package
 Title: Bayesian multistate modeling
-Version: 0.2.8
+Version: 0.2.9
 Authors@R:
     c(person(given = "Juho",
            family = "Timonen",
diff --git a/R/DosingData.R b/R/DosingData.R
index 835f9df..3c2ef1b 100644
--- a/R/DosingData.R
+++ b/R/DosingData.R
@@ -136,7 +136,9 @@ PSSDosingData <- R6::R6Class(
     #' and \code{upper}.
     #' @param subject_df Subject data frame.
     #' @param max_num_subjects Max number of subjects to plot.
-    plot = function(df_fit = NULL, subject_df = NULL, max_num_subjects = 12) {
+    #' @param subject_ids Which subjects to plot?
+    plot = function(df_fit = NULL, subject_df = NULL, max_num_subjects = 12,
+                    subject_ids = NULL) {
       dos <- self$as_data_frame()
       fitcolor <- "steelblue"
       if (is.null(max_num_subjects)) {
@@ -144,7 +146,13 @@ PSSDosingData <- R6::R6Class(
       }
       checkmate::assert_integerish(max_num_subjects, len = 1, lower = 1)
       if (self$num_subjects() > max_num_subjects) {
-        sid <- sample(unique(dos$subject_id), max_num_subjects)
+        if (is.null(subject_ids)) {
+          sid <- sample(unique(dos$subject_id), max_num_subjects)
+        } else {
+          checkmate::assert_character(subject_ids)
+          sid <- subject_ids
+        }
+
         dos <- dos |> dplyr::filter(.data$subject_id %in% sid)
         if (!is.null(df_fit)) {
           df_fit <- df_fit |> dplyr::filter(.data$subject_id %in% sid)
diff --git a/R/JointData.R b/R/JointData.R
index ad84152..bfee466 100644
--- a/R/JointData.R
+++ b/R/JointData.R
@@ -66,13 +66,14 @@ JointData <- R6::R6Class(
     #' @description Plot dosing
     #' @param df_fit Fit data frame.
     #' @param max_num_subjects Max number of subjects to plot.
+    #' @param subject_ids Which subjects to plot?
     #' @return a \code{ggplot}
-    plot_dosing = function(df_fit = NULL, max_num_subjects = 12) {
+    plot_dosing = function(df_fit = NULL, max_num_subjects = 12, subject_ids = NULL) {
       if (is.null(self$dosing)) {
         stop("No dosing data.")
       }
       sdf <- self$paths$subject_df
-      self$dosing$plot(df_fit, sdf, max_num_subjects)
+      self$dosing$plot(df_fit, sdf, max_num_subjects, subject_ids)
     }
   )
 )
diff --git a/R/MultistateModel.R b/R/MultistateModel.R
index 245b054..632e933 100644
--- a/R/MultistateModel.R
+++ b/R/MultistateModel.R
@@ -47,8 +47,8 @@ MultistateModel <- R6::R6Class("MultistateModel",
     categorical = NULL,
     normalizer_locations = NULL,
     normalizer_scales = NULL,
-    auc_normalizer_loc = 300,
-    auc_normalizer_scale = 100,
+    auc_normalizer_loc = 2000,
+    auc_normalizer_scale = 1000,
     n_grid = NULL,
     simulate_log_hazard_multipliers = function(df_subjects, beta) {
       ts <- self$target_states()
diff --git a/R/MultistateModelFit.R b/R/MultistateModelFit.R
index 78f1862..35d3d1f 100644
--- a/R/MultistateModelFit.R
+++ b/R/MultistateModelFit.R
@@ -169,8 +169,10 @@ MultistateModelFit <- R6::R6Class("MultistateModelFit",
     #' @param n_prev number of previous doses to show fit for
     #' @param oos Out-of-sample subjects?
     #' @param ci_alpha Width of central credible interval.
+    #' @param subject_ids Which subjects to plot?
     plot_pk = function(max_num_subjects = 12, oos = FALSE, data = NULL, L = 100,
-                       timescale = 24, n_prev = 3, ci_alpha = 0.9) {
+                       timescale = 24, n_prev = 3, ci_alpha = 0.9,
+                       subject_ids = NULL) {
       pksim <- self$simulate_pk(oos, data, L, timescale, n_prev)
       if (is.null(data)) {
         data <- self$data
@@ -180,16 +182,15 @@ MultistateModelFit <- R6::R6Class("MultistateModelFit",
       } else {
         capt <- paste0("Point-estimate fit")
       }
-      checkmate::assert_number(ci_alpha, lower = 0, upper = 1)
-      av <- (1 - ci_alpha) / 2
-      pksim <- pksim |>
-        dplyr::group_by(.data$subject_id, .data$time) |>
-        dplyr::summarise(q = list(quantile(.data$val, probs = c(av / 2, 0.5, 1 - av / 2))), .groups = "drop") |>
-        tidyr::unnest_wider(q, names_sep = "_")
-      colnames(pksim)[3:5] <- c("lower", "val", "upper")
-
-      data$plot_dosing(df_fit = pksim, max_num_subjects = max_num_subjects) +
-        labs(caption = capt)
+
+      # Slow for some reason
+      pksim <- pksim_to_quantiles(pksim, ci_alpha)
+
+      # Create plot
+      data$plot_dosing(
+        df_fit = pksim, max_num_subjects = max_num_subjects,
+        subject_ids = subject_ids
+      ) + labs(caption = capt)
     },
 
     #' @description Plot baseline hazard distribution
diff --git a/R/PKModel.R b/R/PKModel.R
index 997156d..4136bb4 100644
--- a/R/PKModel.R
+++ b/R/PKModel.R
@@ -102,16 +102,18 @@ PKModel <- R6::R6Class("PKModel",
       conc
     },
 
-    #' @description Compute steady-state area under curve over one dosing interval
+    #' @description Compute steady-state area under concentration curve over
+    #' one dosing interval
     #'
     #' @param theta parameter values
     #' @param dose dose
     #' @return A numeric value
     compute_ss_auc = function(theta, dose) {
       CL <- theta[2]
+      V2 <- theta[3]
       checkmate::assert_number(CL, lower = 0)
       checkmate::assert_number(dose, lower = 0)
-      dose / CL
+      dose / (CL * V2)
     },
 
     #' @description Simulate data with many subjects
diff --git a/R/stan.R b/R/stan.R
index a665556..aff4304 100644
--- a/R/stan.R
+++ b/R/stan.R
@@ -101,8 +101,9 @@ fit_stan <- function(model, data,
   if (set_normalizers) {
     model$set_normalizers(data)
     if (!is.null(data$dosing)) {
-      mu_CL <- exp(-2)
-      aaa <- data$dosing$dose_ss / mu_CL
+      mu_CL <- exp(-2) # should match msm.stan
+      mu_V2 <- exp(-2) # should match msm.stan
+      aaa <- data$dosing$dose_ss / (mu_CL * mu_V2)
       loc <- mean(aaa)
       sca <- stats::sd(aaa)
       model$set_auc_normalizers(loc, sca)
diff --git a/R/utils.R b/R/utils.R
index ec5788a..7210e14 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -1,3 +1,15 @@
+# Helper
+pksim_to_quantiles <- function(sim, ci_alpha) {
+  checkmate::assert_number(ci_alpha, lower = 0, upper = 1)
+  av <- (1 - ci_alpha) / 2
+  sim <- sim |>
+    dplyr::group_by(.data$subject_id, .data$time) |>
+    dplyr::summarise(q = list(quantile(.data$val, probs = c(av / 2, 0.5, 1 - av / 2))), .groups = "drop") |>
+    tidyr::unnest_wider(q, names_sep = "_")
+  colnames(sim)[3:5] <- c("lower", "val", "upper")
+  sim
+}
+
 # Check data frame columns
 check_columns <- function(df, needed_columns) {
   checkmate::assert_data_frame(df)
diff --git a/inst/stan/msm.stan b/inst/stan/msm.stan
index 742c524..357bfa8 100644
--- a/inst/stan/msm.stan
+++ b/inst/stan/msm.stan
@@ -427,8 +427,8 @@ transformed parameters {
       MAX_CONC
     );
 
-    // Auc at steady state
-    ss_auc[1] = dose_ss ./ theta_pk[1][:,2]; // D/CL
+    // Concentration AUC at steady state
+    ss_auc[1] = (dose_ss ./ (theta_pk[1][:,2] .* theta_pk[1][:,3])); // D/(CL*V2)
 
     // Set AUC corresponding to each interval
     for(n in 1:N_int){
diff --git a/man/JointData.Rd b/man/JointData.Rd
index f191aaf..af30df3 100644
--- a/man/JointData.Rd
+++ b/man/JointData.Rd
@@ -78,7 +78,7 @@ Print info
 \subsection{Method \code{plot_dosing()}}{
 Plot dosing
 \subsection{Usage}{
-\if{html}{\out{<div class="r">}}\preformatted{JointData$plot_dosing(df_fit = NULL, max_num_subjects = 12)}\if{html}{\out{</div>}}
+\if{html}{\out{<div class="r">}}\preformatted{JointData$plot_dosing(df_fit = NULL, max_num_subjects = 12, subject_ids = NULL)}\if{html}{\out{</div>}}
 }
 
 \subsection{Arguments}{
@@ -87,6 +87,8 @@ Plot dosing
 \item{\code{df_fit}}{Fit data frame.}
 
 \item{\code{max_num_subjects}}{Max number of subjects to plot.}
+
+\item{\code{subject_ids}}{Which subjects to plot?}
 }
 \if{html}{\out{</div>}}
 }
diff --git a/man/MultistateModelFit.Rd b/man/MultistateModelFit.Rd
index b574e52..0a6e799 100644
--- a/man/MultistateModelFit.Rd
+++ b/man/MultistateModelFit.Rd
@@ -242,7 +242,8 @@ Plot PK fit.
   L = 100,
   timescale = 24,
   n_prev = 3,
-  ci_alpha = 0.9
+  ci_alpha = 0.9,
+  subject_ids = NULL
 )}\if{html}{\out{</div>}}
 }
 
@@ -263,6 +264,8 @@ training data is used.}
 \item{\code{n_prev}}{number of previous doses to show fit for}
 
 \item{\code{ci_alpha}}{Width of central credible interval.}
+
+\item{\code{subject_ids}}{Which subjects to plot?}
 }
 \if{html}{\out{</div>}}
 }
diff --git a/man/PKModel.Rd b/man/PKModel.Rd
index e035db1..a8a270c 100644
--- a/man/PKModel.Rd
+++ b/man/PKModel.Rd
@@ -154,7 +154,8 @@ concentration in the central compartment.
 \if{html}{\out{<a id="method-PKModel-compute_ss_auc"></a>}}
 \if{latex}{\out{\hypertarget{method-PKModel-compute_ss_auc}{}}}
 \subsection{Method \code{compute_ss_auc()}}{
-Compute steady-state area under curve over one dosing interval
+Compute steady-state area under concentration curve over
+one dosing interval
 \subsection{Usage}{
 \if{html}{\out{<div class="r">}}\preformatted{PKModel$compute_ss_auc(theta, dose)}\if{html}{\out{</div>}}
 }
diff --git a/man/PSSDosingData.Rd b/man/PSSDosingData.Rd
index 4811553..6ff85ea 100644
--- a/man/PSSDosingData.Rd
+++ b/man/PSSDosingData.Rd
@@ -116,7 +116,12 @@ a \code{data.frame}
 \subsection{Method \code{plot()}}{
 Plot dosing (and PK) data
 \subsection{Usage}{
-\if{html}{\out{<div class="r">}}\preformatted{PSSDosingData$plot(df_fit = NULL, subject_df = NULL, max_num_subjects = 12)}\if{html}{\out{</div>}}
+\if{html}{\out{<div class="r">}}\preformatted{PSSDosingData$plot(
+  df_fit = NULL,
+  subject_df = NULL,
+  max_num_subjects = 12,
+  subject_ids = NULL
+)}\if{html}{\out{</div>}}
 }
 
 \subsection{Arguments}{
@@ -128,6 +133,8 @@ and \code{upper}.}
 \item{\code{subject_df}}{Subject data frame.}
 
 \item{\code{max_num_subjects}}{Max number of subjects to plot.}
+
+\item{\code{subject_ids}}{Which subjects to plot?}
 }
 \if{html}{\out{</div>}}
 }
diff --git a/vignettes/math.Rmd b/vignettes/math.Rmd
index 3491b56..bf3cdb2 100644
--- a/vignettes/math.Rmd
+++ b/vignettes/math.Rmd
@@ -300,8 +300,7 @@ print(mod)
 ```
 
 Now we can simulate death events with known constant hazard rate `lambda`. For this
-model, we are simulating from a homogeneous Poisson process with rate $\lambda$,
-meaning that the event times should follow an exponential distribution with
+model, the event times should follow an exponential distribution with
 mean $\frac{1}{\lambda}$.
 
 ```{r, fig.width=5}
@@ -360,10 +359,11 @@ mod <- create_msm(tm, t_max = tmax)
 ```
 
 This model has 3 possible transitions. 
-Now we can simulate events with known constant hazard rates `lambda_1`, 
-`lambda_2`, `lambda_3`. For this
-model, we are simulating from independent homogeneous Poisson processes with 
-rates $\lambda_i$, meaning that the transition probabilities should be
+Now we can simulate transition times with known constant hazard rates `lambda_1`, 
+`lambda_2`, `lambda_3`. Simulating the transition times is mathematically equivalent to
+simulating from independent homogeneous Poisson processes with 
+rates $\lambda_i$, and taking the transition that occurred first. This means that 
+the transition probabilities should be
 $$
 P(i) = \frac{\lambda_i}{\sum_{j=1}^ 3 \lambda_j}.
 $$