diff --git a/.github/CONTRIBUTING.md b/.github/CONTRIBUTING.md
index 4a229b60d..592893bd0 100644
--- a/.github/CONTRIBUTING.md
+++ b/.github/CONTRIBUTING.md
@@ -1,6 +1,6 @@
 # Contributing to parameters
 
-This outlines how to propose a change to **parameters**. 
+This outlines how to propose a change to **parameters**.
 
 ## Fixing typos
 
@@ -9,7 +9,7 @@ Small typos or grammatical errors in documentation may be edited directly using
 ## Filing an issue
 
 The easiest way to propose a change or new feature is to file an issue. If you've found a
-bug, you may also create an associated issue. If possible, try to illustrate your proposal or the bug with a minimal [reproducible example](https://www.tidyverse.org/help/#reprex).
+bug, you may also create an associated issue. If possible, try to illustrate your proposal or the bug with a minimal reproducible example.
 
 ## Pull requests
 
diff --git a/.github/SUPPORT.md b/.github/SUPPORT.md
index eaba28234..5686c4ea2 100644
--- a/.github/SUPPORT.md
+++ b/.github/SUPPORT.md
@@ -7,9 +7,7 @@ Start by making a minimal **repr**oducible **ex**ample using the
 [reprex](http://reprex.tidyverse.org/) package. If you haven't heard of or used
 reprex before, you're in for a treat! Seriously, reprex will make all of your
 R-question-asking endeavors easier (which is a pretty insane ROI for the five to
-ten minutes it'll take you to learn what it's all about). For additional reprex
-pointers, check out the [Get help!](https://www.tidyverse.org/help/) resource
-used by the tidyverse team.
+ten minutes it'll take you to learn what it's all about).
 
 Armed with your reprex, the next step is to figure out where to ask:
 
@@ -26,4 +24,4 @@ default, the search will be pre-populated with `is:issue is:open`. You can
 (e.g. `is:pr`, `is:closed`) as needed. For example, you'd simply
 remove `is:open` to search _all_ issues in the repo, open or closed.
 
-Thanks for your help!
\ No newline at end of file
+Thanks for your help!
diff --git a/DESCRIPTION b/DESCRIPTION
index de25c80fe..0bb253371 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Type: Package
 Package: parameters
 Title: Processing of Model Parameters
-Version: 0.28.2.7
+Version: 0.28.2.8
 Authors@R:
     c(person(given = "Daniel",
              family = "Lüdecke",
diff --git a/NAMESPACE b/NAMESPACE
index 7b593c795..b6e5885bc 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -114,6 +114,8 @@ S3method(ci,systemfit)
 S3method(ci,varest)
 S3method(ci,zerocount)
 S3method(ci,zeroinfl)
+S3method(ci_kenward,default)
+S3method(ci_kenward,glmmTMB)
 S3method(cluster_discrimination,cluster_analysis)
 S3method(cluster_discrimination,default)
 S3method(cluster_performance,dbscan)
@@ -138,7 +140,6 @@ S3method(display,parameters_pca_summary)
 S3method(display,parameters_sem)
 S3method(display,parameters_simulate)
 S3method(display,parameters_standardized)
-S3method(dof_satterthwaite,lmerMod)
 S3method(equivalence_test,MixMod)
 S3method(equivalence_test,estimate_contrasts)
 S3method(equivalence_test,estimate_means)
@@ -568,7 +569,6 @@ S3method(p_value,wbm)
 S3method(p_value,zcpglm)
 S3method(p_value,zerocount)
 S3method(p_value,zeroinfl)
-S3method(p_value_kenward,lmerMod)
 S3method(plot,cluster_analysis)
 S3method(plot,cluster_analysis_summary)
 S3method(plot,compare_parameters)
@@ -660,6 +660,8 @@ S3method(reduce_parameters,lm)
 S3method(reduce_parameters,merMod)
 S3method(reshape_loadings,data.frame)
 S3method(reshape_loadings,parameters_efa)
+S3method(se_kenward,default)
+S3method(se_kenward,glmmTMB)
 S3method(se_satterthwaite,default)
 S3method(select_parameters,lm)
 S3method(select_parameters,merMod)
diff --git a/NEWS.md b/NEWS.md
index 6e74de38f..4126ad5d6 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -4,6 +4,11 @@
 
 * Added support for package *lcmm*.
 
+* Added `ci_method` options `"kenward-roger"` and `"satterthwaite"` for models
+  from package *glmmTMB*. Consequently, `se_kenward()`, `se_satterthwaite()`,
+  `ci_kenward()`, `ci_satterthwaite()`, `p_value_kenward()` and
+  `p_value_satterthwaite()` can now be used with `glmmTMB` models.
+
 # parameters 0.28.2
 
 ## Bug fixes
diff --git a/R/ci_generic.R b/R/ci_generic.R
index 9a9341df0..cc94de6ea 100644
--- a/R/ci_generic.R
+++ b/R/ci_generic.R
@@ -1,19 +1,22 @@
 # generic function for CI calculation
-.ci_generic <- function(model,
-                        ci = 0.95,
-                        method = "wald",
-                        dof = NULL,
-                        effects = "fixed",
-                        component = "all",
-                        vcov = NULL,
-                        vcov_args = NULL,
-                        verbose = TRUE,
-                        ...) {
+.ci_generic <- function(
+  model,
+  ci = 0.95,
+  method = "wald",
+  dof = NULL,
+  effects = "fixed",
+  component = "all",
+  vcov = NULL,
+  vcov_args = NULL,
+  verbose = TRUE,
+  ...
+) {
   # check method
   if (is.null(method)) {
     method <- "wald"
   }
   method <- tolower(method)
+  # fmt: skip
   method <- insight::validate_argument(
     method,
     c(
@@ -23,6 +26,7 @@
   )
 
   effects <- insight::validate_argument(effects, c("fixed", "random", "all"))
+  # fmt: skip
   component <- insight::validate_argument(
     component,
     c(
@@ -31,7 +35,8 @@
     )
   )
 
-  if (method == "ml1") { # nolint
+  if (method == "ml1") {
+    # nolint
     return(ci_ml1(model, ci = ci))
   } else if (method == "betwithin") {
     return(ci_betwithin(model, ci = ci))
diff --git a/R/ci_kenward.R b/R/ci_kenward.R
index 855bf5868..57e1644b5 100644
--- a/R/ci_kenward.R
+++ b/R/ci_kenward.R
@@ -1,7 +1,15 @@
 #' @rdname p_value_kenward
 #' @export
-ci_kenward <- function(model, ci = 0.95) {
-  .check_REML_fit(model)
+ci_kenward <- function(model, ci = 0.95, ...) {
+  UseMethod("ci_kenward")
+}
+
+#' @export
+ci_kenward.default <- function(model, ci = 0.95, ...) {
+  if (!.check_REML_fit(model)) {
+    model <- stats::update(model, . ~ ., REML = TRUE)
+  }
+
   df_kr <- dof_kenward(model)
   out <- lapply(ci, function(i) {
     .ci_dof(
@@ -19,15 +27,43 @@ ci_kenward <- function(model, ci = 0.95) {
   out
 }
 
+#' @export
+ci_kenward.glmmTMB <- function(model, ci = 0.95, ...) {
+  if (!.check_REML_fit(model)) {
+    model <- stats::update(model, . ~ ., REML = TRUE)
+  }
+
+  df_kr <- dof_kenward(model)
+  out <- lapply(ci, function(i) {
+    .ci_dof(
+      model = model,
+      ci = i,
+      dof = df_kr,
+      effects = "fixed",
+      component = "conditional", # for glmmTMB, only conditional
+      method = "kenward",
+      se = attr(df_kr, "se", exact = TRUE)
+    )
+  })
+  out <- do.call(rbind, out)
+  row.names(out) <- NULL
+  out
+}
+
 
 .ci_kenward_dof <- function(model, ci = 0.95, df_kr) {
+  if (inherits(model, "glmmTMB")) {
+    component <- "conditional"
+  } else {
+    component <- "all"
+  }
   out <- lapply(ci, function(i) {
     .ci_dof(
       model = model,
       ci = i,
       dof = df_kr$df_error,
       effects = "fixed",
-      component = "all",
+      component = component,
       method = "kenward",
       se = df_kr$SE
     )
diff --git a/R/ci_satterthwaite.R b/R/ci_satterthwaite.R
index 2f5ebe423..89631643e 100644
--- a/R/ci_satterthwaite.R
+++ b/R/ci_satterthwaite.R
@@ -2,13 +2,18 @@
 #' @export
 ci_satterthwaite <- function(model, ci = 0.95, ...) {
   df_satter <- dof_satterthwaite(model)
+  if (inherits(model, "glmmTMB")) {
+    component <- "conditional"
+  } else {
+    component <- "all"
+  }
   out <- lapply(ci, function(i) {
     .ci_dof(
       model = model,
       ci = i,
       dof = df_satter,
       effects = "fixed",
-      component = "all",
+      component = component,
       method = "satterthwaite",
       ...
     )
diff --git a/R/dof.R b/R/dof.R
index cb84690aa..bb3d9293f 100644
--- a/R/dof.R
+++ b/R/dof.R
@@ -154,7 +154,7 @@ dof <- degrees_of_freedom
       ))
     }
     return(FALSE)
-}
+  }
 
   if (!info$is_linear && method %in% c("satterthwaite", "kenward", "kr")) {
     if (verbose) {
diff --git a/R/dof_kenward.R b/R/dof_kenward.R
index 1cb951186..6ab28c089 100644
--- a/R/dof_kenward.R
+++ b/R/dof_kenward.R
@@ -1,293 +1,5 @@
 #' @rdname p_value_kenward
 #' @export
 dof_kenward <- function(model) {
-  parameters <- insight::find_parameters(model, effects = "fixed", flatten = TRUE)
-  L <- as.data.frame(diag(rep(1, n_parameters(model, effects = "fixed"))))
-  krvcov <- .vcov_kenward_ajusted(model)
-
-  dof <- stats::setNames(sapply(L, .kenward_adjusted_ddf, model = model, adjusted_vcov = krvcov), parameters)
-
-  attr(dof, "vcov") <- krvcov
-  attr(dof, "se") <- abs(as.vector(sqrt(diag(as.matrix(krvcov)))))
-  dof
-}
-
-
-# The following code was taken from the "pbkrtest" package and slightly modified
-#' @author Søren Højsgaard, \email{sorenh@@math.aau.dk}
-.kenward_adjusted_ddf <- function(model, linear_coef, adjusted_vcov) {
-  .adjusted_ddf(adjusted_vcov, linear_coef, stats::vcov(model))
-}
-
-
-.adjusted_ddf <- function(adjusted_vcov, linear_coef, unadjusted_vcov = adjusted_vcov) {
-  insight::check_if_installed("Matrix")
-
-  if (!is.matrix(linear_coef)) {
-    linear_coef <- matrix(linear_coef, ncol = 1)
-  }
-  vlb <- sum(linear_coef * (unadjusted_vcov %*% linear_coef))
-  theta <- Matrix::Matrix(as.numeric(outer(linear_coef, linear_coef) / vlb), nrow = length(linear_coef))
-  P <- attr(adjusted_vcov, "P")
-  W <- attr(adjusted_vcov, "W")
-
-  A1 <- A2 <- 0
-  theta_unadjusted_vcov <- theta %*% unadjusted_vcov
-  n.ggamma <- length(P)
-  for (ii in 1:n.ggamma) {
-    for (jj in ii:n.ggamma) {
-      if (ii == jj) {
-        e <- 1
-      } else {
-        e <- 2
-      }
-      ui <- as.matrix(theta_unadjusted_vcov %*% P[[ii]] %*% unadjusted_vcov)
-      uj <- as.matrix(theta_unadjusted_vcov %*% P[[jj]] %*% unadjusted_vcov)
-      A1 <- A1 + e * W[ii, jj] * (sum(diag(ui)) * sum(diag(uj)))
-      A2 <- A2 + e * W[ii, jj] * sum(ui * t(uj))
-    }
-  }
-
-  B <- (A1 + 6 * A2) / 2
-  g <- (2 * A1 - 5 * A2) / (3 * A2)
-  c1 <- g / (3 + 2 * (1 - g))
-  c2 <- (1 - g) / (3 + 2 * (1 - g))
-  c3 <- (3 - g) / (3 + 2 * (1 - g))
-  EE <- 1 + A2
-  VV <- 2 * (1 + B)
-  EEstar <- 1 / (1 - A2)
-  VVstar <- 2 * ((1 + c1 * B) / ((1 - c2 * B)^2 * (1 - c3 * B)))
-  V0 <- 1 + c1 * B
-  V1 <- 1 - c2 * B
-  V2 <- 1 - c3 * B
-  V0 <- ifelse(abs(V0) < 1e-10, 0, V0)
-  rho <- (.divZero(1 - A2, V1))^2 * V0 / V2
-  df2 <- 4 + 3 / (rho - 1)
-  df2
-}
-
-
-.divZero <- function(x, y, tol = 1e-14) {
-  ## ratio x/y is set to 1 if both |x| and |y| are below tol
-  if (abs(x) < tol && abs(y) < tol) {
-    1
-  } else {
-    x / y
-  }
-}
-
-
-.vcov_kenward_ajusted <- function(model) {
-  insight::check_if_installed("lme4")
-
-  if (!(lme4::getME(model, "is_REML"))) {
-    model <- stats::update(model, . ~ ., REML = TRUE)
-  }
-  .vcovAdj16_internal(stats::vcov(model), .get_SigmaG(model), lme4::getME(model, "X"))
-}
-
-
-.get_SigmaG <- function(model) {
-  insight::check_if_installed("lme4")
-  insight::check_if_installed("Matrix")
-
-  GGamma <- lme4::VarCorr(model)
-  SS <- .shgetME(model)
-
-  ## Put covariance parameters for the random effects into a vector:
-  ## TODO: It is a bit ugly to throw everything into one long vector here; a list would be more elegant
-  ggamma <- NULL
-  for (ii in 1:(SS$n.RT)) {
-    Lii <- GGamma[[ii]]
-    ggamma <- c(ggamma, Lii[lower.tri(Lii, diag = TRUE)])
-  }
-  ggamma <- c(ggamma, stats::sigma(model)^2) ## Extend ggamma by the residuals variance
-  n.ggamma <- length(ggamma)
-
-  ## Find G_r:
-  G <- NULL
-  Zt <- lme4::getME(model, "Zt")
-  for (ss in 1:SS$n.RT) {
-    ZZ <- .shget_Zt_group(ss, Zt, SS$Gp)
-    n.lev <- SS$n.lev.by.RT2[ss] ## ; cat(sprintf("n.lev=%i\n", n.lev))
-    Ig <- Matrix::sparseMatrix(1:n.lev, 1:n.lev, x = 1)
-    for (rr in 1:SS$n.parm.by.RT[ss]) {
-      ## This is takes care of the case where there is random regression and several matrices have to be constructed.
-      ## FIXME: I am not sure this is correct if there is a random quadratic term. The '2' below looks suspicious.
-      ii.jj <- .index2UpperTriEntry(rr, SS$n.comp.by.RT[ss]) ## ; cat("ii.jj:"); print(ii.jj)
-      ii.jj <- unique(ii.jj)
-      if (length(ii.jj) == 1) {
-        EE <- Matrix::sparseMatrix(
-          ii.jj,
-          ii.jj,
-          x = 1,
-          dims = rep(SS$n.comp.by.RT[ss], 2)
-        )
-      } else {
-        EE <- Matrix::sparseMatrix(ii.jj, ii.jj[2:1], dims = rep(SS$n.comp.by.RT[ss], 2))
-      }
-      EE <- Ig %x% EE ## Kronecker product
-      G <- c(G, list(t(ZZ) %*% EE %*% ZZ))
-    }
-  }
-
-  ## Extend by the indentity for the residual
-  n.obs <- insight::n_obs(model)
-  G <- c(G, list(Matrix::sparseMatrix(1:n.obs, 1:n.obs, x = 1)))
-
-  Sigma <- ggamma[1] * G[[1]]
-  for (ii in 2:n.ggamma) {
-    Sigma <- Sigma + ggamma[ii] * G[[ii]]
-  }
-
-  list(Sigma = Sigma, G = G, n.ggamma = n.ggamma)
-}
-
-
-.index2UpperTriEntry <- function(k, N) {
-  ## inverse of indexSymmat2vec
-  ## result: index pair (i,j) with i>=j
-  ## k: element in the vector of upper triangular elements
-  ## example: N=3: k=1 -> (1,1), k=2 -> (1,2), k=3 -> (1,3), k=4 -> (2,2)
-  aa <- cumsum(N:1)
-  aaLow <- c(0, aa[-length(aa)])
-  i <- which(aaLow < k & k <= aa)
-  j <- k - N * i + N - i * (3 - i) / 2 + i
-  c(i, j)
-}
-
-
-.vcovAdj16_internal <- function(Phi, SigmaG, X) {
-  insight::check_if_installed("MASS")
-  insight::check_if_installed("Matrix")
-
-  SigmaInv <- chol2inv(chol(Matrix::forceSymmetric(as.matrix(SigmaG$Sigma))))
-  n.ggamma <- SigmaG$n.ggamma
-  TT <- as.matrix(SigmaInv %*% X)
-  HH <- OO <- vector("list", n.ggamma)
-
-  for (ii in 1:n.ggamma) {
-    HH[[ii]] <- as.matrix(SigmaG$G[[ii]] %*% SigmaInv)
-    OO[[ii]] <- as.matrix(HH[[ii]] %*% X)
-  }
-
-  ## Finding PP, QQ
-  PP <- QQ <- NULL
-  for (rr in 1:n.ggamma) {
-    OrTrans <- t(OO[[rr]])
-    PP <- c(PP, list(Matrix::forceSymmetric(-1 * OrTrans %*% TT)))
-    for (ss in rr:n.ggamma) {
-      QQ <- c(QQ, list(OrTrans %*% SigmaInv %*% OO[[ss]]))
-    }
-  }
-
-  PP <- as.matrix(PP)
-  QQ <- as.matrix(QQ)
-
-  Ktrace <- matrix(NA, nrow = n.ggamma, ncol = n.ggamma)
-  for (rr in 1:n.ggamma) {
-    HrTrans <- t(HH[[rr]])
-    for (ss in rr:n.ggamma) {
-      Ktrace[rr, ss] <- Ktrace[ss, rr] <- sum(HrTrans * HH[[ss]])
-    }
-  }
-
-  ## Finding information matrix
-  IE2 <- matrix(NA, nrow = n.ggamma, ncol = n.ggamma)
-  for (ii in 1:n.ggamma) {
-    Phi.P.ii <- Phi %*% PP[[ii]]
-    for (jj in ii:n.ggamma) {
-      www <- .indexSymmat2vec(ii, jj, n.ggamma)
-      IE2[ii, jj] <- IE2[jj, ii] <- Ktrace[ii, jj] -
-        2 * sum(Phi * QQ[[www]]) + sum(Phi.P.ii * (PP[[jj]] %*% Phi))
-    }
-  }
-
-  eigenIE2 <- eigen(IE2, only.values = TRUE)$values
-  condi <- min(abs(eigenIE2))
-
-  WW <- if (condi > 1e-10) {
-    as.matrix(Matrix::forceSymmetric(2 * solve(IE2)))
-  } else {
-    as.matrix(Matrix::forceSymmetric(2 * MASS::ginv(IE2)))
-  }
-
-  UU <- matrix(0, nrow = ncol(X), ncol = ncol(X))
-  for (ii in 1:(n.ggamma - 1)) {
-    for (jj in (ii + 1):n.ggamma) {
-      www <- .indexSymmat2vec(ii, jj, n.ggamma)
-      UU <- UU + WW[ii, jj] * (QQ[[www]] - PP[[ii]] %*% Phi %*% PP[[jj]])
-    }
-  }
-
-  UU <- as.matrix(UU)
-  UU <- UU + t(UU)
-  for (ii in 1:n.ggamma) {
-    www <- .indexSymmat2vec(ii, ii, n.ggamma)
-    UU <- UU + WW[ii, ii] * (QQ[[www]] - PP[[ii]] %*% Phi %*% PP[[ii]])
-  }
-
-  GGAMMA <- Phi %*% UU %*% Phi
-  PhiA <- Phi + 2 * GGAMMA
-  attr(PhiA, "P") <- PP
-  attr(PhiA, "W") <- WW
-  attr(PhiA, "condi") <- condi
-  PhiA
-}
-
-
-.indexSymmat2vec <- function(i, j, N) {
-  ## S[i,j] symetric N times N matrix
-  ## r the vector of upper triangular element  in row major order:
-  ## r= c(S[1,1],S[1,2]...,S[1,j], S[1,N], S[2,2],...S[N,N]
-  ## Result: k: index of k-th element of r
-  k <- if (i <= j) {
-    (i - 1) * (N - i / 2) + j
-  } else {
-    (j - 1) * (N - j / 2) + i
-  }
-}
-
-
-.shgetME <- function(model) {
-  insight::check_if_installed("lme4")
-
-  Gp <- lme4::getME(model, "Gp")
-  n.RT <- length(Gp) - 1 ## Number of random terms (i.e. of (|)'s)
-  n.lev.by.RT <- sapply(lme4::getME(model, "flist"), nlevels)
-  n.comp.by.RT <- .get.RT.dim.by.RT(model)
-  n.parm.by.RT <- (n.comp.by.RT + 1) * n.comp.by.RT / 2
-  n.RE.by.RT <- diff(Gp)
-
-  n.lev.by.RT2 <- n.RE.by.RT / n.comp.by.RT ## Same as n.lev.by.RT2 ???
-
-  list(
-    Gp = Gp, ## group.index
-    n.RT = n.RT, ## n.groupFac
-    n.lev.by.RT = n.lev.by.RT, ## nn.groupFacLevelsNew
-    n.comp.by.RT = n.comp.by.RT, ## nn.GGamma
-    n.parm.by.RT = n.parm.by.RT, ## mm.GGamma
-    n.RE.by.RT = n.RE.by.RT, ## ... Not returned before
-    n.lev.by.RT2 = n.lev.by.RT2, ## nn.groupFacLevels
-    n_rtrms = lme4::getME(model, "n_rtrms")
-  )
-}
-
-
-## Alternative to .get_Zt_group
-.shget_Zt_group <- function(ii.group, Zt, Gp, ...) {
-  zIndex.sub <- (Gp[ii.group] + 1):Gp[ii.group + 1]
-  as.matrix(Zt[zIndex.sub, ])
-}
-
-
-.get.RT.dim.by.RT <- function(model) {
-  insight::check_if_installed("lme4")
-
-  ## output: dimension (no of columns) of covariance matrix for random term ii
-  if (inherits(model, "mer")) {
-    vapply(model@ST, nrow, numeric(1))
-  } else {
-    lengths(lme4::getME(model, "cnms"))
-  }
+  insight::get_df(model, "kenward")
 }
diff --git a/R/dof_satterthwaite.R b/R/dof_satterthwaite.R
index c7c1915e4..c0b49c94f 100644
--- a/R/dof_satterthwaite.R
+++ b/R/dof_satterthwaite.R
@@ -1,17 +1,5 @@
 #' @rdname p_value_satterthwaite
 #' @export
 dof_satterthwaite <- function(model) {
-  UseMethod("dof_satterthwaite")
-}
-
-
-#' @export
-dof_satterthwaite.lmerMod <- function(model) {
-  insight::check_if_installed("lmerTest")
-
-  parameters <- insight::find_parameters(model, effects = "fixed", flatten = TRUE)
-  lmerTest_model <- lmerTest::as_lmerModLmerTest(model)
-  s <- summary(lmerTest_model)
-
-  stats::setNames(as.vector(s$coefficients[, 3]), parameters)
+  insight::get_df(model, "satterthwaite")
 }
diff --git a/R/methods_glmmTMB.R b/R/methods_glmmTMB.R
index 75031197a..9f3fbd9a0 100644
--- a/R/methods_glmmTMB.R
+++ b/R/methods_glmmTMB.R
@@ -1,9 +1,7 @@
 # Package glmmTMB
 
-
 # model_parameters -----
 
-
 #' @title Parameters from Mixed Models
 #' @name model_parameters.glmmTMB
 #'
@@ -167,31 +165,40 @@
 #' }
 #' @return A data frame of indices related to the model's parameters.
 #' @export
-model_parameters.glmmTMB <- function(model,
-                                     ci = 0.95,
-                                     ci_method = "wald",
-                                     ci_random = NULL,
-                                     bootstrap = FALSE,
-                                     iterations = 1000,
-                                     standardize = NULL,
-                                     effects = "all",
-                                     component = "all",
-                                     group_level = FALSE,
-                                     exponentiate = FALSE,
-                                     p_adjust = NULL,
-                                     vcov = NULL,
-                                     vcov_args = NULL,
-                                     wb_component = FALSE,
-                                     include_info = getOption("parameters_mixed_info", FALSE),
-                                     include_sigma = FALSE,
-                                     keep = NULL,
-                                     drop = NULL,
-                                     verbose = TRUE,
-                                     ...) {
+model_parameters.glmmTMB <- function(
+  model,
+  ci = 0.95,
+  ci_method = "wald",
+  ci_random = NULL,
+  bootstrap = FALSE,
+  iterations = 1000,
+  standardize = NULL,
+  effects = "all",
+  component = "all",
+  group_level = FALSE,
+  exponentiate = FALSE,
+  p_adjust = NULL,
+  vcov = NULL,
+  vcov_args = NULL,
+  wb_component = FALSE,
+  include_info = getOption("parameters_mixed_info", FALSE),
+  include_sigma = FALSE,
+  keep = NULL,
+  drop = NULL,
+  verbose = TRUE,
+  ...
+) {
   insight::check_if_installed("glmmTMB")
 
   # p-values, CI and se might be based on different df-methods
-  ci_method <- .check_df_method(ci_method)
+  ci_method <- insight::validate_argument(
+    ci_method,
+    # fmt: skip
+    c(
+      "wald", "normal", "residual", "ml1", "betwithin", "satterthwaite",
+      "kenward", "kr", "boot", "profile", "uniroot", "robust"
+    )
+  )
 
   # which components to return?
   effects <- insight::validate_argument(
@@ -241,6 +248,11 @@ model_parameters.glmmTMB <- function(model,
     component <- "conditional"
   }
 
+  # for ci_method kenward or satterthwaite, only conditional component
+  if (ci_method %in% c("satterthwaite", "kenward", "kr") && component != "conditional") {
+    component <- "conditional"
+  }
+
   # initialize
   params <- att <- NULL
   dispersion_param <- FALSE
@@ -297,7 +309,6 @@ model_parameters.glmmTMB <- function(model,
     att <- attributes(params)
   }
 
-
   # add random effects, either group level or re variances
   # ======================================================
 
@@ -357,10 +368,16 @@ model_parameters.glmmTMB <- function(model,
 
 # helper -----------------------------------------
 
-
 # this functions adds the dispersion parameter, if it is not already
 # present in the output
-.add_dispersion_param_glmmTMB <- function(model, params, effects, component, ci, verbose) {
+.add_dispersion_param_glmmTMB <- function(
+  model,
+  params,
+  effects,
+  component,
+  ci,
+  verbose
+) {
   dispersion_param <- FALSE
   if (
     # must be glmmTMB
@@ -410,16 +427,18 @@ model_parameters.glmmTMB <- function(model,
 # whether group level estimates or random effects variances are requested,
 # the related parameters are returned. It also correctly deals with the
 # dispersion parameter, if present in random effects
-.add_random_effects_glmmTMB <- function(model,
-                                        params,
-                                        ci,
-                                        ci_method,
-                                        ci_random,
-                                        effects,
-                                        component,
-                                        dispersion_param,
-                                        group_level,
-                                        verbose = TRUE) {
+.add_random_effects_glmmTMB <- function(
+  model,
+  params,
+  ci,
+  ci_method,
+  ci_random,
+  effects,
+  component,
+  dispersion_param,
+  group_level,
+  verbose = TRUE
+) {
   params_random <- params_variance <- NULL
   random_effects <- insight::find_random(model, flatten = TRUE)
 
@@ -496,18 +515,24 @@ model_parameters.glmmTMB <- function(model,
 # ci -----
 
 #' @export
-ci.glmmTMB <- function(x,
-                       ci = 0.95,
-                       dof = NULL,
-                       method = "wald",
-                       component = "all",
-                       verbose = TRUE,
-                       ...) {
-  method <- tolower(method)
+ci.glmmTMB <- function(
+  x,
+  ci = 0.95,
+  dof = NULL,
+  method = "wald",
+  component = "all",
+  verbose = TRUE,
+  ...
+) {
   method <- insight::validate_argument(
-    method,
-    c("wald", "normal", "ml1", "betwithin", "profile", "uniroot", "robust", "residual")
+    tolower(method),
+    # fmt: skip
+    c(
+      "wald", "normal", "residual", "ml1", "betwithin", "satterthwaite",
+      "kenward", "kr", "boot", "profile", "uniroot", "robust"
+    )
   )
+
   component <- insight::validate_argument(
     component,
     c("all", "conditional", "zi", "zero_inflated", "dispersion")
@@ -524,16 +549,27 @@ ci.glmmTMB <- function(x,
     } else {
       pp <- NULL
     }
-    out <- lapply(ci, function(i) .ci_profile_glmmTMB(x, ci = i, profiled = pp, component = component, ...))
+    out <- lapply(ci, function(i) {
+      .ci_profile_glmmTMB(x, ci = i, profiled = pp, component = component, ...)
+    })
     do.call(rbind, out)
 
     # uniroot CIs
   } else if (method == "uniroot") {
-    out <- lapply(ci, function(i) .ci_uniroot_glmmTMB(x, ci = i, component = component, ...))
+    out <- lapply(ci, function(i) {
+      .ci_uniroot_glmmTMB(x, ci = i, component = component, ...)
+    })
     do.call(rbind, out)
   } else {
     # all other
-    .ci_generic(model = x, ci = ci, dof = dof, method = method, component = component, ...)
+    .ci_generic(
+      model = x,
+      ci = ci,
+      dof = dof,
+      method = method,
+      component = component,
+      ...
+    )
   }
 }
 
@@ -541,45 +577,54 @@ ci.glmmTMB <- function(x,
 # standard_error -----
 
 #' @export
-standard_error.glmmTMB <- function(model,
-                                   effects = "fixed",
-                                   component = "all",
-                                   vcov = NULL,
-                                   vcov_args = NULL,
-                                   verbose = TRUE,
-                                   ...) {
+standard_error.glmmTMB <- function(
+  model,
+  effects = "fixed",
+  component = "all",
+  method = NULL,
+  vcov = NULL,
+  vcov_args = NULL,
+  verbose = TRUE,
+  ...
+) {
   component <- insight::validate_argument(
     component,
     c("all", "conditional", "zi", "zero_inflated", "dispersion")
   )
-  effects <- insight::validate_argument(
-    effects,
-    c("fixed", "random")
-  )
+  effects <- insight::validate_argument(effects, c("fixed", "random"))
 
   if (effects == "random") {
     .se_random_effects_glmmTMB(model)
   } else if (!is.null(vcov)) {
     .se_robust_glmmTMB(model, component, vcov, vcov_args, verbose, ...)
   } else {
-    .se_fixed_effects_glmmTMB(model, component, verbose)
+    .se_fixed_effects_glmmTMB(model, component, method, verbose)
   }
 }
 
 
 # helper --------------------------------------------------------------------
 
-
 # extract standard errors for fixed effects parameters
-.se_fixed_effects_glmmTMB <- function(model, component, verbose = TRUE) {
+.se_fixed_effects_glmmTMB <- function(model, component, method = NULL, verbose = TRUE) {
   if (is.null(.check_component(model, component, verbose = verbose))) {
     return(NULL)
   }
 
+  # kenward approx
+  if (!is.null(method) && method %in% c("kenward", "kr")) {
+    return(se_kenward(model, component = "conditional"))
+  }
+
   cs <- insight::compact_list(stats::coef(summary(model)))
   x <- lapply(names(cs), function(i) {
     .data_frame(
-      Parameter = insight::find_parameters(model, effects = "fixed", component = i, flatten = TRUE),
+      Parameter = insight::find_parameters(
+        model,
+        effects = "fixed",
+        component = i,
+        flatten = TRUE
+      ),
       SE = as.vector(cs[[i]][, 2]),
       Component = i
     )
@@ -595,18 +640,15 @@ standard_error.glmmTMB <- function(model,
 
 
 # extract robust standard errors for fixed effects parameters
-.se_robust_glmmTMB <- function(model,
-                               component = "all",
-                               vcov,
-                               vcov_args = NULL,
-                               verbose = TRUE,
-                               ...) {
-  fun_args <- list(
-    model,
-    component = component,
-    vcov = vcov,
-    vcov_args = vcov_args
-  )
+.se_robust_glmmTMB <- function(
+  model,
+  component = "all",
+  vcov,
+  vcov_args = NULL,
+  verbose = TRUE,
+  ...
+) {
+  fun_args <- list(model, component = component, vcov = vcov, vcov_args = vcov_args)
   fun_args <- c(fun_args, list(...))
   do.call("standard_error.default", fun_args)
 }
@@ -635,13 +677,14 @@ standard_error.glmmTMB <- function(model,
 
 # simulate model -----
 
-
 #' @export
-simulate_model.glmmTMB <- function(model,
-                                   iterations = 1000,
-                                   component = "all",
-                                   verbose = FALSE,
-                                   ...) {
+simulate_model.glmmTMB <- function(
+  model,
+  iterations = 1000,
+  component = "all",
+  verbose = FALSE,
+  ...
+) {
   component <- insight::validate_argument(
     component,
     c("all", "conditional", "zi", "zero_inflated", "dispersion")
@@ -656,7 +699,6 @@ simulate_model.glmmTMB <- function(model,
   has_zeroinflated <- !is.null(info) && isTRUE(info$is_zero_inflated)
   has_dispersion <- !is.null(info) && isTRUE(info$is_dispersion)
 
-
   # check component-argument ----
 
   if (component == "all") {
@@ -688,8 +730,9 @@ simulate_model.glmmTMB <- function(model,
     insight::format_error("No dispersion model found.")
   }
 
-
-  if (is.null(iterations)) iterations <- 1000
+  if (is.null(iterations)) {
+    iterations <- 1000
+  }
 
   if (all(component == c("conditional", "zero_inflated"))) {
     d1 <- .simulate_model(model, iterations, component = "conditional", ...)
@@ -725,13 +768,15 @@ simulate_model.glmmTMB <- function(model,
 # simulate_parameters -----
 
 #' @export
-simulate_parameters.glmmTMB <- function(model,
-                                        iterations = 1000,
-                                        centrality = "median",
-                                        ci = 0.95,
-                                        ci_method = "quantile",
-                                        test = "p-value",
-                                        ...) {
+simulate_parameters.glmmTMB <- function(
+  model,
+  iterations = 1000,
+  centrality = "median",
+  ci = 0.95,
+  ci_method = "quantile",
+  test = "p-value",
+  ...
+) {
   sim_data <- simulate_model(model, iterations = iterations, ...)
   out <- .summary_bootstrap(
     data = sim_data,
diff --git a/R/p_value_kenward.R b/R/p_value_kenward.R
index 388b16c9f..bb1c5017e 100644
--- a/R/p_value_kenward.R
+++ b/R/p_value_kenward.R
@@ -37,32 +37,42 @@
 #'   fixed effects from restricted maximum likelihood. Biometrics, 983-997.
 #' @export
 p_value_kenward <- function(model, dof = NULL) {
-  UseMethod("p_value_kenward")
-}
-
+  if (!.check_REML_fit(model)) {
+    model <- stats::update(model, . ~ ., REML = TRUE)
+  }
 
-#' @export
-p_value_kenward.lmerMod <- function(model, dof = NULL) {
   if (is.null(dof)) {
     dof <- dof_kenward(model)
   }
   .p_value_dof(model, dof, method = "kenward")
 }
 
-
 # helper ------------------------------
 
-.p_value_dof <- function(model,
-                         dof,
-                         method = NULL,
-                         statistic = NULL,
-                         se = NULL,
-                         component = c("all", "conditional", "zi", "zero_inflated", "dispersion", "precision", "scale", "smooth_terms", "full", "marginal"),
-                         effects = c("fixed", "random", "all"),
-                         verbose = TRUE,
-                         vcov = NULL,
-                         vcov_args = NULL,
-                         ...) {
+.p_value_dof <- function(
+  model,
+  dof,
+  method = NULL,
+  statistic = NULL,
+  se = NULL,
+  component = c(
+    "all",
+    "conditional",
+    "zi",
+    "zero_inflated",
+    "dispersion",
+    "precision",
+    "scale",
+    "smooth_terms",
+    "full",
+    "marginal"
+  ),
+  effects = c("fixed", "random", "all"),
+  verbose = TRUE,
+  vcov = NULL,
+  vcov_args = NULL,
+  ...
+) {
   component <- match.arg(component)
   effects <- match.arg(effects)
 
@@ -87,7 +97,8 @@ p_value_kenward.lmerMod <- function(model, dof = NULL) {
       se <- se_kenward(model)$SE
     }
   } else if (!is.null(vcov)) {
-    se <- standard_error(model,
+    se <- standard_error(
+      model,
       vcov = vcov,
       vcov_args = vcov_args,
       component = component,
@@ -106,14 +117,17 @@ p_value_kenward.lmerMod <- function(model, dof = NULL) {
   }
 
   p <- 2 * stats::pt(abs(statistic), df = dof, lower.tail = FALSE)
-  out <- .data_frame(
-    Parameter = params$Parameter,
-    p = unname(p)
-  )
+  out <- .data_frame(Parameter = params$Parameter, p = unname(p))
 
-  if ("Component" %in% names(params)) out$Component <- params$Component
-  if ("Effects" %in% names(params) && effects != "fixed") out$Effects <- params$Effects
-  if ("Response" %in% names(params)) out$Response <- params$Response
+  if ("Component" %in% names(params)) {
+    out$Component <- params$Component
+  }
+  if ("Effects" %in% names(params) && effects != "fixed") {
+    out$Effects <- params$Effects
+  }
+  if ("Response" %in% names(params)) {
+    out$Response <- params$Response
+  }
 
   out
 }
@@ -124,21 +138,28 @@ p_value_kenward.lmerMod <- function(model, dof = NULL) {
     params$SE <- NULL
   }
   params <- merge(params, dof, by = "Parameter")
-  p <- 2 * stats::pt(abs(params$Estimate / params$SE), df = params$df_error, lower.tail = FALSE)
+  p <- 2 *
+    stats::pt(abs(params$Estimate / params$SE), df = params$df_error, lower.tail = FALSE)
 
-  .data_frame(
-    Parameter = params$Parameter,
-    p = unname(p)
-  )
+  .data_frame(Parameter = params$Parameter, p = unname(p))
 }
 
 
 # helper -------------------------
 
-.check_REML_fit <- function(model) {
-  insight::check_if_installed("lme4")
+.check_REML_fit <- function(model, verbose = TRUE) {
+  if (inherits(model, "glmmTMB")) {
+    is_reml <- isTRUE(model$modelInfo$REML)
+  } else {
+    insight::check_if_installed("lme4")
+    is_reml <- lme4::getME(model, "is_REML")
+  }
 
-  if (!(lme4::getME(model, "is_REML"))) {
-    insight::format_warning("Model was not fitted by REML. Re-fitting model now, but p-values, df, etc. still might be unreliable.")
+  if (!is_reml && verbose) {
+    insight::format_warning(
+      "Model was not fitted by REML. Re-fitting model now, but p-values, df, etc. still might be unreliable."
+    )
   }
+
+  is_reml
 }
diff --git a/R/standard_error_kenward.R b/R/standard_error_kenward.R
index 379b3276b..67f53be9a 100644
--- a/R/standard_error_kenward.R
+++ b/R/standard_error_kenward.R
@@ -1,12 +1,35 @@
 #' @rdname p_value_kenward
 #' @export
-se_kenward <- function(model) {
-  .check_REML_fit(model)
-  vcov_adj <- .vcov_kenward_ajusted(model)
+se_kenward <- function(model, ...) {
+  UseMethod("se_kenward")
+}
+
+
+#' @export
+se_kenward.default <- function(model, ...) {
+  if (!.check_REML_fit(model)) {
+    model <- stats::update(model, . ~ ., REML = TRUE)
+  }
+
+  vcov_adjusted <- insight::get_varcov(model, vcov = "kenward-roger")
   params <- insight::get_parameters(model, effects = "fixed")
 
+  .data_frame(Parameter = params$Parameter, SE = abs(sqrt(diag(vcov_adjusted))))
+}
+
+
+#' @export
+se_kenward.glmmTMB <- function(model, component = "conditional", ...) {
+  if (!.check_REML_fit(model)) {
+    model <- stats::update(model, . ~ ., REML = TRUE)
+  }
+
+  vcov_adjusted <- insight::get_varcov(model, vcov = "kenward-roger")
+  params <- insight::get_parameters(model, effects = "fixed", component = component)
+
   .data_frame(
     Parameter = params$Parameter,
-    SE = abs(as.vector(sqrt(diag(as.matrix(vcov_adj)))))
+    SE = abs(sqrt(diag(vcov_adjusted))),
+    Component = component
   )
 }
diff --git a/inst/WORDLIST b/inst/WORDLIST
index 90c8832df..01f2ffa79 100644
--- a/inst/WORDLIST
+++ b/inst/WORDLIST
@@ -291,6 +291,7 @@ labelled
 lakens
 lavaan
 lavaSearch
+lcmm
 lesslikely
 lm
 lme
diff --git a/man/p_value_kenward.Rd b/man/p_value_kenward.Rd
index effcc2037..2d5e230bd 100644
--- a/man/p_value_kenward.Rd
+++ b/man/p_value_kenward.Rd
@@ -8,19 +8,24 @@
 \alias{se_kenward}
 \title{Kenward-Roger approximation for SEs, CIs and p-values}
 \usage{
-ci_kenward(model, ci = 0.95)
+ci_kenward(model, ci = 0.95, ...)
 
 dof_kenward(model)
 
 p_value_kenward(model, dof = NULL)
 
-se_kenward(model)
+se_kenward(model, ...)
 }
 \arguments{
 \item{model}{A statistical model.}
 
 \item{ci}{Confidence Interval (CI) level. Default to \code{0.95} (\verb{95\%}).}
 
+\item{...}{Additional arguments passed down to the underlying functions.
+E.g., arguments like \code{vcov} or \code{vcov_args} can be used to compute confidence
+intervals using a specific variance-covariance matrix for the standard
+errors.}
+
 \item{dof}{Degrees of Freedom.}
 }
 \value{
diff --git a/tests/testthat/test-model_parameters_df_method.R b/tests/testthat/test-model_parameters_df_method.R
index 6513488ab..a3eb2d36a 100644
--- a/tests/testthat/test-model_parameters_df_method.R
+++ b/tests/testthat/test-model_parameters_df_method.R
@@ -1,59 +1,56 @@
 skip_if_not_installed("lmerTest")
 skip_if_not_installed("pbkrtest")
 skip_if_not_installed("lme4")
+skip_if_not_installed("glmmTMB")
 
 mtcars$cyl <- as.factor(mtcars$cyl)
-model <- suppressMessages(lme4::lmer(mpg ~ as.factor(gear) * hp + as.factor(am) + wt + (1 | cyl), data = mtcars))
-model2 <- suppressMessages(lmerTest::lmer(mpg ~ as.factor(gear) * hp + as.factor(am) + wt + (1 | cyl), data = mtcars))
+model <- suppressMessages(lme4::lmer(
+  mpg ~ as.factor(gear) * hp + as.factor(am) + wt + (1 | cyl),
+  data = mtcars
+))
+model2 <- suppressMessages(lmerTest::lmer(
+  mpg ~ as.factor(gear) * hp + as.factor(am) + wt + (1 | cyl),
+  data = mtcars
+))
+model3 <- suppressMessages(glmmTMB::glmmTMB(
+  mpg ~ as.factor(gear) * hp + as.factor(am) + wt + (1 | cyl),
+  data = mtcars,
+  REML = TRUE
+))
+model4 <- suppressMessages(glmmTMB::glmmTMB(
+  mpg ~ as.factor(gear) * hp + as.factor(am) + wt + (1 | cyl),
+  data = mtcars,
+  REML = FALSE
+))
 
 mp0 <- model_parameters(model, digits = 5, effects = "fixed")
 mp1 <- model_parameters(model, digits = 5, ci_method = "normal", effects = "fixed")
 mp2 <- model_parameters(model, digits = 5, ci_method = "s", effects = "fixed")
 mp3 <- model_parameters(model, digits = 5, ci_method = "kr", effects = "fixed")
 mp4 <- model_parameters(model, digits = 5, ci_method = "wald", effects = "fixed")
+mp5 <- model_parameters(model3, digits = 5, ci_method = "kr", effects = "fixed")
+mp6 <- model_parameters(
+  model3,
+  digits = 5,
+  ci_method = "satterthwaite",
+  effects = "fixed"
+)
 
 test_that("model_parameters, ci_method default (residual)", {
   expect_equal(
     mp0$SE,
-    c(
-      2.77457,
-      3.69574,
-      3.521,
-      0.01574,
-      1.58514,
-      0.86316,
-      0.02973,
-      0.01668
-    ),
+    c(2.77457, 3.69574, 3.521, 0.01574, 1.58514, 0.86316, 0.02973, 0.01668),
     tolerance = 1e-3
   )
   expect_equal(mp0$df, c(22, 22, 22, 22, 22, 22, 22, 22), tolerance = 1e-3)
   expect_equal(
     mp0$p,
-    c(
-      0,
-      0.00258,
-      0.14297,
-      0.17095,
-      0.84778,
-      0.00578,
-      0.00151,
-      0.32653
-    ),
+    c(0, 0.00258, 0.14297, 0.17095, 0.84778, 0.00578, 0.00151, 0.32653),
     tolerance = 1e-3
   )
   expect_equal(
     mp0$CI_low,
-    c(
-      24.54722,
-      4.89698,
-      -1.95317,
-      -0.05493,
-      -2.97949,
-      -4.42848,
-      -0.16933,
-      -0.05133
-    ),
+    c(24.54722, 4.89698, -1.95317, -0.05493, -2.97949, -4.42848, -0.16933, -0.05133),
     tolerance = 1e-3
   )
 })
@@ -61,23 +58,10 @@ test_that("model_parameters, ci_method default (residual)", {
 test_that("model_parameters, ci_method normal", {
   expect_equal(
     mp1$SE,
-    c(
-      2.77457,
-      3.69574,
-      3.521,
-      0.01574,
-      1.58514,
-      0.86316,
-      0.02973,
-      0.01668
-    ),
-    tolerance = 1e-3
-  )
-  expect_equal(
-    mp1$df,
-    c(22, 22, 22, 22, 22, 22, 22, 22),
+    c(2.77457, 3.69574, 3.521, 0.01574, 1.58514, 0.86316, 0.02973, 0.01668),
     tolerance = 1e-3
   )
+  expect_equal(mp1$df, c(22, 22, 22, 22, 22, 22, 22, 22), tolerance = 1e-3)
   expect_equal(
     mp1$p,
     c(0, 0.00068, 0.12872, 0.15695, 0.846, 0.00224, 0.00029, 0.31562),
@@ -85,16 +69,7 @@ test_that("model_parameters, ci_method normal", {
   )
   expect_equal(
     mp1$CI_low,
-    c(
-      24.86326,
-      5.31796,
-      -1.5521,
-      -0.05313,
-      -2.79893,
-      -4.33015,
-      -0.16595,
-      -0.04943
-    ),
+    c(24.86326, 5.31796, -1.5521, -0.05313, -2.79893, -4.33015, -0.16595, -0.04943),
     tolerance = 1e-3
   )
 })
@@ -102,150 +77,75 @@ test_that("model_parameters, ci_method normal", {
 test_that("model_parameters, ci_method satterthwaite", {
   expect_equal(
     mp2$SE,
-    c(
-      2.77457,
-      3.69574,
-      3.521,
-      0.01574,
-      1.58514,
-      0.86316,
-      0.02973,
-      0.01668
-    ),
+    c(2.77457, 3.69574, 3.521, 0.01574, 1.58514, 0.86316, 0.02973, 0.01668),
     tolerance = 1e-3
   )
   expect_equal(mp2$df, c(24, 24, 24, 24, 24, 24, 24, 24), tolerance = 1e-3)
   expect_equal(
     mp2$p,
-    c(
-      0,
-      0.00236,
-      0.14179,
-      0.16979,
-      0.84763,
-      0.00542,
-      0.00136,
-      0.32563
-    ),
+    c(0, 0.00236, 0.14179, 0.16979, 0.84763, 0.00542, 0.00136, 0.32563),
     tolerance = 1e-3
   )
   expect_equal(
     mp2$CI_low,
-    c(
-      24.57489,
-      4.93385,
-      -1.91805,
-      -0.05477,
-      -2.96368,
-      -4.41987,
-      -0.16904,
-      -0.05117
-    ),
+    c(24.57489, 4.93385, -1.91805, -0.05477, -2.96368, -4.41987, -0.16904, -0.05117),
     tolerance = 1e-3
   )
+  expect_equal(mp2$SE, mp6$SE, tolerance = 1e-3)
+  expect_equal(mp2$df_error, mp6$df_error, tolerance = 1e-3)
 })
 
 test_that("model_parameters, ci_method kenward", {
   expect_equal(
     mp3$SE,
-    c(
-      2.97608,
-      6.10454,
-      3.98754,
-      0.02032,
-      1.60327,
-      0.91599,
-      0.05509,
-      0.01962
-    ),
+    c(2.97608, 6.10454, 3.98754, 0.02032, 1.60327, 0.91599, 0.05509, 0.01962),
     tolerance = 1e-3
   )
   expect_equal(
     mp3$df,
-    c(
-      19.39553,
-      5.27602,
-      23.57086,
-      8.97297,
-      22.7421,
-      23.76299,
-      2.72622,
-      22.82714
-    ),
+    c(19.39553, 5.27602, 23.57086, 8.97297, 22.7421, 23.76299, 2.72622, 22.82714),
     tolerance = 1e-3
   )
   expect_equal(
     mp3$p,
-    c(
-      0,
-      0.09176,
-      0.19257,
-      0.30147,
-      0.84942,
-      0.00828,
-      0.15478,
-      0.40248
-    ),
+    c(0, 0.09176, 0.19257, 0.30147, 0.84942, 0.00828, 0.15478, 0.40248),
     tolerance = 1e-3
   )
   expect_equal(
     mp3$CI_low,
-    c(
-      24.08091,
-      -2.887,
-      -2.88887,
-      -0.06828,
-      -3.01082,
-      -4.5299,
-      -0.29339,
-      -0.05735
-    ),
+    c(24.08091, -2.887, -2.88887, -0.06828, -3.01082, -4.5299, -0.29339, -0.05735),
     tolerance = 1e-3
   )
+
+  expect_equal(mp5$SE, mp3$SE, tolerance = 1e-3)
+  expect_equal(mp5$df_error, mp3$df_error, tolerance = 1e-3)
+
+  expect_warning(expect_warning(expect_warning(
+    {
+      mp7 <- model_parameters(model4, digits = 5, ci_method = "kr", effects = "fixed")
+    },
+    regex = "Model was not fitted"
+  )))
+
+  expect_equal(mp5$SE, mp7$SE, tolerance = 1e-3)
+  expect_equal(mp5$df_error, mp7$df_error, tolerance = 1e-3)
 })
 
 test_that("model_parameters, ci_method wald (t)", {
   expect_equal(
     mp4$SE,
-    c(
-      2.77457,
-      3.69574,
-      3.521,
-      0.01574,
-      1.58514,
-      0.86316,
-      0.02973,
-      0.01668
-    ),
+    c(2.77457, 3.69574, 3.521, 0.01574, 1.58514, 0.86316, 0.02973, 0.01668),
     tolerance = 1e-3
   )
   expect_equal(mp4$df, c(22, 22, 22, 22, 22, 22, 22, 22), tolerance = 1e-3)
   expect_equal(
     mp4$p,
-    c(
-      0,
-      0.00258,
-      0.14297,
-      0.17095,
-      0.84778,
-      0.00578,
-      0.00151,
-      0.32653
-    ),
+    c(0, 0.00258, 0.14297, 0.17095, 0.84778, 0.00578, 0.00151, 0.32653),
     tolerance = 1e-3
   )
   expect_equal(
     mp4$CI_low,
-    c(
-      24.54722,
-      4.89698,
-      -1.95317,
-      -0.05493,
-      -2.97949,
-      -4.42848,
-      -0.16933,
-      -0.05133
-    ),
+    c(24.54722, 4.89698, -1.95317, -0.05493, -2.97949, -4.42848, -0.16933, -0.05133),
     tolerance = 1e-3
   )
 })
@@ -268,8 +168,10 @@ test_that("model_parameters, satterthwaite Conf Int-1", {
 
 test_that("model_parameters, satterthwaite Conf Int-2", {
   coef.table <- as.data.frame(summary(model2)$coefficients)
-  coef.table$CI_low <- coef.table$Estimate - (coef.table$"Std. Error" * qt(0.975, df = coef.table$df))
-  coef.table$CI_high <- coef.table$Estimate + (coef.table$"Std. Error" * qt(0.975, df = coef.table$df))
+  coef.table$CI_low <- coef.table$Estimate -
+    (coef.table$"Std. Error" * qt(0.975, df = coef.table$df))
+  coef.table$CI_high <- coef.table$Estimate +
+    (coef.table$"Std. Error" * qt(0.975, df = coef.table$df))
 
   expect_equal(mp2$CI_low, coef.table$CI_low, tolerance = 1e-4)
   expect_equal(mp2$CI_high, coef.table$CI_high, tolerance = 1e-4)