diff --git a/article/article.Rmd b/article/article.Rmd
index 1d7a1757..591f03b9 100644
--- a/article/article.Rmd
+++ b/article/article.Rmd
@@ -803,7 +803,7 @@ In principle all models provide an estimate for "no-effect" toxicity concentrati
 averaging approach, can yield robust estimates of N(S)EC and of their uncertainty within a single
 analysis framework [@fisher2023ieam]. Both NEC and NSEC can be calculated from fitted models using the functions \code{nec} and \code{nsec}.
 
-ECx estimates can be equally obtained from both \code{"nec"} and \code{"ecx"} models. ECx estimates will usually be lower (more conservative) for \code{"ecx"} models fitted to the same data as \code{"nec"} models. There is ambiguity in the definition of ECx estimates from hormesis models---these allow an initial increase in the response [see @Mattson2008] and include models with the string **horme** in their name---as well as those that have no natural lower bound on the scale of the response (models with the string **lin** in their name, in the case of Gaussian response data). For this reason the \code{ecx} function has arguments \code{hormesis_def} and \code{type}, both character vectors indicating the desired behaviour. For \code{hormesis_def = "max"}, ECx values are calculated as a decline from the maximum estimates (i.e., the peak at $\eta = \text{NEC}$); and \code{hormesis_def = "control"} (the default) indicates that ECx values should be calculated relative to the control, which is assumed to be the lowest observed concentration. For \code{type = "relative"} ECx is calculated as the percentage decrease from the maximum predicted value of the response ($\tau = \text{top}$) to the minimum predicted value of the response (i.e., *relative* to the observed data). For \code{type = "absolute"} (the default) ECx is calculated as the percentage decrease from the maximum value of the response ($\tau = \text{top}$) to 0 (or $\delta = \text{bottom}$ for models with that parameter). For \code{type = "direct"}, a direct interpolation of the response on the predictor is obtained.
+ECx estimates can be equally obtained from both \code{"nec"} and \code{"ecx"} models. ECx estimates will usually be lower (more conservative) for \code{"ecx"} models fitted to the same data as \code{"nec"} models. There is ambiguity in the definition of ECx estimates from hormesis models---these allow an initial increase in the response [see @Mattson2008] and include models with the string **horme** in their name---as well as those that have no natural lower bound on the scale of the response (models with the string **lin** in their name, in the case of Gaussian response data). For this reason the \code{ecx} function has arguments \code{hormesis_def} and \code{type}, both character vectors indicating the desired behaviour. For \code{hormesis_def = "max"}, ECx values are calculated as a decline from the maximum estimates (i.e., the peak at $\eta = \text{NEC}$); and \code{hormesis_def = "control"} (the default) indicates that ECx values should be calculated relative to the control, which is assumed to be the lowest observed concentration. For \code{type = "relative"} ECx is calculated as the percentage decrease from the maximum predicted value of the response ($\tau = \text{top}$) to the minimum predicted value of the response (i.e., *relative* to the observed data). For \code{type = "absolute"} (the default) ECx is calculated as the percentage decrease from the maximum value of the response ($\tau = \text{top}$) to 0. For \code{type = "direct"}, a direct interpolation of the response on the predictor is obtained.
 
 ## Model suitability for response types
 
@@ -819,6 +819,8 @@ Considerations of analytical reproducibility are particularly relevant to CR mod
 
 To help with reproducibility \pkg{bayesnec} now allows a seed to be passed to \pkg{brm} and \proglang{Stan}. If a seed is used in the \code{bnec} call, it will also be used internally to generate initial values. Although in \proglang{R} seeds are consistent across versions and operational systems, and therefore the initial values will be the same across different users for a given seed, the underlying \proglang{Stan} model fitting mechanism may yield slightly different parameter estimates for known reasons relating to floating point operations [see chapter 20 in @stan2021]. A potentially better strategy for ensuring reproducibility is to build a docker (https://docs.docker.com/get-docker/) container, an approach representing one strategy towards overcoming the reproducibility crisis [@Baker2016]. Also note that while setting a seed can be useful to obtain consistent outputs it might be worth examining how robust the inference is across different seeds.
 
+Due to the fact that the underlying \code{brmsfit} model fitted using \pkg{cmdstanr} does not retain initial values as part of the returned model object, reproducibility may be reduced when using \pkg{cmdstanr}.
+
 ## Computational details {short-title="Computational detail" #compdetails}
 
 All computations in this paper were performed using rmarkdown [@allaire2023] with the following
diff --git a/vignettes/example2b.Rmd b/vignettes/example2b.Rmd
index e23d1ff0..f17e842d 100644
--- a/vignettes/example2b.Rmd
+++ b/vignettes/example2b.Rmd
@@ -102,7 +102,7 @@ All models provide an estimate for the No-Effect-Concentration (*NEC*). For mode
 
 *EC~x~* estimates can be equally obtained from both `"nec"` and `"ecx"` models. *EC~x~* estimates will usually be lower (more conservative) for `"ecx"` models fitted to the same data as `"nec"` models (see the [Comparing posterior predictions][e4]) vignette for an example. However, we recommend using `"all"` models where *EC~x~* estimation is required because `"nec"` models can fit some datasets better than `"ecx"` models and the model averaging approach will place the greatest weight for the outcome that best fits the supplied data. This approach will yield *EC~x~* estimates that are the most representative of the underlying relationship in the dataset.
 
-There is ambiguity in the definition of *EC~x~* estimates from hormesis models---these allow an initial increase in the response [see @Mattson2008] and include models with the character string `horme` in their name---as well as those that have no natural lower bound on the scale of the response (models with the string **lin** in their name, in the case of Gaussian response data). For this reason the `ecx` function has arguments `hormesis_def` and `type`, both character vectors indicating the desired behaviour. For `hormesis_def = "max"`, *EC~x~* values are calculated as a decline from the maximum estimates (i.e., the peak at $\eta = \text{NEC}$); and `hormesis_def = "control"` (the default) indicates that *EC~x~* values should be calculated relative to the control, which is assumed to be the lowest observed concentration. For `type = "relative"` *EC~x~* is calculated as the percentage decrease from the maximum predicted value of the response ($\tau = \text{top}$) to the minimum predicted value of the response (i.e., *relative* to the observed data). For `type = "absolute"` (the default) *EC~x~* is calculated as the percentage decrease from the maximum value of the response ($\tau = \text{top}$) to 0 (or $\delta = \text{bottom}$ for models with that parameter). For `type = "direct"`, a direct interpolation of the response on the predictor is obtained.
+There is ambiguity in the definition of *EC~x~* estimates from hormesis models---these allow an initial increase in the response [see @Mattson2008] and include models with the character string `horme` in their name---as well as those that have no natural lower bound on the scale of the response (models with the string **lin** in their name, in the case of Gaussian response data). For this reason the `ecx` function has arguments `hormesis_def` and `type`, both character vectors indicating the desired behaviour. For `hormesis_def = "max"`, *EC~x~* values are calculated as a decline from the maximum estimates (i.e., the peak at $\eta = \text{NEC}$); and `hormesis_def = "control"` (the default) indicates that *EC~x~* values should be calculated relative to the control, which is assumed to be the lowest observed concentration. For `type = "relative"` *EC~x~* is calculated as the percentage decrease from the maximum predicted value of the response ($\tau = \text{top}$) to the minimum predicted value of the response (i.e., *relative* to the observed data). For `type = "absolute"` (the default) *EC~x~* is calculated as the percentage decrease from the maximum value of the response ($\tau = \text{top}$). For `type = "direct"`, a direct interpolation of the response on the predictor is obtained.
 
 # Model suitability for response types
 
diff --git a/vignettes/example2b.Rmd.orig b/vignettes/example2b.Rmd.orig
index 92c38934..8f61a26d 100644
--- a/vignettes/example2b.Rmd.orig
+++ b/vignettes/example2b.Rmd.orig
@@ -81,7 +81,7 @@ All models provide an estimate for the No-Effect-Concentration (*NEC*). For mode
 
 *EC~x~* estimates can be equally obtained from both `"nec"` and `"ecx"` models. *EC~x~* estimates will usually be lower (more conservative) for `"ecx"` models fitted to the same data as `"nec"` models (see the [Comparing posterior predictions][e4]) vignette for an example. However, we recommend using `"all"` models where *EC~x~* estimation is required because `"nec"` models can fit some datasets better than `"ecx"` models and the model averaging approach will place the greatest weight for the outcome that best fits the supplied data. This approach will yield *EC~x~* estimates that are the most representative of the underlying relationship in the dataset.
 
-There is ambiguity in the definition of *EC~x~* estimates from hormesis models---these allow an initial increase in the response [see @Mattson2008] and include models with the character string `horme` in their name---as well as those that have no natural lower bound on the scale of the response (models with the string **lin** in their name, in the case of Gaussian response data). For this reason the `ecx` function has arguments `hormesis_def` and `type`, both character vectors indicating the desired behaviour. For `hormesis_def = "max"`, *EC~x~* values are calculated as a decline from the maximum estimates (i.e., the peak at $\eta = \text{NEC}$); and `hormesis_def = "control"` (the default) indicates that *EC~x~* values should be calculated relative to the control, which is assumed to be the lowest observed concentration. For `type = "relative"` *EC~x~* is calculated as the percentage decrease from the maximum predicted value of the response ($\tau = \text{top}$) to the minimum predicted value of the response (i.e., *relative* to the observed data). For `type = "absolute"` (the default) *EC~x~* is calculated as the percentage decrease from the maximum value of the response ($\tau = \text{top}$) to 0 (or $\delta = \text{bottom}$ for models with that parameter). For `type = "direct"`, a direct interpolation of the response on the predictor is obtained.
+There is ambiguity in the definition of *EC~x~* estimates from hormesis models---these allow an initial increase in the response [see @Mattson2008] and include models with the character string `horme` in their name---as well as those that have no natural lower bound on the scale of the response (models with the string **lin** in their name, in the case of Gaussian response data). For this reason the `ecx` function has arguments `hormesis_def` and `type`, both character vectors indicating the desired behaviour. For `hormesis_def = "max"`, *EC~x~* values are calculated as a decline from the maximum estimates (i.e., the peak at $\eta = \text{NEC}$); and `hormesis_def = "control"` (the default) indicates that *EC~x~* values should be calculated relative to the control, which is assumed to be the lowest observed concentration. For `type = "relative"` *EC~x~* is calculated as the percentage decrease from the maximum predicted value of the response ($\tau = \text{top}$) to the minimum predicted value of the response (i.e., *relative* to the observed data). For `type = "absolute"` (the default) *EC~x~* is calculated as the percentage decrease from the maximum value of the response ($\tau = \text{top}$) to 0. For `type = "direct"`, a direct interpolation of the response on the predictor is obtained.
 
 # Model suitability for response types