diff --git a/DESCRIPTION b/DESCRIPTION
index 53c0dad..8a70858 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -21,7 +21,8 @@ Roxygen: list(markdown = TRUE)
 RoxygenNote: 7.2.3
 Imports: 
     ggplot2 (>= 2.2.1),
-    stats
+    stats,
+    withr
 Suggests: 
     testthat (>= 3.0.0),
     nlme,
diff --git a/R/circa_single.R b/R/circa_single.R
index 0c01cb0..e83d603 100644
--- a/R/circa_single.R
+++ b/R/circa_single.R
@@ -151,7 +151,7 @@ circa_single <- function(x,
       )
     if (data_rhythmic) {
       fig_out <- p +
-        ggplot2::stat_function(fun = eq, size = 1) +
+        ggplot2::stat_function(fun = eq, linewidth = 1) +
         ggplot2::labs(subtitle = "Data is rhythmic", x = "time (hours)")
     } else {
       fig_out <- p +
diff --git a/R/circa_single_mixed.R b/R/circa_single_mixed.R
index 96c0dad..c713b4e 100644
--- a/R/circa_single_mixed.R
+++ b/R/circa_single_mixed.R
@@ -23,6 +23,7 @@
 #' @export
 #'
 #' @examples
+#' \donttest{
 #' set.seed(42)
 #' mixed_data <- function(n) {
 #'   counter <- 1
@@ -50,7 +51,7 @@
 #'   x = df, col_time = "time", col_outcome = "measure",
 #'   col_id = "id", randomeffects = c("k"), weights = sw
 #' )
-#'
+#' }
 circa_single_mixed <- function(x,
                                col_time,
                                col_outcome,
@@ -190,7 +191,7 @@ circa_single_mixed <- function(x,
   if (return_figure) {
     if (data_rhythmic) {
       fig_out <- ggplot2::ggplot(x, ggplot2::aes(time, measure)) +
-        ggplot2::stat_function(fun = eq, size = 1) +
+        ggplot2::stat_function(fun = eq, linewidth = 1) +
         ggplot2::geom_point() +
         ggplot2::xlim(
           min(floor(x$time / period) * period),
diff --git a/R/circacompare.R b/R/circacompare.R
index 6aa79c4..c1d4ba9 100644
--- a/R/circacompare.R
+++ b/R/circacompare.R
@@ -199,8 +199,8 @@ circacompare <- function(x,
   eval(parse(text = eq_expression$g2))
 
   fig_out <- ggplot2::ggplot(x, ggplot2::aes(time, measure)) +
-    ggplot2::stat_function(ggplot2::aes(colour = group_1_text), fun = eq_1, size = 1) +
-    ggplot2::stat_function(ggplot2::aes(colour = group_2_text), fun = eq_2, size = 1) +
+    ggplot2::stat_function(ggplot2::aes(colour = group_1_text), fun = eq_1, linewidth = 1) +
+    ggplot2::stat_function(ggplot2::aes(colour = group_2_text), fun = eq_2, linewidth = 1) +
     ggplot2::geom_point(ggplot2::aes(colour = group)) +
     ggplot2::scale_colour_manual(
       breaks = c(group_1_text, group_2_text),
diff --git a/R/circacompare_mixed.R b/R/circacompare_mixed.R
index 495f9f7..c63e64a 100644
--- a/R/circacompare_mixed.R
+++ b/R/circacompare_mixed.R
@@ -24,6 +24,7 @@
 #'
 #' @examples
 #' # Generate some data with within-id correlation for phase-shift (phi1)
+#' \donttest{
 #' set.seed(99)
 #' phi1_in <- 3.15
 #'
@@ -62,6 +63,7 @@
 #'   control = list(grouped_params = c("phi"), random_params = c("phi1")),
 #'   weights = sw
 #' )
+#' }
 #'
 circacompare_mixed <- function(x,
                                col_time,
@@ -265,8 +267,8 @@ circacompare_mixed <- function(x,
   eval(parse(text = eq_expression$g2))
 
   fig_out <- ggplot2::ggplot(x, ggplot2::aes(time, measure)) +
-    ggplot2::stat_function(ggplot2::aes(colour = group_1_text), fun = eq_1, size = 1) +
-    ggplot2::stat_function(ggplot2::aes(colour = group_2_text), fun = eq_2, size = 1) +
+    ggplot2::stat_function(ggplot2::aes(colour = group_1_text), fun = eq_1, linewidth = 1) +
+    ggplot2::stat_function(ggplot2::aes(colour = group_2_text), fun = eq_2, linewidth = 1) +
     ggplot2::geom_point(ggplot2::aes(colour = group)) +
     ggplot2::scale_colour_manual(
       breaks = c(group_1_text, group_2_text),
diff --git a/R/make_data.R b/R/make_data.R
index ba341fc..04a3fd6 100644
--- a/R/make_data.R
+++ b/R/make_data.R
@@ -20,7 +20,19 @@
 #' data <- make_data(k1 = 3, alpha1 = 4, phi1 = 6)
 make_data <- function(k = 0, k1 = 3, alpha = 10, alpha1 = 4, phi = 0, phi1 = 3.15, tau = 24, hours = 48, noise_sd = 0.1, seed = NULL) {
   if (!is.null(seed)) {
-    set.seed(seed)
+    withr::with_seed(seed, {
+      return(make_data(
+        k = k,
+        k1 = k1,
+        alpha = alpha,
+        alpha1 = alpha1,
+        phi = phi,
+        phi1 = phi1,
+        tau = tau,
+        hours = hours,
+        noise_sd = noise_sd
+      ))
+    })
   }
   g1 <- data.frame(
     time = rep(NA, hours),
diff --git a/man/circa_single_mixed.Rd b/man/circa_single_mixed.Rd
index c0eb397..b223d18 100644
--- a/man/circa_single_mixed.Rd
+++ b/man/circa_single_mixed.Rd
@@ -58,6 +58,7 @@ list
 \code{circa_single_mixed} is similar to \code{circa_single} but allows for some simple, user-specified random-effects on the rhythmic parameters of choice.
 }
 \examples{
+\donttest{
 set.seed(42)
 mixed_data <- function(n) {
   counter <- 1
@@ -85,5 +86,5 @@ out2 <- circa_single_mixed(
   x = df, col_time = "time", col_outcome = "measure",
   col_id = "id", randomeffects = c("k"), weights = sw
 )
-
+}
 }
diff --git a/man/circacompare_mixed.Rd b/man/circacompare_mixed.Rd
index 63e3e4f..549a035 100644
--- a/man/circacompare_mixed.Rd
+++ b/man/circacompare_mixed.Rd
@@ -59,6 +59,7 @@ list
 }
 \examples{
 # Generate some data with within-id correlation for phase-shift (phi1)
+\donttest{
 set.seed(99)
 phi1_in <- 3.15
 
@@ -97,5 +98,6 @@ out2 <- circacompare_mixed(
   control = list(grouped_params = c("phi"), random_params = c("phi1")),
   weights = sw
 )
+}
 
 }
diff --git a/tests/testthat/test-circa_single.R b/tests/testthat/test-circa_single.R
index 8a3d8af..86ca1a2 100644
--- a/tests/testthat/test-circa_single.R
+++ b/tests/testthat/test-circa_single.R
@@ -1,28 +1,26 @@
 test_that("circa_single works", {
-  set.seed(42)
   tau <- 15
-  data_rhythmic <- make_data(k1 = 0, alpha1 = 0, phi = pi, phi1 = 0, noise_sd = 1)
-  out_rhythmic <- circa_single(x = data_rhythmic, col_time = "time", col_outcome = "measure")
-
-  data_rhythmic$time <- data_rhythmic$time / 24 * tau
-  out_rhythmic_free_tau <-
-    circa_single(
-      x = data_rhythmic, col_time = "time", col_outcome = "measure",
-      period = NA,
-      control = list(
-        main_params = c("k", "alpha", "phi", "tau"),
-        period_min = tau - 5,
-        period_max = tau + 5
+  withr::with_seed(42, {
+
+    data_rhythmic <- make_data(k1 = 0, alpha1 = 0, phi = pi, phi1 = 0, noise_sd = 1)
+    out_rhythmic <- circa_single(x = data_rhythmic, col_time = "time", col_outcome = "measure")
+
+    data_rhythmic$time <- data_rhythmic$time / 24 * tau
+    out_rhythmic_free_tau <-
+      circa_single(
+        x = data_rhythmic, col_time = "time", col_outcome = "measure",
+        period = NA,
+        control = list(
+          main_params = c("k", "alpha", "phi", "tau"),
+          period_min = tau - 5,
+          period_max = tau + 5
+        )
       )
-    )
-
-
-  # out_rhythmic_free_tau$plot + ggplot2::geom_vline(xintercept=out_rhythmic_free_tau$summary[5, 'value'])
-
-  data_arrhythmic <- make_data(alpha = 0)
-  data_arrhythmic <- data_arrhythmic[data_arrhythmic$group == "g1", ]
-  out_arrhythmic <- circa_single(x = data_arrhythmic, col_time = "time", col_outcome = "measure")
 
+    data_arrhythmic <- make_data(alpha = 0)
+    data_arrhythmic <- data_arrhythmic[data_arrhythmic$group == "g1", ]
+    out_arrhythmic <- circa_single(x = data_arrhythmic, col_time = "time", col_outcome = "measure")
+  })
   expect_true(class(out_rhythmic) == "list") # no errors when running circa_single()
   expect_true(out_rhythmic$summary[1, 2] < 0.01) # amplitude_p for rhythmic data is small
   expect_true(out_arrhythmic$summary[1, 2] > 0.05) # amplitude_p for arrhythmic data is large.
@@ -37,22 +35,23 @@ test_that("circa_single works", {
   # assess whether decay on amplitude per-hour is modelled well when period is parameterized
   alpha_decay_in <- 0.01
   tau_in <- 16
-  df <- make_data(k = 5, k1 = 0, alpha = 20, alpha1 = 0, phi = 2, phi1 = 0, hours = 96, noise_sd = 1)
-  df$time <- df$time / 24 * tau_in
-  df$measure <- df$measure * exp(-alpha_decay_in * (df$time))
-  out_alpha_decay <- circa_single(
-    x = df,
-    col_time = "time",
-    col_outcome = "measure",
-    period = NA,
-    control = list(
-      main_params = c("k", "alpha", "phi", "tau"),
-      decay_params = c("alpha"),
-      period_min = 12,
-      period_max = 20
+  withr::with_seed(1, {
+    df <- make_data(k = 5, k1 = 0, alpha = 20, alpha1 = 0, phi = 2, phi1 = 0, hours = 96, noise_sd = 1)
+    df$time <- df$time / 24 * tau_in
+    df$measure <- df$measure * exp(-alpha_decay_in * (df$time))
+    out_alpha_decay <- circa_single(
+      x = df,
+      col_time = "time",
+      col_outcome = "measure",
+      period = NA,
+      control = list(
+        main_params = c("k", "alpha", "phi", "tau"),
+        decay_params = c("alpha"),
+        period_min = 12,
+        period_max = 20
+      )
     )
-  )
-  out_alpha_decay
+  })
 
   fit_alpha_decay <- extract_model_coefs(out_alpha_decay$fit)["alpha_decay", ]
   alpha_decay_est <- fit_alpha_decay["estimate"]
@@ -80,22 +79,22 @@ test_that("suppress_all works", {
   class = "data.frame", row.names = c(NA, -16L)
   )
 
-  set.seed(1)
-  output <- capture.output(
-    {
-      circa_single(x = y, col_time = "time", col_outcome = "value", return_figure = FALSE)
-    },
-    type = "message"
-  )
-  expect_true(length(output) > 1)
+  withr::with_seed(1, {
+    output <- capture.output(
+      circa_single(x = y, col_time = "time", col_outcome = "value", return_figure = FALSE),
+      type = "message"
+    )
 
-  output <- capture.output(
-    {
-      circa_single(x = y, col_time = "time", col_outcome = "value", return_figure = FALSE, suppress_all = TRUE)
-    },
-    type = "message"
-  )
-  expect_true(length(output) == 0)
+
+    output_suppressed <- capture.output(
+      circa_single(x = y, col_time = "time", col_outcome = "value", return_figure = FALSE, suppress_all = TRUE),
+      type = "message"
+    )
+
+  })
+
+  expect_true(length(output) > 1)
+  expect_true(length(output_suppressed) == 0)
 })
 
 ### make test that weights are used correctly and malformatted weights are detected
diff --git a/tests/testthat/test-circa_single_mixed.R b/tests/testthat/test-circa_single_mixed.R
index 1659051..5801432 100644
--- a/tests/testthat/test-circa_single_mixed.R
+++ b/tests/testthat/test-circa_single_mixed.R
@@ -1,5 +1,4 @@
 test_that("circa_single_mixed() works", {
-  set.seed(42)
   tau_in <- 14
   mixed_data <- function(n) {
     counter <- 1
@@ -15,24 +14,29 @@ test_that("circa_single_mixed() works", {
     }
     return(res)
   }
-  df <- mixed_data(n = 12)
+  withr::with_seed(42, {
+    df <- mixed_data(n = 12)
 
-  out <- circa_single_mixed(
-    x = df, col_time = "time", col_outcome = "measure",
-    col_id = "id", randomeffects = c("k"),
-    control = list(decay_params = c("k"))
-  )
+    out <- circa_single_mixed(
+      x = df, col_time = "time", col_outcome = "measure",
+      col_id = "id", randomeffects = c("k"),
+      control = list(decay_params = c("k"))
+    )
+  })
 
   expect_true(class(out) == "list")
   expect_true(round(extract_model_coefs(out$fit)["alpha", "estimate"]) == 10)
 
   df$time <- df$time / 24 * tau_in
-  out_tau_adjusted <-
-    circa_single_mixed(
-      x = df, col_time = "time", col_outcome = "measure",
-      col_id = "id", randomeffects = c("k"), period = NA,
-      control = list(period_param = T, period_min = tau_in - 4, period_max = tau_in + 4)
-    )
+
+  withr::with_seed(42, {
+    out_tau_adjusted <-
+      circa_single_mixed(
+        x = df, col_time = "time", col_outcome = "measure",
+        col_id = "id", randomeffects = c("k"), period = NA,
+        control = list(period_param = T, period_min = tau_in - 4, period_max = tau_in + 4)
+      )
+  })
 
   fit_tau <- extract_model_coefs(out_tau_adjusted$fit)["tau", ]
   tau_est <- fit_tau["estimate"]
@@ -43,7 +47,6 @@ test_that("circa_single_mixed() works", {
 
 ### make test that weights are used correctly and malformatted weights are detected
 test_that("weights work", {
-  set.seed(42)
   mixed_data <- function(n) {
     counter <- 1
     for (i in 1:n) {
@@ -59,49 +62,53 @@ test_that("weights work", {
     return(res)
   }
 
-  df <- mixed_data(n = 50)
-
-  # no weights used (= all weights are 1), hence fit$apVar should not be populated
-  out <- circa_single_mixed(
-    x = df, col_time = "time", col_outcome = "measure",
-    col_id = "id", randomeffects = c("k")
-  )
-  expect_true(is(out$fit$apVar, "character"))
-
-  # when weights are not all 1 then fit$apVar should be a matrix
-  sw <- runif(n = nrow(df))
-  out2 <- circa_single_mixed(
-    x = df, col_time = "time", col_outcome = "measure",
-    col_id = "id", randomeffects = c("k"), weights = sw
-  )
-  expect_true(is(out2$fit$apVar, "matrix"))
+  withr::with_seed(42, {
+    df <- mixed_data(n = 50)
 
-  # weights must be same length as nrow(x)
-  sw2 <- c(sw, 1)
-  expect_error(
-    circa_single_mixed(
+    out <- circa_single_mixed(
       x = df, col_time = "time", col_outcome = "measure",
-      col_id = "id", randomeffects = c("k"), weights = sw2
+      col_id = "id", randomeffects = c("k")
     )
-  )
 
-  # weights must not contain NA
-  sw3 <- sw
-  sw3[1] <- NA
-  expect_error(
-    circa_single_mixed(
+    sw <- runif(n = nrow(df))
+    out2 <- circa_single_mixed(
       x = df, col_time = "time", col_outcome = "measure",
-      col_id = "id", randomeffects = c("k"), weights = sw3
+      col_id = "id", randomeffects = c("k"), weights = sw
     )
-  )
 
-  # weights must not be negative
-  sw4 <- sw
-  sw4[1] <- -1
-  expect_error(
-    circa_single_mixed(
-      x = df, col_time = "time", col_outcome = "measure",
-      col_id = "id", randomeffects = c("k"), weights = sw4
+    # weights must be same length as nrow(x)
+    sw2 <- c(sw, 1)
+    expect_error(
+      circa_single_mixed(
+        x = df, col_time = "time", col_outcome = "measure",
+        col_id = "id", randomeffects = c("k"), weights = sw2
+      )
+    )
+
+    # weights must not contain NA
+    sw3 <- sw
+    sw3[1] <- NA
+    expect_error(
+      circa_single_mixed(
+        x = df, col_time = "time", col_outcome = "measure",
+        col_id = "id", randomeffects = c("k"), weights = sw3
+      )
     )
-  )
+
+    # weights must not be negative
+    sw4 <- sw
+    sw4[1] <- -1
+    expect_error(
+      circa_single_mixed(
+        x = df, col_time = "time", col_outcome = "measure",
+        col_id = "id", randomeffects = c("k"), weights = sw4
+      )
+    )
+  })
+
+  # no weights used (= all weights are 1), hence fit$apVar should not be populated
+  expect_true(is(out$fit$apVar, "character"))
+
+  # when weights are not all 1 then fit$apVar should be a matrix
+  expect_true(is(out2$fit$apVar, "matrix"))
 })
diff --git a/tests/testthat/test-circacompare.R b/tests/testthat/test-circacompare.R
index 13a2001..ab30452 100644
--- a/tests/testthat/test-circacompare.R
+++ b/tests/testthat/test-circacompare.R
@@ -1,20 +1,22 @@
 test_that("circacompare() fits a good model to generated data", {
-  set.seed(42)
   tau_in <- 15
   phi1_in <- 12
-  df <- make_data(phi1 = (phi1_in / 24) * (2 * pi), noise_sd = 2)
-  out <- circacompare(x = df, col_time = "time", col_group = "group", col_outcome = "measure")
 
-  df$time <- df$time / 24 * tau_in
-  out_tau_adjusted <- circacompare(
-    x = df, col_time = "time", col_group = "group", col_outcome = "measure",
-    period = NA,
-    control = list(
-      main_params = c("k", "alpha", "tau", "phi"),
-      grouped_params = c("k", "alpha", "phi"),
-      period_min = tau_in - 4, period_max = tau_in + 4
+  withr::with_seed(42, {
+    df <- make_data(phi1 = (phi1_in / 24) * (2 * pi), noise_sd = 2)
+    out <- circacompare(x = df, col_time = "time", col_group = "group", col_outcome = "measure")
+
+    df$time <- df$time / 24 * tau_in
+    out_tau_adjusted <- circacompare(
+      x = df, col_time = "time", col_group = "group", col_outcome = "measure",
+      period = NA,
+      control = list(
+        main_params = c("k", "alpha", "tau", "phi"),
+        grouped_params = c("k", "alpha", "phi"),
+        period_min = tau_in - 4, period_max = tau_in + 4
+      )
     )
-  )
+  })
 
   both_groups_rhythmic <- as.logical(out$summary[1, "value"] < 0.05 & out$summary[2, "value"] < 0.05)
   phase_shift_estimated_within_2hours <- abs(abs(out$summary[13, "value"]) - phi1_in) < 2
@@ -31,22 +33,25 @@ test_that("circacompare() fits a good model to generated data", {
 
   # create some longer time period data and keep all parameters the same except amplitude
   # create some decay in one group for amplitude and test whether it's well estimated by the model.
-  df <- make_data(k1 = 0, alpha1 = 10, phi1 = 0, seed = 42, hours = 96, noise_sd = 2)
-  df$time <- df$time / 24 * tau_in
-  alpha_decay1_in <- 0.05
-  # note that decay is on a scale of time in radians, not time in hours.
-  df$measure[df$group == "g2"] <- df$measure[df$group == "g2"] * exp(-alpha_decay1_in * df$time[df$group == "g2"])
 
-  out_alpha_decay <-
-    circacompare(
-      x = df, "time", "group", "measure", period = NA,
-      control = list(
-        main_params = c("k", "alpha", "phi", "tau"),
-        decay_params = c("alpha"),
-        grouped_params = c("alpha", "alpha_decay"),
-        period_min = tau_in - 4, period_max = tau_in + 4
+  withr::with_seed(42, {
+    df <- make_data(k1 = 0, alpha1 = 10, phi1 = 0, seed = 42, hours = 96, noise_sd = 2)
+    df$time <- df$time / 24 * tau_in
+    alpha_decay1_in <- 0.05
+    # note that decay is on a scale of time in radians, not time in hours.
+    df$measure[df$group == "g2"] <- df$measure[df$group == "g2"] * exp(-alpha_decay1_in * df$time[df$group == "g2"])
+
+    out_alpha_decay <-
+      circacompare(
+        x = df, "time", "group", "measure", period = NA,
+        control = list(
+          main_params = c("k", "alpha", "phi", "tau"),
+          decay_params = c("alpha"),
+          grouped_params = c("alpha", "alpha_decay"),
+          period_min = tau_in - 4, period_max = tau_in + 4
+        )
       )
-    )
+  })
 
   fit_alpha_decay1 <- extract_model_coefs(out_alpha_decay$fit)["alpha_decay1", ]
   alpha_decay1_est <- fit_alpha_decay1["estimate"]
diff --git a/tests/testthat/test-circacompare_mixed.R b/tests/testthat/test-circacompare_mixed.R
index a43b435..b216d20 100644
--- a/tests/testthat/test-circacompare_mixed.R
+++ b/tests/testthat/test-circacompare_mixed.R
@@ -1,5 +1,4 @@
 test_that("circacompare_mixed() works", {
-  set.seed(99)
   phi1_in <- 3.15
   mixed_data <- function(n) {
     counter <- 1
@@ -16,16 +15,18 @@ test_that("circacompare_mixed() works", {
     return(res)
   }
 
-  df <- mixed_data(20)
+  withr::with_seed(99, {
+    df <- mixed_data(20)
 
-  out <- circacompare_mixed(
-    x = df,
-    col_time = "time",
-    col_group = "group",
-    col_outcome = "measure",
-    col_id = "id",
-    control = list(grouped_params = c("phi"), random_params = c("phi1"))
-  )
+    out <- circacompare_mixed(
+      x = df,
+      col_time = "time",
+      col_group = "group",
+      col_outcome = "measure",
+      col_id = "id",
+      control = list(grouped_params = c("phi"), random_params = c("phi1"))
+    )
+  })
 
   phi1_fit <- extract_model_coefs(out$fit)["phi1", ]
   phi1_est <- phi1_fit["estimate"]
@@ -36,7 +37,6 @@ test_that("circacompare_mixed() works", {
 
 ### make test that weights are used correctly and malformatted weights are detected
 test_that("weights work", {
-  set.seed(99)
   phi1_in <- 3.15
   mixed_data <- function(n) {
     counter <- 1
@@ -53,82 +53,76 @@ test_that("weights work", {
     return(res)
   }
 
-  df <- mixed_data(20)
-
-  out <- circacompare_mixed(
-    x = df,
-    col_time = "time",
-    col_group = "group",
-    col_outcome = "measure",
-    col_id = "id",
-    control = list(grouped_params = c("phi"), random_params = c("phi1"))
-  )
+  withr::with_seed(99, {
+    df <- mixed_data(20)
 
-  # no weights used (= all weights are 1), hence fit$apVar should not be populated
-  out <- circacompare_mixed(
-    x = df,
-    col_time = "time",
-    col_group = "group",
-    col_outcome = "measure",
-    col_id = "id",
-    control = list(grouped_params = c("phi"), random_params = c("phi1"))
-  )
-  expect_true(is(out$fit$apVar, "character"))
-
-  # when weights are not all 1 then fit$apVar should be a matrix
-  sw <- runif(n = nrow(df))
-  out2 <- circacompare_mixed(
-    x = df,
-    col_time = "time",
-    col_group = "group",
-    col_outcome = "measure",
-    col_id = "id",
-    control = list(grouped_params = c("phi"), random_params = c("phi1")),
-    weights = sw
-  )
-  expect_true(is(out2$fit$apVar, "matrix"))
-
-  # weights must be same length as nrow(x)
-  sw2 <- c(sw, 1)
-  expect_error(
-    circacompare_mixed(
+    # no weights used (= all weights are 1), hence fit$apVar should not be populated
+    out <- circacompare_mixed(
       x = df,
       col_time = "time",
       col_group = "group",
       col_outcome = "measure",
       col_id = "id",
-      control = list(grouped_params = c("phi"), random_params = c("phi1")),
-      weights = sw2
+      control = list(grouped_params = c("phi"), random_params = c("phi1"))
     )
-  )
 
-  # weights must not contain NA
-  sw3 <- sw
-  sw3[1] <- NA
-  expect_error(
-    circacompare_mixed(
+    # when weights are not all 1 then fit$apVar should be a matrix
+    sw <- runif(n = nrow(df))
+    out2 <- circacompare_mixed(
       x = df,
       col_time = "time",
       col_group = "group",
       col_outcome = "measure",
       col_id = "id",
       control = list(grouped_params = c("phi"), random_params = c("phi1")),
-      weights = sw3
+      weights = sw
     )
-  )
 
-  # weights must not be negative
-  sw4 <- sw
-  sw4[1] <- -1
-  expect_error(
-    circacompare_mixed(
-      x = df,
-      col_time = "time",
-      col_group = "group",
-      col_outcome = "measure",
-      col_id = "id",
-      control = list(grouped_params = c("phi"), random_params = c("phi1")),
-      weights = sw4
+    sw2 <- c(sw, 1)
+    expect_error(
+      circacompare_mixed(
+        x = df,
+        col_time = "time",
+        col_group = "group",
+        col_outcome = "measure",
+        col_id = "id",
+        control = list(grouped_params = c("phi"), random_params = c("phi1")),
+        weights = sw2
+      )
+    )
+
+    # weights must not contain NA
+    sw3 <- sw
+    sw3[1] <- NA
+    expect_error(
+      circacompare_mixed(
+        x = df,
+        col_time = "time",
+        col_group = "group",
+        col_outcome = "measure",
+        col_id = "id",
+        control = list(grouped_params = c("phi"), random_params = c("phi1")),
+        weights = sw3
+      )
     )
-  )
+
+    # weights must not be negative
+    sw4 <- sw
+    sw4[1] <- -1
+    expect_error(
+      circacompare_mixed(
+        x = df,
+        col_time = "time",
+        col_group = "group",
+        col_outcome = "measure",
+        col_id = "id",
+        control = list(grouped_params = c("phi"), random_params = c("phi1")),
+        weights = sw4
+      )
+    )
+  })
+
+  expect_true(is(out$fit$apVar, "character"))
+  expect_true(is(out2$fit$apVar, "matrix"))
+
 })
diff --git a/vignettes/circacompare-vignette.Rmd b/vignettes/circacompare-vignette.Rmd
index 57bd8a2..c9c3343 100644
--- a/vignettes/circacompare-vignette.Rmd
+++ b/vignettes/circacompare-vignette.Rmd
@@ -52,7 +52,6 @@ head(data_single)
 In the case that you have data from two groups and you're wishing to determine the differences in mesor, amplitude, or phase between them, you will need an additional column (with two possible values) representing the groups.
 
 ```{r}
-set.seed(42)
 data_grouped <- make_data(phi1 = 6, noise_sd = 1)
 head(data_grouped)
 tail(data_grouped)
@@ -169,8 +168,7 @@ If you have repeated measures data, then it may be inappropriate to use a standa
 
 Here, I have some data that has some simulated within-id correlation in terms of the amount of phase shift between groups. 
 
-```{r}
-set.seed(42)
+```{r, warning=FALSE, message=FALSE}
 phi1_in <- 3
 mixed_data <- function(n) {
   counter <- 1
@@ -194,7 +192,8 @@ out <- circacompare_mixed(
   col_outcome = "measure",
   col_id = "id",
   control = list(grouped_params = c("alpha", "phi"), random_params = c("phi1")),
-  period = 24
+  period = 24,
+  suppress_all = TRUE
 )
 
 ggplot(data = df[df$id %in% c(1:6), ], aes(time, measure)) +
@@ -253,7 +252,6 @@ This produces the following model:
 `measure~k+((alpha+alpha1*x_group)*exp(-(alpha_decay+alpha_decay1*x_group)*time_r))*cos((24/(tau))*time_r-(phi))`
 
 ```{r}
-set.seed(42)
 tau_in <- runif(1, 8, 20)
 alpha_decay1_in <- runif(1, 0.02, 0.05)