diff --git a/R/HMMplot.R b/R/HMMplot.R
index 7ff1414f..18c2ad3f 100644
--- a/R/HMMplot.R
+++ b/R/HMMplot.R
@@ -95,7 +95,7 @@ HMMplot <- function(x, layout = "horizontal", pie = TRUE,
   # Vertex label distances
   if (is.character(vertex.label.dist)) {
     ind <- pmatch(vertex.label.dist, "auto")
-    stopifnot(
+    stopifnot_(
       !is.na(ind),
       "{.arg vertex.label.dist} only accepts the value {.val 'auto'} or a 
       numerical vector."
diff --git a/R/average_marginal_prediction.R b/R/average_marginal_prediction.R
index 44dc21bd..e45dc1e1 100644
--- a/R/average_marginal_prediction.R
+++ b/R/average_marginal_prediction.R
@@ -61,7 +61,7 @@ average_marginal_prediction <- function(
       "Can't find time variable {.var {variable}} in {.arg newdata}."
     )
   } else {
-    stopifnot(
+    stopifnot_(
       !is.null(model$data),
       "Model does not contain original data and argument {.arg newdata} is 
       {.var NULL}."
diff --git a/R/plot.hmm.R b/R/plot.hmm.R
index ba40a9e1..bcd02bf7 100644
--- a/R/plot.hmm.R
+++ b/R/plot.hmm.R
@@ -328,7 +328,7 @@ plot.hmm <- function(x, layout = "horizontal", pie = TRUE,
   # Vertex label distances
   if (is.character(vertex.label.dist)) {
     ind <- pmatch(vertex.label.dist, "auto")
-    stopifnot(
+    stopifnot_(
       !is.na(ind),
       "{.arg vertex.label.dist} only accepts the value {.val 'auto'} or a 
       numerical vector."
diff --git a/R/predict.R b/R/predict.R
index 6b25ad5d..2fdfd362 100644
--- a/R/predict.R
+++ b/R/predict.R
@@ -1,6 +1,6 @@
 #' Predict method for non-homogeneous hidden Markov models
 #' 
-#' @param model A Hidden Markov Model of class `nhmm` or `mnhmm`.
+#' @param object A Hidden Markov Model of class `nhmm` or `mnhmm`.
 #' @param newdata Optional data frame which is used for prediction.
 #' @param nsim Non-negative integer defining the number of samples from the
 #' normal approximation of the model coefficients.
@@ -9,7 +9,7 @@
 #' quantiles.
 #' @export
 predict.nhmm <- function(
-    model, newdata = NULL, nsim = 0, 
+    object, newdata = NULL, nsim = 0, 
     probs = c(0.025, 0.5, 0.975), return_samples = FALSE) {
   
   stopifnot_(
@@ -17,8 +17,8 @@ predict.nhmm <- function(
     "Argument {.arg nsim} should be a single non-negative integer."
   )
   if (!is.null(newdata)) {
-    time <- model$time_variable
-    id <- model$id_variable
+    time <- object$time_variable
+    id <- object$id_variable
     stopifnot_(
       is.data.frame(newdata), 
       "Argument {.arg newdata} must be a {.cls data.frame} object."
@@ -32,44 +32,44 @@ predict.nhmm <- function(
       "Can't find time index variable {.var {time}} in {.arg newdata}."
     )
   } else {
-    stopifnot(
-      !is.null(model$data),
+    stopifnot_(
+      !is.null(object$data),
       "Model does not contain original data and argument {.arg newdata} is 
       {.var NULL}."
     )
-    model <- update(model, newdata = newdata)
+    object <- update(object, newdata = newdata)
   }
   
   beta_i_raw <- stan_to_cpp_initial(
-    model$estimation_results$parameters$beta_i_raw
+    object$estimation_results$parameters$beta_i_raw
   )
   beta_s_raw <- stan_to_cpp_transition(
-    model$estimation_results$parameters$beta_s_raw
+    object$estimation_results$parameters$beta_s_raw
   )
   beta_o_raw <- stan_to_cpp_emission(
-    model$estimation_results$parameters$beta_o_raw,
+    object$estimation_results$parameters$beta_o_raw,
     1,
-    model$n_channels > 1
+    object$n_channels > 1
   )
-  X_initial <- t(model$X_initial)
-  X_transition <- aperm(model$X_transition, c(3, 1, 2))
-  X_emission <- aperm(model$X_emission, c(3, 1, 2))
+  X_initial <- t(object$X_initial)
+  X_transition <- aperm(object$X_transition, c(3, 1, 2))
+  X_emission <- aperm(object$X_emission, c(3, 1, 2))
   out <- list()
   out$pi <- get_pi(beta_i_raw, X_initial, 0)
   out$A <- get_A(beta_s_raw, X_transition, 0)
-  out$B <- if (model$n_channels == 1) {
+  out$B <- if (object$n_channels == 1) {
     get_B(beta_o_raw, X_emission, 0)
   } else {
     get_multichannel_B(
       beta_o_raw, 
       X_emission, 
-      model$n_states, 
-      model$n_channels, 
-      model$n_symbols, 
+      object$n_states, 
+      object$n_channels, 
+      object$n_symbols, 
       0, 0)
   }
   if (nsim > 0) {
-    samples <- sample_parameters(model, nsim, probs, return_samples)
+    samples <- sample_parameters(object, nsim, probs, return_samples)
     if (return_samples) {
       out$samples <- samples
     } else {
@@ -80,7 +80,7 @@ predict.nhmm <- function(
 }
 #' @export
 predict.mnhmm <- function(
-    model, newdata = NULL, nsim = 0, 
+    object, newdata = NULL, nsim = 0, 
     probs = c(0.025, 0.5, 0.975), return_samples = FALSE) {
   
   stopifnot_(
@@ -88,8 +88,8 @@ predict.mnhmm <- function(
     "Argument {.arg nsim} should be a single non-negative integer."
   )
   if (!is.null(newdata)) {
-    time <- model$time_variable
-    id <- model$id_variable
+    time <- object$time_variable
+    id <- object$id_variable
     stopifnot_(
       is.data.frame(newdata), 
       "Argument {.arg newdata} must be a {.cls data.frame} object."
@@ -103,48 +103,48 @@ predict.mnhmm <- function(
       "Can't find time index variable {.var {time}} in {.arg newdata}."
     )
   } else {
-    stopifnot(
-      !is.null(model$data),
+    stopifnot_(
+      !is.null(object$data),
       "Model does not contain original data and argument {.arg newdata} is 
       {.var NULL}."
     )
-    model <- update(model, newdata = newdata)
+    object <- update(object, newdata = newdata)
   }
   
   beta_i_raw <- stan_to_cpp_initial(
-    model$estimation_results$parameters$beta_i_raw
+    object$estimation_results$parameters$beta_i_raw
   )
   beta_s_raw <- stan_to_cpp_transition(
-    model$estimation_results$parameters$beta_s_raw
+    object$estimation_results$parameters$beta_s_raw
   )
   beta_o_raw <- stan_to_cpp_emission(
-    model$estimation_results$parameters$beta_o_raw,
+    object$estimation_results$parameters$beta_o_raw,
     1,
-    model$n_channels > 1
+    object$n_channels > 1
   )
-  X_initial <- t(model$X_initial)
-  X_transition <- aperm(model$X_transition, c(3, 1, 2))
-  X_emission <- aperm(model$X_emission, c(3, 1, 2))
-  X_cluster <- t(model$X_cluster)
+  X_initial <- t(object$X_initial)
+  X_transition <- aperm(object$X_transition, c(3, 1, 2))
+  X_emission <- aperm(object$X_emission, c(3, 1, 2))
+  X_cluster <- t(object$X_cluster)
   out <- list()
   out$pi <- get_pi(beta_i_raw, X_initial, 0)
   out$A <- get_A(beta_s_raw, X_transition, 0)
-  out$B <- if (model$n_channels == 1) {
+  out$B <- if (object$n_channels == 1) {
     get_B(beta_o_raw, X_emission, 0)
   } else {
     get_multichannel_B(
       beta_o_raw, 
       X_emission, 
-      model$n_states, 
-      model$n_channels, 
-      model$n_symbols, 
+      object$n_states, 
+      object$n_channels, 
+      object$n_symbols, 
       0, 0)
   }
   out$omega <- get_omega(
-    model$estimation_results$parameters$theta_raw, X_cluster, 0
+    object$estimation_results$parameters$theta_raw, X_cluster, 0
   )
   if (nsim > 0) {
-    samples <- sample_parameters(model, nsim, probs, return_samples)
+    samples <- sample_parameters(object, nsim, probs, return_samples)
     if (return_samples) {
       out$samples <- samples
     } else {
diff --git a/man/predict.nhmm.Rd b/man/predict.nhmm.Rd
index c6dd5881..eec0739d 100644
--- a/man/predict.nhmm.Rd
+++ b/man/predict.nhmm.Rd
@@ -5,7 +5,7 @@
 \title{Predict method for non-homogeneous hidden Markov models}
 \usage{
 \method{predict}{nhmm}(
-  model,
+  object,
   newdata = NULL,
   nsim = 0,
   probs = c(0.025, 0.5, 0.975),
@@ -13,7 +13,7 @@
 )
 }
 \arguments{
-\item{model}{A Hidden Markov Model of class \code{nhmm} or \code{mnhmm}.}
+\item{object}{A Hidden Markov Model of class \code{nhmm} or \code{mnhmm}.}
 
 \item{newdata}{Optional data frame which is used for prediction.}