diff --git a/CITATION.cff b/CITATION.cff
index 20b9a6b7e..6222c2131 100644
--- a/CITATION.cff
+++ b/CITATION.cff
@@ -8,7 +8,7 @@ message: 'To cite package "multinma" in publications use:'
 type: software
 license: GPL-3.0-only
 title: 'multinma: Bayesian Network Meta-Analysis of Individual and Aggregate Data'
-version: 0.6.0
+version: 0.6.1
 doi: 10.5281/zenodo.3904454
 abstract: Network meta-analysis and network meta-regression models for aggregate data,
   individual patient data, and mixtures of both individual and aggregate data using
@@ -28,7 +28,7 @@ preferred-citation:
     email: david.phillippo@bristol.ac.uk
     orcid: https://orcid.org/0000-0003-2672-7841
   year: '2023'
-  notes: R package version 0.6.0
+  notes: R package version 0.6.1
   url: https://dmphillippo.github.io/multinma/
   doi: 10.5281/zenodo.3904454
 repository: https://CRAN.R-project.org/package=multinma
diff --git a/DESCRIPTION b/DESCRIPTION
index 0a0c2eeb9..f6441312d 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: multinma
 Title: Bayesian Network Meta-Analysis of Individual and Aggregate Data
-Version: 0.6.0
+Version: 0.6.1
 Authors@R: 
     person(given = c("David", "M."),
            family = "Phillippo",
diff --git a/NEWS.md b/NEWS.md
index 1e80270b7..dbf9fb614 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,8 @@
+# multinma 0.6.1
+
+* Fix: Piecewise exponential hazard models no longer give errors during set-up. 
+Calculation of RW1 prior weights needed to be handled as a special case.
+
 # multinma 0.6.0
 
 ## Feature: Survival/time-to-event models are now supported
diff --git a/R/integration.R b/R/integration.R
index 8006f5563..ede3ebf4a 100644
--- a/R/integration.R
+++ b/R/integration.R
@@ -515,14 +515,19 @@ distr <- function(qfun, ...) {
 
 #' The Bernoulli Distribution
 #'
-#' The quantile function `qbern` for a Bernoulli distribution, with success
-#' probability `prob`. This is equivalent to `qbinom(p, 1, prob)`.
+#' The density function `dbern()`, distribution function `pbern()`, and quantile
+#' function `qbern()` for a Bernoulli distribution, with success probability
+#' `prob`. These are equivalent to `dbinom(p, 1, prob)`, `pbinom(p, 1, prob)`
+#' and `qbinom(p, 1, prob)`.
 #'
 #' @param x,q vector of quantiles
 #' @param p vector of probabilities
 #' @param prob probability of success
 #' @param lower.tail,log.p,log see [stats::Binomial]
 #'
+#' @return Numeric vector of length equal to the maximum of the lengths of the
+#'   input arguments.
+#'
 #' @export
 #' @rdname Bernoulli
 #' @aliases qbern
@@ -556,6 +561,9 @@ dbern <- function(x, prob, log = FALSE) {
 #' @param mean,sd mean and standard deviation, overriding `shape` and
 #'   `rate` or `scale` if specified
 #'
+#' @return Numeric vector of length equal to the maximum of the lengths of the
+#'   input arguments.
+#'
 #' @export
 #' @rdname GammaDist
 #' @aliases qgamma
@@ -606,6 +614,9 @@ pgamma <- function(q, shape, rate = 1, scale = 1/rate, lower.tail = TRUE,
 #' @param mean,sd mean and standard deviation, overriding `mu` and `sigma` if
 #'   specified
 #'
+#' @return Numeric vector of length equal to the maximum of the lengths of the
+#'   input arguments.
+#'
 #' @export
 #' @rdname logitNormal
 #' @aliases qlogitnorm
diff --git a/R/mspline.R b/R/mspline.R
index b04dede6d..9e2202dd4 100644
--- a/R/mspline.R
+++ b/R/mspline.R
@@ -268,8 +268,12 @@ rw1_prior_weights <- function(basis) {
   iknots <- attr(basis, "knots")
   bknots <- attr(basis, "Boundary.knots")
   knots <- c(rep(bknots[1], ord), iknots, rep(bknots[2], ord))
-  wts <- 1/(ord - 1) * (knots[(ord + 1):(nscoef + ord - 1)] - knots[2:nscoef])
-  return(sqrt(wts / sum(wts)))
+  if (ord == 1) {
+    wts <- (knots[2:nscoef] - knots[1:(nscoef - 1)]) / (knots[nscoef] - bknots[1])
+  } else {
+    wts <- (knots[(ord + 1):(nscoef + ord - 1)] - knots[2:nscoef]) / ((ord - 1) * (bknots[2] - bknots[1]))
+  }
+  return(sqrt(wts))
 }
 
 #' softmax transform
diff --git a/R/nma_data-class.R b/R/nma_data-class.R
index 691511469..fb29f50a0 100644
--- a/R/nma_data-class.R
+++ b/R/nma_data-class.R
@@ -66,6 +66,8 @@ NULL
 #' @param ... other options (not used)
 #' @param n number of studies of each type to print
 #'
+#' @return `x` is returned invisibly.
+#'
 #' @export
 print.nma_data <- function(x, ..., n = 10) {
   cwidth <- getOption("width")
diff --git a/R/nma_nodesplit-class.R b/R/nma_nodesplit-class.R
index b5137377a..f6b43c605 100644
--- a/R/nma_nodesplit-class.R
+++ b/R/nma_nodesplit-class.R
@@ -33,6 +33,8 @@ NULL
 #' @seealso The summary method [summary.nma_nodesplit_df()] summarises the
 #'   node-splitting results.
 #'
+#' @return `x` is returned invisibly.
+#'
 #' @export
 print.nma_nodesplit_df <- function(x, ...) {
   n_ns <- nrow(dplyr::filter(x, !is.na(.data$trt1) & !is.na(.data$trt2)))
diff --git a/R/predict.R b/R/predict.R
index a049990e7..c8fd19b7e 100644
--- a/R/predict.R
+++ b/R/predict.R
@@ -2089,12 +2089,11 @@ make_surv_predict <- function(eta, aux, times, likelihood,
   if (type %in% c("survival", "hazard", "cumhaz") && n_eta == 1) { # Multiple times, single linear predictor
     if (!is.null(aux)) aux <- matrix(aux, ncol = dim(aux)[3], dimnames = list(NULL, dimnames(aux)[[3]]))
     for (i in 1:length(times)) {
-      out[, , i] <- do.call(surv_predfuns[[likelihood]][[type]],
-                            args = list(times = times[i],
-                                        eta = as.vector(eta),
-                                        aux = aux,
-                                        quantiles = quantiles,
-                                        basis = basis))
+      out[, , i] <- surv_predfuns[[likelihood]][[type]](times = times[i],
+                                                        eta = as.vector(eta),
+                                                        aux = aux,
+                                                        quantiles = quantiles,
+                                                        basis = basis)
     }
   } else if (n_eta > 1) { # Single/multiple times, multiple linear predictors
     if (type == "quantile") {
@@ -2118,22 +2117,20 @@ make_surv_predict <- function(eta, aux, times, likelihood,
         }
       }
       out[ , , ((i-1)*iinc+1):(i*iinc)] <-
-        do.call(surv_predfuns[[likelihood]][[type]],
-                args = list(times = ti,
-                            eta = as.vector(eta[ , , i]),
-                            aux = auxi,
-                            quantiles = quantiles,
-                            basis = basis))
+        surv_predfuns[[likelihood]][[type]](times = ti,
+                                            eta = as.vector(eta[ , , i]),
+                                            aux = auxi,
+                                            quantiles = quantiles,
+                                            basis = basis)
     }
   } else { # Single time, single linear predictor
     if (!is.null(aux)) aux <- matrix(aux, ncol = dim(aux)[3], dimnames = list(NULL, dimnames(aux)[[3]]))
     if (type == "quantile") quantiles <- rep(quantiles, each = length(eta))
-    out <- array(do.call(surv_predfuns[[likelihood]][[type]],
-                         args = list(times = times,
-                                     eta = as.vector(eta),
-                                     aux = aux,
-                                     quantiles = quantiles,
-                                     basis = basis)),
+    out <- array(surv_predfuns[[likelihood]][[type]](times = times,
+                                                     eta = as.vector(eta),
+                                                     aux = aux,
+                                                     quantiles = quantiles,
+                                                     basis = basis),
                  dim = d_out, dimnames = dn_out)
   }
 
@@ -2284,11 +2281,10 @@ rmst_Sbar <- function(times, eta, weights, aux, likelihood, basis, start = 0) {
 pSbar <- function(times, eta, weights, aux, likelihood, basis) {
   Sbar <- numeric(length(times))
   for (i in 1:length(times)) {
-    S <- do.call(surv_predfuns[[likelihood]][["survival"]],
-                       args = list(times = times[i],
-                                   eta = eta,
-                                   aux = aux,
-                                   basis = basis))
+    S <- surv_predfuns[[likelihood]][["survival"]](times = times[i],
+                                                   eta = eta,
+                                                   aux = aux,
+                                                   basis = basis)
 
     Sbar[i] <- weighted.mean(S, weights)
   }
@@ -2346,11 +2342,10 @@ qSbar <- function(times, p, eta, weights, aux, likelihood, basis) {
         auxi <- as.vector(aux[ , , i])
       }
     }
-    S[ , , i] <- do.call(surv_predfuns[[likelihood]][["survival"]],
-                         args = list(times = times,
-                                     eta = as.vector(eta[ , , i]),
-                                     aux = auxi,
-                                     basis = basis))
+    S[ , , i] <- surv_predfuns[[likelihood]][["survival"]](times = times,
+                                                           eta = as.vector(eta[ , , i]),
+                                                           aux = auxi,
+                                                           basis = basis)
   }
 
   Sbar <- apply(S, MARGIN = 1:2, FUN = weighted.mean, w = weights)
diff --git a/R/stan_nma-class.R b/R/stan_nma-class.R
index 630d366d3..95bfc3465 100644
--- a/R/stan_nma-class.R
+++ b/R/stan_nma-class.R
@@ -40,6 +40,8 @@ NULL
 #' @param x A [stan_nma] object
 #' @param ... Further arguments passed to [print.stanfit()]
 #'
+#' @return `x` is returned invisibly.
+#'
 #' @export
 print.stan_nma <- function(x, ...) {
   if (inherits(x$network, "mlnmr_data")) type <- "ML-NMR"
diff --git a/README.md b/README.md
index 096b55d9a..7723194af 100644
--- a/README.md
+++ b/README.md
@@ -78,7 +78,7 @@ papers:
 The `multinma` package can be cited as follows:
 
 > Phillippo, D. M. (2024). *multinma: Bayesian Network Meta-Analysis of
-> Individual and Aggregate Data*. R package version 0.6.0, doi:
+> Individual and Aggregate Data*. R package version 0.6.1, doi:
 > [10.5281/zenodo.3904454](https://doi.org/10.5281/zenodo.3904454).
 
 When fitting ML-NMR models, please cite the methods paper:
diff --git a/cran-comments.md b/cran-comments.md
index 789429581..a1e1d7e5b 100644
--- a/cran-comments.md
+++ b/cran-comments.md
@@ -1,12 +1,35 @@
-This release ensures compatibility with the forthcoming update to rstan and
-StanHeaders, as well as providing extensive new features, improvements, and
-fixes.
+This release addresses an "additional checks" UBSAN failure, caused by a memory
+allocation bug in Stan (https://github.com/stan-dev/rstan/issues/1111).
+StanHeaders has now been patched which resolves the issue; no changes are
+necessary in multinma. Verified passing UBSAN checks with rocker/r-devel-san.
+
+In response to CRAN request, I have also:
+
+* Added the few missing \value tags to documentation
+
+Regarding the other issues raised:
+
+* \dontrun{} is used in man/pairs.stan_nma.Rd, as the example is designed to
+  raise warning messages. These are divergent transitions in Stan, which we then
+  demonstrate for the user how to examine using the pairs() plot.
+  
+* The global environment is never modified. rm(list = ls()) in some large test
+  files is never run on CRAN or by the user, these lines only run when testing 
+  on GitHub Actions to avoid running out of memory. These tests are skipped on 
+  CRAN to keep check time manageable.
+  
+* Similarly, .GlobalEnv is never modified. In the .Rd-files starting with 
+  aa_example_, the assign() to .GlobalEnv only runs and is only needed when 
+  building the package website using pkgdown, i.e. when pkgdown::in_pkgdown() is
+  TRUE. This code is never run otherwise.
+  
 
 ## Test environments
-* local R installation (Windows 10), R 4.3.2
+* local rocker/r-devel-san image with sanitizers
+* local R installation (Windows 10), R 4.3.3
 * win-builder (release, devel)
-* Ubuntu 22.04.3 on GitHub Actions (release, devel, oldrel)
-* Mac OS 12.7.2 on GitHub Actions (release)
+* Ubuntu 22.04.4 on GitHub Actions (release, devel, oldrel)
+* Mac OS 12.7.3 on GitHub Actions (release)
 * Windows Server 2022 10.0.20348 on GitHub Actions (release)
 
 ## R CMD check results
diff --git a/man/Bernoulli.Rd b/man/Bernoulli.Rd
index acc87885c..b3389a6db 100644
--- a/man/Bernoulli.Rd
+++ b/man/Bernoulli.Rd
@@ -21,7 +21,13 @@ dbern(x, prob, log = FALSE)
 
 \item{x, q}{vector of quantiles}
 }
+\value{
+Numeric vector of length equal to the maximum of the lengths of the
+input arguments.
+}
 \description{
-The quantile function \code{qbern} for a Bernoulli distribution, with success
-probability \code{prob}. This is equivalent to \code{qbinom(p, 1, prob)}.
+The density function \code{dbern()}, distribution function \code{pbern()}, and quantile
+function \code{qbern()} for a Bernoulli distribution, with success probability
+\code{prob}. These are equivalent to \code{dbinom(p, 1, prob)}, \code{pbinom(p, 1, prob)}
+and \code{qbinom(p, 1, prob)}.
 }
diff --git a/man/GammaDist.Rd b/man/GammaDist.Rd
index 6f5e69630..d2ff4164e 100644
--- a/man/GammaDist.Rd
+++ b/man/GammaDist.Rd
@@ -40,6 +40,10 @@ pgamma(
 
 \item{x, q}{vector of quantiles}
 }
+\value{
+Numeric vector of length equal to the maximum of the lengths of the
+input arguments.
+}
 \description{
 We provide convenient extensions of the \verb{[dpq]gamma} functions, which allow
 the distribution to be specified in terms of its mean and standard deviation,
diff --git a/man/logitNormal.Rd b/man/logitNormal.Rd
index a7f16b12a..3a6309924 100644
--- a/man/logitNormal.Rd
+++ b/man/logitNormal.Rd
@@ -22,6 +22,10 @@ specified}
 
 \item{q}{vector of probabilities}
 }
+\value{
+Numeric vector of length equal to the maximum of the lengths of the
+input arguments.
+}
 \description{
 We provide convenient extensions of the \verb{[dpq]logitnorm} functions in the
 package \link[logitnorm:logitnorm-package]{logitnorm}, which allow the
diff --git a/man/print.nma_data.Rd b/man/print.nma_data.Rd
index 71110cf00..a67a7857f 100644
--- a/man/print.nma_data.Rd
+++ b/man/print.nma_data.Rd
@@ -16,6 +16,9 @@
 
 \item{n}{number of studies of each type to print}
 }
+\value{
+\code{x} is returned invisibly.
+}
 \description{
 Print details of networks stored as \link{nma_data} objects, as created by
 \code{\link[=set_ipd]{set_ipd()}}, \code{\link[=set_agd_arm]{set_agd_arm()}}, \code{\link[=set_agd_contrast]{set_agd_contrast()}}, \code{\link[=set_agd_surv]{set_agd_surv()}}, or
diff --git a/man/print.nma_nodesplit_df.Rd b/man/print.nma_nodesplit_df.Rd
index b362f8b4c..d9b02d730 100644
--- a/man/print.nma_nodesplit_df.Rd
+++ b/man/print.nma_nodesplit_df.Rd
@@ -14,6 +14,9 @@
 
 \item{...}{Further arguments passed to \code{\link[=print.stanfit]{print.stanfit()}}}
 }
+\value{
+\code{x} is returned invisibly.
+}
 \description{
 Print \code{nma_nodesplit_df} objects
 }
diff --git a/man/print.stan_nma.Rd b/man/print.stan_nma.Rd
index f63332147..c71255c9c 100644
--- a/man/print.stan_nma.Rd
+++ b/man/print.stan_nma.Rd
@@ -11,6 +11,9 @@
 
 \item{...}{Further arguments passed to \code{\link[=print.stanfit]{print.stanfit()}}}
 }
+\value{
+\code{x} is returned invisibly.
+}
 \description{
 Print \code{stan_nma} objects
 }