diff --git a/DESCRIPTION b/DESCRIPTION
index 7be5b7a22..ba2b05ba1 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Type: Package
 Package: effectsize
 Title: Indices of Effect Size and Standardized Parameters
-Version: 0.5.0.2
+Version: 0.5.0.9
 Authors@R: 
     c(person(given = "Mattan S.",
              family = "Ben-Shachar",
@@ -53,9 +53,9 @@ Depends:
 Imports:
     bayestestR (>= 0.10.5),
     insight (>= 0.14.3),
-    parameters (>= 0.15.0),
+    parameters (>= 0.15.0.1),
     performance (>= 0.8.0),
-    datawizard (>= 0.2.0),
+    datawizard (>= 0.2.1.9001),
     stats,
     utils
 Suggests:
@@ -85,6 +85,9 @@ Suggests:
     spelling,
     testthat,
     tidymodels
+Remotes:
+    easystats/parameters,
+    easystats/datawizard
 VignetteBuilder: 
     knitr
 Encoding: UTF-8
diff --git a/NAMESPACE b/NAMESPACE
index a575560ce..b5c3f1a3f 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -40,19 +40,20 @@ S3method(rb_to_cles,numeric)
 S3method(standardize,Surv)
 S3method(standardize,bcplm)
 S3method(standardize,clm2)
-S3method(standardize,coxme)
-S3method(standardize,coxph)
 S3method(standardize,default)
 S3method(standardize,mediate)
-S3method(standardize,mlm)
 S3method(standardize,wbgee)
 S3method(standardize,wbm)
+S3method(standardize_info,default)
 S3method(standardize_parameters,bootstrap_model)
 S3method(standardize_parameters,bootstrap_parameters)
 S3method(standardize_parameters,default)
 S3method(standardize_parameters,mediate)
 S3method(standardize_parameters,model_fit)
 S3method(standardize_parameters,parameters_model)
+export(.es_aov_simple)
+export(.es_aov_strata)
+export(.es_aov_table)
 export(F_to_d)
 export(F_to_epsilon2)
 export(F_to_eta2)
diff --git a/NEWS.md b/NEWS.md
index cd48be779..1e5a54e3d 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,4 +1,4 @@
-# effectsize 0.5.0.1
+# effectsize 0.5.0.9
 
 ## Breaking Changes
 
@@ -6,6 +6,7 @@
 
 ## New features
 
+- `eta_squared()` family now supports `afex::mixed()` models.
 - `cles()` for estimating common language effect sizes.
 - `rb_to_cles()` for converting rank-biserial correlation to Probability of superiority.
 
@@ -15,7 +16,9 @@
 
 ## Bug fixes
 
-- `eta_squared()`: fixed a bug that caused `afex_aov` models with more than 2 within-subject factros to return incorrect effect sizes for the lower level factors ( #389 ).
+- `standardize()` for multivariate models standardizes the (multivariate) response.
+- `standardize()` for models with offsets standardizes offset variables according to `include_response` and `two_sd` ( #396 ).
+- `eta_squared()`: fixed a bug that caused `afex_aov` models with more than 2 within-subject factors to return incorrect effect sizes for the lower level factors ( #389 ).
 
 # effectsize 0.5.0
 
diff --git a/R/docs_extra.R b/R/docs_extra.R
index 1dfbbd387..bc392b9f4 100644
--- a/R/docs_extra.R
+++ b/R/docs_extra.R
@@ -141,3 +141,12 @@
 #' @rdname effectsize_CIs
 #' @name effectsize_CIs
 NULL
+
+
+#' `effectsize` API
+#'
+#' Read the [*Support functions for model extensions*](https://easystats.github.io/effectsize/articles/effectsize_API.html) vignette.
+#'
+#' @rdname effectsize_API
+#' @name effectsize_API
+NULL
\ No newline at end of file
diff --git a/R/effectsize.BFBayesFactor.R b/R/effectsize.BFBayesFactor.R
index c6ce17fe6..ec41142e5 100644
--- a/R/effectsize.BFBayesFactor.R
+++ b/R/effectsize.BFBayesFactor.R
@@ -55,7 +55,7 @@ effectsize.BFBayesFactor <- function(model, type = NULL, verbose = TRUE, test =
       mu <- 0
       D <- samps$delta
     } else {
-      mu <- as.numeric(gsub("Null, mu=", "", model@denominator@shortName))
+      mu <- model@numerator[[1]]@prior$mu
       D <- (samps$mu - mu) / sqrt(samps$sig2)
     }
 
diff --git a/R/effectsize.htest.R b/R/effectsize.htest.R
index d6bb78288..f3432b049 100644
--- a/R/effectsize.htest.R
+++ b/R/effectsize.htest.R
@@ -115,7 +115,7 @@ effectsize.htest <- function(model, type = NULL, verbose = TRUE, ...) {
     return(out)
   } else if (grepl("One-way", model$method)) {
     # one way anove ----
-    if (approx <- grepl("not assuming", model$method, fixed = TRUE) && verbose) {
+    if ((approx <- grepl("not assuming", model$method, fixed = TRUE)) && verbose) {
       warning("`var.equal = FALSE` - effect size is an approximation.", call. = FALSE)
     }
 
diff --git a/R/eta_squared.R b/R/eta_squared.R
index ffdd5411b..503c59a2d 100644
--- a/R/eta_squared.R
+++ b/R/eta_squared.R
@@ -361,100 +361,50 @@ cohens_f_squared <- function(model, partial = TRUE, ci = 0.95, alternative = "gr
 
 # Get ES ------------------------------------------------------------------
 
-#' @keywords internal
-.anova_es <-
-  function(model,
-           type = c("eta", "omega", "epsilon"),
-           partial = TRUE,
-           generalized = FALSE,
-           ci = 0.95, alternative = "greater",
-           verbose = TRUE,
-           ...) {
-    UseMethod(".anova_es")
-  }
-
-#' @keywords internal
-#' @importFrom stats anova
-.anova_es.default <- function(model,
-                              type = c("eta", "omega", "epsilon"),
-                              partial = TRUE,
-                              generalized = FALSE,
-                              ci = 0.95, alternative = "greater",
-                              verbose = TRUE,
-                              ...) {
-  .anova_es.anova(
-    stats::anova(model),
-    type = type,
-    partial = partial,
-    generalized = generalized,
-    ci = ci,
-    alternative = alternative,
-    verbose = verbose
-  )
-}
-
-#' @keywords internal
-#' @importFrom parameters model_parameters
-#' @importFrom stats anova
-.anova_es.aov <- function(model,
-                          type = c("eta", "omega", "epsilon"),
-                          partial = TRUE,
-                          generalized = FALSE,
-                          ci = 0.95, alternative = "greater",
-                          verbose = TRUE,
-                          ...) {
-  if (!inherits(model, c("Gam", "anova"))) {
-    # Pass to ANOVA table method
-    res <- .anova_es.anova(
-      stats::anova(model),
-      type = type,
-      partial = partial,
-      generalized = generalized,
-      ci = ci, alternative = alternative,
-      verbose = verbose,
-      ...
-    )
-    return(res)
-  }
-
+#' @param aov_table Input data frame
+#' @param type Which effect size to compute?
+#' @param include_intercept Should the intercept (`(Intercept)`) be included?
+#' @param partial,generalized,ci,alternative,verbose See [eta_squared()].
+#'
+#' @rdname effectsize_API
+#' @export
+.es_aov_simple <- function(aov_table,
+                           type = c("eta", "omega", "epsilon"),
+                           partial = TRUE,
+                           generalized = FALSE,
+                           ci = 0.95, alternative = "greater",
+                           verbose = TRUE,
+                           include_intercept = FALSE) {
   type <- match.arg(type)
+  aov_table <- as.data.frame(aov_table)
 
-  params <- parameters::model_parameters(model, verbose = verbose, effects = "fixed")
-  out <- .es_aov(as.data.frame(params), type, partial, generalized, ci, alternative, verbose = verbose, ...)
-  if (is.null(attr(out, "anova_type"))) attr(out, "anova_type") <- attr(params, "anova_type")
-  out
-}
+  # Clean up data ---
+  if (!"Mean_Square" %in% colnames(aov_table)) {
+    aov_table[["Mean_Square"]] <- aov_table[["Sum_Squares"]] / aov_table[["df"]]
+  }
 
-#' @keywords internal
-.es_aov <- function(params,
-                    type = c("eta", "omega", "epsilon"),
-                    partial = TRUE,
-                    generalized = FALSE,
-                    ci = 0.95, alternative = "greater",
-                    verbose = TRUE,
-                    include_intercept = FALSE,
-                    ...) {
-
-  if (!"Residuals" %in% params$Parameter) {
-    stop(
-      "No residuals data found - ",
-      type,
-      " squared can only be computed for simple `aov` models."
-    )
+  if (!"Residuals" %in% aov_table$Parameter) {
+    stop(insight::format_message("No residuals data found - cannot compute effect size."))
   }
 
+
+  # Include intercept? ---
   if (include_intercept) {
-    values <- .values_aov(params[params$Parameter != "(Intercept)", ])
+    values <- .values_aov(aov_table[aov_table$Parameter != "(Intercept)", ])
   } else {
-    params <- params[params$Parameter != "(Intercept)", ]
-    values <- .values_aov(params)
+    aov_table <- aov_table[aov_table$Parameter != "(Intercept)", ]
+    values <- .values_aov(aov_table)
   }
 
-  df_error <- params$df[params$Parameter == "Residuals"]
-  params <- params[params$Parameter != "Residuals", , drop = FALSE]
+  # Get error df ---
+  df_error <- aov_table$df[aov_table$Parameter == "Residuals"]
+  aov_table <- aov_table[aov_table$Parameter != "Residuals", , drop = FALSE]
 
-  if (nrow(params) == 1L &&
-    (partial || isTRUE(generalized) || is.character(generalized))) {
+
+  # Validate anova type (1,2,3) and partial ---
+  anova_type <- NULL
+  if (nrow(aov_table) == 1L &&
+      (partial || isTRUE(generalized) || is.character(generalized))) {
     if (verbose) {
       txt_type <- ifelse(isTRUE(generalized) || is.character(generalized), "generalized", "partial")
       message(
@@ -464,91 +414,90 @@ cohens_f_squared <- function(model, partial = TRUE, ci = 0.95, alternative = "gr
     }
     partial <- FALSE
     anova_type <- NA
-  } else {
-    anova_type <- NULL
   }
 
 
-
+  # Estimate effect size ---
   if (type == "eta") {
     if (isTRUE(generalized) || is.character(generalized)) {
       ## copied from afex
-      obs <- logical(nrow(params))
+      obs <- logical(nrow(aov_table))
       if (is.character(generalized)) {
         for (o in generalized) {
-          oi <- grepl(paste0("\\b", o, "\\b"), params$Parameter)
+          oi <- grepl(paste0("\\b", o, "\\b"), aov_table$Parameter)
 
           if (!any(oi)) stop("Observed variable not in data: ", o, call. = FALSE)
 
           obs <- obs | oi
         }
       }
-      obs_SSn1 <- sum(params$Sum_Squares * obs)
-      obs_SSn2 <- params$Sum_Squares * obs
+      obs_SSn1 <- sum(aov_table$Sum_Squares * obs)
+      obs_SSn2 <- aov_table$Sum_Squares * obs
 
-      params$Eta2_generalized <- params$Sum_Squares /
-        (params$Sum_Squares + values$Sum_Squares_residuals + obs_SSn1 - obs_SSn2)
+      aov_table$Eta2_generalized <- aov_table$Sum_Squares /
+        (aov_table$Sum_Squares + values$Sum_Squares_residuals + obs_SSn1 - obs_SSn2)
     } else if (!isTRUE(partial)) {
-      params$Eta2 <- params$Sum_Squares /
+      aov_table$Eta2 <- aov_table$Sum_Squares /
         values$Sum_Squares_total
     } else {
-      params$Eta2_partial <-
-        params$Sum_Squares /
-          (params$Sum_Squares + values$Sum_Squares_residuals)
+      aov_table$Eta2_partial <-
+        aov_table$Sum_Squares /
+        (aov_table$Sum_Squares + values$Sum_Squares_residuals)
     }
   } else if (type == "omega") {
     if (!isTRUE(partial)) {
-      params$Omega2 <-
-        (params$Sum_Squares - params$df * values$Mean_Square_residuals) /
-          (values$Sum_Squares_total + values$Mean_Square_residuals)
+      aov_table$Omega2 <-
+        (aov_table$Sum_Squares - aov_table$df * values$Mean_Square_residuals) /
+        (values$Sum_Squares_total + values$Mean_Square_residuals)
     } else {
-      params$Omega2_partial <-
-        (params$Sum_Squares - params$df * values$Mean_Square_residuals) /
-          (params$Sum_Squares + (values$n - params$df) * values$Mean_Square_residuals)
+      aov_table$Omega2_partial <-
+        (aov_table$Sum_Squares - aov_table$df * values$Mean_Square_residuals) /
+        (aov_table$Sum_Squares + (values$n - aov_table$df) * values$Mean_Square_residuals)
     }
   } else if (type == "epsilon") {
     if (!isTRUE(partial)) {
-      params$Epsilon2 <-
-        (params$Sum_Squares - params$df * values$Mean_Square_residuals) /
-          values$Sum_Squares_total
+      aov_table$Epsilon2 <-
+        (aov_table$Sum_Squares - aov_table$df * values$Mean_Square_residuals) /
+        values$Sum_Squares_total
     } else {
-      params$Epsilon2_partial <-
-        (params$Sum_Squares - params$df * values$Mean_Square_residuals) /
-          (params$Sum_Squares + values$Sum_Squares_residuals)
+      aov_table$Epsilon2_partial <-
+        (aov_table$Sum_Squares - aov_table$df * values$Mean_Square_residuals) /
+        (aov_table$Sum_Squares + values$Sum_Squares_residuals)
     }
   }
 
+  out <- aov_table
 
-
-  out <- params
-
+  # Add CIs ---
   if (is.numeric(ci)) {
     # based on MBESS::ci.R2
     ES <- pmax(0, out[[ncol(out)]])
     f <- (ES / out$df) / ((1 - ES) / df_error)
 
-    out <- cbind(
-      out,
-      # This really is a generic F_to_R2
+    CI_tab <- # This really is a generic F_to_R2
       F_to_eta2(f,
-        out$df,
-        df_error,
-        ci = ci, alternative = alternative
+                out$df,
+                df_error,
+                ci = ci, alternative = alternative
       )[-1]
-    )
+
+    out[c("CI", "CI_low", "CI_high")] <- CI_tab[c("CI", "CI_low", "CI_high")]
   } else {
     alternative <- NULL
   }
 
 
+  # Clean up output ---
   out <- out[, colnames(out) %in% c(
-    "Group", "Response", "Parameter",
+    "Parameter",
     "Eta2", "Eta2_partial", "Eta2_generalized",
     "Omega2", "Omega2_partial",
     "Epsilon2", "Epsilon2_partial",
-    "CI", "CI_low", "CI_high"
+    if (!is.null(ci)) c("CI", "CI_low", "CI_high")
   ), drop = FALSE]
+  rownames(out) <- NULL
 
+  # Set attributes ---
   attr(out, "partial") <- partial
   attr(out, "generalized") <- generalized
   attr(out, "ci") <- ci
@@ -558,166 +507,135 @@ cohens_f_squared <- function(model, partial = TRUE, ci = 0.95, alternative = "gr
   out
 }
 
-.anova_es.lm <- .anova_es.aov
-
-.anova_es.glm <- .anova_es.aov
-
-.anova_es.manova <- .anova_es.aov
-
-#' @keywords internal
-#' @importFrom parameters model_parameters
-#' @importFrom insight find_predictors
-.anova_es.aovlist <- function(model,
-                              type = c("eta", "omega", "epsilon"),
-                              partial = TRUE,
-                              generalized = FALSE,
-                              ci = 0.95, alternative = "greater",
-                              verbose = TRUE,
-                              include_intercept = FALSE,
-                              ...) {
-  params <- parameters::model_parameters(model, verbose = verbose, effects = "fixed")
-  anova_type <- attr(params, "anova_type")
-  params <- as.data.frame(params)
-
-  IVs <- insight::find_predictors(model)[[1]]
-
-  out <-
-    .es_aovlist(
-      params,
-      IVs = IVs,
-      type = type,
-      partial = partial,
-      generalized = generalized,
-      ci = ci, alternative = alternative,
-      verbose = verbose,
-      include_intercept = include_intercept
-    )
-  attr(out, "anova_type") <- anova_type
-  out
-}
-
-#' @keywords internal
-.es_aovlist <- function(params, IVs,
-                        type = c("eta", "omega", "epsilon"),
-                        partial = TRUE,
-                        generalized = FALSE,
-                        ci = 0.95, alternative = "greater",
-                        verbose = TRUE,
-                        include_intercept = FALSE) {
+#' @param DV_names A character vector with the names of all the predictors,
+#'   including the grouping variable (e.g., `"Subject"`).
+#'
+#' @rdname effectsize_API
+#' @export
+.es_aov_strata <- function(aov_table, DV_names,
+                           type = c("eta", "omega", "epsilon"),
+                           partial = TRUE,
+                           generalized = FALSE,
+                           ci = 0.95, alternative = "greater",
+                           verbose = TRUE,
+                           include_intercept = FALSE) {
   type <- match.arg(type)
+  aov_table <- as.data.frame(aov_table)
 
-  if (!"Residuals" %in% params$Parameter) {
-    stop(
-      "No residuals data found - ",
-      type,
-      " squared can only be computed for simple `aov` models."
-    )
+  # Clean up data ---
+  if (!"Mean_Square" %in% colnames(aov_table)) {
+    aov_table[["Mean_Square"]] <- aov_table[["Sum_Squares"]] / aov_table[["df"]]
   }
 
+  if (!"Residuals" %in% aov_table$Parameter) {
+    stop(insight::format_message("No residuals data found - cannot compute effect size."))
+  }
+
+
+  # Include intercept? ---
   if (include_intercept) {
-    values <- .values_aov(params[params$Parameter != "(Intercept)", ])
+    values <- .values_aov(aov_table[aov_table$Parameter != "(Intercept)", ], group = TRUE)
   } else {
-    params <- params[params$Parameter != "(Intercept)", ]
-    values <- .values_aov(params)
+    aov_table <- aov_table[aov_table$Parameter != "(Intercept)", ]
+    values <- .values_aov(aov_table, group = TRUE)
   }
 
-  params <- params[params$Parameter != "Residuals" & !is.na(params$`F`), , drop = FALSE]
+
+  # Get all the correct SSs... ---
+  aov_table <- aov_table[aov_table$Parameter != "Residuals", , drop = FALSE]
   Sum_Squares_total <- sum(sapply(values, "[[", "Sum_Squares_total"))
-  Sum_Squares_residuals <- sapply(values[params$Group], "[[", "Sum_Squares_residuals")
-  Mean_Square_residuals <- sapply(values[params$Group], "[[", "Mean_Square_residuals")
-  df_residuals <- sapply(values[params$Group], "[[", "df_residuals")
-  ns <- sapply(values[params$Group], "[[", "n")
+  Sum_Squares_residuals <- sapply(values[aov_table$Group], "[[", "Sum_Squares_residuals")
+  Mean_Square_residuals <- sapply(values[aov_table$Group], "[[", "Mean_Square_residuals")
+  df_residuals <- sapply(values[aov_table$Group], "[[", "df_residuals")
+  ns <- sapply(values[aov_table$Group], "[[", "n")
 
 
+  # Estimate effect size ---
   if (type == "eta") {
     if (isTRUE(generalized) || is.character(generalized)) {
       ## copied from afex
-      obs <- logical(nrow(params))
+      obs <- logical(nrow(aov_table))
       if (is.character(generalized)) {
         for (o in generalized) {
-          oi <- grepl(paste0("\\b", o, "\\b"), params$Parameter)
+          oi <- grepl(paste0("\\b", o, "\\b"), aov_table$Parameter)
 
           if (!any(oi)) stop("Observed variable not in data: ", o, call. = FALSE)
 
           obs <- obs | oi
         }
       }
-      obs_SSn1 <- sum(params$Sum_Squares * obs)
-      obs_SSn2 <- params$Sum_Squares * obs
+      obs_SSn1 <- sum(aov_table$Sum_Squares * obs)
+      obs_SSn2 <- aov_table$Sum_Squares * obs
 
-      params$Eta2_generalized <- params$Sum_Squares /
-        (params$Sum_Squares + sum(sapply(values, "[[", "Sum_Squares_residuals")) +
+      aov_table$Eta2_generalized <- aov_table$Sum_Squares /
+        (aov_table$Sum_Squares + sum(sapply(values, "[[", "Sum_Squares_residuals")) +
            obs_SSn1 - obs_SSn2)
     } else if (!isTRUE(partial)) {
-      params$Eta2 <- params$Sum_Squares / Sum_Squares_total
+      aov_table$Eta2 <- aov_table$Sum_Squares / Sum_Squares_total
     } else {
-      params$Eta2_partial <-
-        params$Sum_Squares /
-        (params$Sum_Squares + Sum_Squares_residuals)
+      aov_table$Eta2_partial <-
+        aov_table$Sum_Squares /
+        (aov_table$Sum_Squares + Sum_Squares_residuals)
     }
   } else if (type == "omega") {
-    SSS_values <- values[[which(names(values) %in% IVs)]]
-    is_within <- !params$Group %in% IVs
+    SSS_values <- values[[which(names(values) %in% DV_names)]]
+    is_within <- !aov_table$Group %in% DV_names
     Sum_Squares_Subjects <- SSS_values$Sum_Squares_residuals
     Mean_Squares_Subjects <- SSS_values$Mean_Square_residuals
 
     # implemented from https://www.jasonfinley.com/tools/OmegaSquaredQuickRef_JRF_3-31-13.pdf
     if (!isTRUE(partial)) {
-      params$Omega2 <-
-        (params$Sum_Squares - params$df * Mean_Square_residuals) /
+      aov_table$Omega2 <-
+        (aov_table$Sum_Squares - aov_table$df * Mean_Square_residuals) /
         (Sum_Squares_total + Mean_Squares_Subjects)
     } else {
-      params$Omega2_partial <-
-        (params$Sum_Squares - params$df * Mean_Square_residuals) /
-        (params$Sum_Squares + is_within * Sum_Squares_residuals +
+      aov_table$Omega2_partial <-
+        (aov_table$Sum_Squares - aov_table$df * Mean_Square_residuals) /
+        (aov_table$Sum_Squares + is_within * Sum_Squares_residuals +
            Sum_Squares_Subjects + Mean_Squares_Subjects)
     }
   } else if (type == "epsilon") {
     if (!isTRUE(partial)) {
-      params$Epsilon2 <-
-        (params$Sum_Squares - params$df * Mean_Square_residuals) /
+      aov_table$Epsilon2 <-
+        (aov_table$Sum_Squares - aov_table$df * Mean_Square_residuals) /
         Sum_Squares_total
     } else {
-      params$Epsilon2_partial <-
-        (params$Sum_Squares - params$df * Mean_Square_residuals) /
-        (params$Sum_Squares + Sum_Squares_residuals)
+      aov_table$Epsilon2_partial <-
+        (aov_table$Sum_Squares - aov_table$df * Mean_Square_residuals) /
+        (aov_table$Sum_Squares + Sum_Squares_residuals)
     }
   }
 
-  out <- params
+  out <- aov_table
+
 
+  # Add CIs ---
   if (!is.null(ci)) {
     # based on MBESS::ci.R2
     ES <- pmax(0, out[[ncol(out)]])
     f <- (ES / out$df) / ((1 - ES) / df_residuals)
 
-    out <- cbind(
-      out,
-      # This really is a generic F_to_R2
+    CI_tab <- # This really is a generic F_to_R2
       F_to_eta2(f,
                 out$df,
                 df_residuals,
                 ci = ci, alternative = alternative
       )[-1]
-    )
+
+    out[c("CI", "CI_low", "CI_high")] <- CI_tab[c("CI", "CI_low", "CI_high")]
   } else {
     alternative <- NULL
   }
 
+
+  # Clean up output ---
   out <- out[, colnames(out) %in% c(
     "Group",
-    "Response",
     "Parameter",
-    "Eta2_generalized",
-    "Eta2_partial",
-    "Omega2_partial",
-    "Epsilon2_partial",
-    "Eta2",
-    "Omega2",
-    "Epsilon2",
-    "CI",
-    "CI_low",
-    "CI_high"
+    "Eta2", "Eta2_generalized", "Eta2_partial",
+    "Omega2", "Omega2_partial",
+    "Epsilon2", "Epsilon2_partial",
+    if (!is.null(ci)) c("CI", "CI_low", "CI_high")
   ), drop = FALSE]
   rownames(out) <- NULL
 
@@ -729,6 +647,180 @@ cohens_f_squared <- function(model, partial = TRUE, ci = 0.95, alternative = "gr
   out
 }
 
+#' @rdname effectsize_API
+#' @export
+.es_aov_table <- function(aov_table,
+                          type = c("eta", "omega", "epsilon"),
+                          partial = TRUE,
+                          generalized = FALSE,
+                          ci = 0.95, alternative = "greater",
+                          verbose = TRUE,
+                          include_intercept = FALSE) {
+  aov_table <- as.data.frame(aov_table)
+
+  # Get correct function ---
+  type <- match.arg(type)
+  es_fun <- switch(type,
+                   eta = F_to_eta2,
+                   omega = F_to_omega2,
+                   epsilon = F_to_epsilon2)
+
+
+  # Non-Partial / Generalized -> BAD ---
+  if (verbose) {
+    if (!isTRUE(partial))
+      warning(
+        "Currently only supports partial ",
+        type,
+        " squared for this class of objects.",
+        call. = FALSE
+      )
+
+    if (isTRUE(generalized) || is.character(generalized))
+      warning(
+        "generalized ", type, " squared ",
+        "is not supported for this class of object.",
+        call. = FALSE
+      )
+  }
+
+
+  # Turn ts to Fs (if needed) ---
+  if (!"F" %in% colnames(aov_table))
+    if ("t" %in% colnames(aov_table)) {
+      aov_table[["F"]] <- aov_table[["t"]]^2
+      aov_table[["df"]] <- 1
+    } else {
+      stop(insight::format_message("ANOVA table does not have F values - cannot compute effect size."))
+    }
+
+
+  # include_intercept? ---
+  if (!include_intercept)
+    aov_table <- aov_table[aov_table$Parameter != "(Intercept)", , drop = FALSE]
+
+  ES_tab <- es_fun(aov_table[["F"]],
+                   aov_table[["df"]],
+                   aov_table[["df_error"]],
+                   ci = ci, alternative = alternative)
+
+  out <- cbind(Parameter = aov_table[["Parameter"]], ES_tab)
+  rownames(out) <- NULL
+
+  # Set attributes ---
+  attr(out, "partial") <- TRUE
+  attr(out, "generalized") <- FALSE
+  attr(out, "ci") <- if ("CI" %in% colnames(out)) ci
+  attr(out, "alternative") <- if (!is.null(attr(out, "ci"))) alternative
+  attr(out, "anova_type") <- NULL
+  attr(out, "approximate") <- NULL
+  out
+}
+
+# Clean up wrappers -------------------------------------------------------
+
+
+
+#' @keywords internal
+.anova_es <-
+  function(model,
+           type = c("eta", "omega", "epsilon"),
+           partial = TRUE,
+           generalized = FALSE,
+           ci = 0.95, alternative = "greater",
+           verbose = TRUE,
+           ...) {
+    UseMethod(".anova_es")
+  }
+
+#' @keywords internal
+#' @importFrom stats anova
+.anova_es.default <- function(model,
+                              type = c("eta", "omega", "epsilon"),
+                              partial = TRUE,
+                              generalized = FALSE,
+                              ci = 0.95, alternative = "greater",
+                              verbose = TRUE,
+                              ...) {
+  .anova_es.anova(
+    stats::anova(model),
+    type = type,
+    partial = partial,
+    generalized = generalized,
+    ci = ci,
+    alternative = alternative,
+    verbose = verbose
+  )
+}
+
+#' @keywords internal
+#' @importFrom parameters model_parameters
+#' @importFrom stats anova
+.anova_es.aov <- function(model,
+                          type = c("eta", "omega", "epsilon"),
+                          partial = TRUE,
+                          generalized = FALSE,
+                          ci = 0.95, alternative = "greater",
+                          verbose = TRUE,
+                          ...) {
+  if (!inherits(model, c("Gam", "anova"))) {
+    # Pass to ANOVA table method
+    res <- .anova_es.anova(
+      stats::anova(model),
+      type = type,
+      partial = partial,
+      generalized = generalized,
+      ci = ci, alternative = alternative,
+      verbose = verbose,
+      ...
+    )
+    return(res)
+  }
+
+  params <- parameters::model_parameters(model, verbose = verbose, effects = "fixed")
+  out <- .es_aov_simple(as.data.frame(params), type, partial, generalized, ci, alternative, verbose = verbose, ...)
+  if (is.null(attr(out, "anova_type"))) attr(out, "anova_type") <- attr(params, "anova_type")
+  out
+}
+
+.anova_es.lm <- .anova_es.aov
+
+.anova_es.glm <- .anova_es.aov
+
+.anova_es.manova <- .anova_es.aov
+
+#' @keywords internal
+#' @importFrom parameters model_parameters
+#' @importFrom insight find_predictors
+.anova_es.aovlist <- function(model,
+                              type = c("eta", "omega", "epsilon"),
+                              partial = TRUE,
+                              generalized = FALSE,
+                              ci = 0.95, alternative = "greater",
+                              verbose = TRUE,
+                              include_intercept = FALSE,
+                              ...) {
+  params <- parameters::model_parameters(model, verbose = verbose, effects = "fixed")
+  anova_type <- attr(params, "anova_type")
+  params <- as.data.frame(params)
+
+  DV_names <- insight::find_predictors(model)[[1]]
+
+  out <-
+    .es_aov_strata(
+      params,
+      DV_names = DV_names,
+      type = type,
+      partial = partial,
+      generalized = generalized,
+      ci = ci, alternative = alternative,
+      verbose = verbose,
+      include_intercept = include_intercept
+    )
+  attr(out, "anova_type") <- anova_type
+  out
+}
+
 #' @keywords internal
 #' @importFrom stats aov
 #' @importFrom utils packageVersion
@@ -743,9 +835,8 @@ cohens_f_squared <- function(model, partial = TRUE, ci = 0.95, alternative = "gr
   params <- parameters::model_parameters(model, verbose = verbose, effects = "fixed")
   anova_type <- attr(params, "anova_type")
 
-  params <- as.data.frame(params)
   params <- split(params, params$Response)
-  params <- lapply(params, .es_aov,
+  params <- lapply(params, .es_aov_simple,
     type = type,
     partial = partial,
     generalized = generalized,
@@ -753,6 +844,10 @@ cohens_f_squared <- function(model, partial = TRUE, ci = 0.95, alternative = "gr
     verbose = verbose,
     ...
   )
+
+  params <- lapply(names(params), function(nm) {
+    cbind(Response = nm, params[[nm]])
+  })
   out <- do.call("rbind", params)
   rownames(out) <- NULL
 
@@ -776,18 +871,12 @@ cohens_f_squared <- function(model, partial = TRUE, ci = 0.95, alternative = "gr
                             verbose = TRUE,
                             include_intercept = FALSE,
                             ...) {
-  type <- match.arg(type)
-  es_fun <- switch(type,
-    eta = F_to_eta2,
-    omega = F_to_omega2,
-    epsilon = F_to_epsilon2
-  )
-
-  F_val <- c("F value", "approx F", "F-value")
-  numDF <- c("NumDF", "num Df", "numDF", "npar")
-  denDF <- c("DenDF", "den Df", "denDF", "df_error")
+  F.nm <- c("F value", "approx F", "F-value")
+  df.nm <- c("NumDF", "num Df", "numDF", "npar")
+  df_error.nm <- c("DenDF", "den Df", "denDF", "df_error")
 
-  if (!any(denDF %in% colnames(model))) {
+  # If there is no df_error *or* is there IS a residuals row...
+  if (!any(df_error.nm %in% colnames(model))) {
     # Pass to AOV method
     res <- .anova_es.aov(model,
       partial = partial,
@@ -801,51 +890,30 @@ cohens_f_squared <- function(model, partial = TRUE, ci = 0.95, alternative = "gr
     return(res)
   }
 
-  anova_type <- tryCatch(attr(parameters::model_parameters(model, verbose = FALSE, effects = "fixed"), "anova_type"),
-                         error = function(...) 1)
-
-  if (!include_intercept) model <- model[rownames(model) != "(Intercept)", , drop = FALSE]
-  model <- model[rownames(model) != "Residuals", , drop = FALSE]
-
-  F_val <- F_val[F_val %in% colnames(model)]
-  numDF <- numDF[numDF %in% colnames(model)]
-  denDF <- denDF[denDF %in% colnames(model)]
-
-
-  if (verbose && !isTRUE(partial)) {
-    warning(
-      "Currently only supports partial ",
-      type,
-      " squared for this class of objects.",
-      call. = FALSE
-    )
-    partial <- TRUE
-  }
-
-  if (verbose && (isTRUE(generalized) || is.character(generalized))) {
-    warning(
-      "generalized ", type, " squared ",
-      "is not supported for this class of object."
-    )
-    generalized <- FALSE
-  }
-
-  par_table <- as.data.frame(model)
+  # Clean up table ---
+  par_table <- data.frame(
+    Parameter = rownames(model),
+    F = model[,F.nm[F.nm %in% colnames(model)]],
+    df = model[,df.nm[df.nm %in% colnames(model)]],
+    df_error = model[,df_error.nm[df_error.nm %in% colnames(model)]]
+  )
+  par_table <- par_table[!par_table[["Parameter"]] %in% "Residuals",]
 
-  out <- cbind(
-    Parameter = rownames(par_table),
-    es_fun(par_table[[F_val]],
-      par_table[[numDF]],
-      par_table[[denDF]],
-      ci = ci, alternative = alternative
+  out <-
+    .es_aov_table(
+      par_table,
+      type = type,
+      partial = partial,
+      generalized = generalized,
+      ci = ci,
+      alternative = alternative,
+      verbose = verbose,
+      include_intercept = include_intercept
     )
-  )
 
-  attr(out, "partial") <- partial
-  attr(out, "generalized") <- generalized
-  attr(out, "ci") <- ci
-  attr(out, "anova_type") <- anova_type
-  attr(out, "alternative") <- if (is.numeric(ci)) alternative
+  attr(out, "anova_type") <- tryCatch(attr(parameters::model_parameters(model, verbose = FALSE, effects = "fixed"), "anova_type"),
+                                      error = function(...) 1)
+  attr(out, "approximate") <- TRUE
   out
 }
 
@@ -860,71 +928,62 @@ cohens_f_squared <- function(model, partial = TRUE, ci = 0.95, alternative = "gr
                                        generalized = FALSE,
                                        ci = 0.95, alternative = "greater",
                                        verbose = TRUE,
+                                       by_response = TRUE,
                                        ...) {
-  if (verbose && !isTRUE(partial)) {
-    warning(
-      "Currently only supports partial ",
-      type,
-      " squared for this class of objects.",
-      call. = FALSE
-    )
-    partial <- TRUE
+  if (by_response && "Response" %in% colnames(model)) {
+    out <- split(model, model[["Response"]])
+    out <- lapply(out, .anova_es.parameters_model,
+                  type = type, partial = partial, generalized = generalized,
+                  ci = ci, alternative = alternative,
+                  verbose = verbose,
+                  by_response = FALSE,
+                  ...)
+    saved_attr <- attributes(out[[1]])
+    out <- mapply(out, names(out),
+                  FUN = function(x, nm) cbind(Response = nm, x),
+                  SIMPLIFY = FALSE)
+    out <- do.call(rbind, out)
+
+    # Set attributes ---
+    attr(out, "partial") <- saved_attr$partial
+    attr(out, "generalized") <- saved_attr$generalized
+    attr(out, "ci") <- saved_attr$ci
+    attr(out, "alternative") <- saved_attr$alternative
+    attr(out, "anova_type") <- attr(model, "anova_type")
+    attr(out, "approximate") <- saved_attr$approximate
+    return(out)
   }
 
-  if (verbose && (isTRUE(generalized) || is.character(generalized))) {
-    warning(
-      "generalized ", type, " squared ",
-      "is not supported for this class of object."
-    )
-    generalized <- FALSE
-  }
-
-
-  type <- match.arg(type)
 
-  if ("Group" %in% colnames(model) && sum(model$Parameter == "Residuals") > 1) {
-    x <- split(model, model$Group)
-    out <- do.call(rbind, lapply(x, function(i) {
-      f <- i[["F"]]
-      df_num <- i[["df"]][!is.na(f)]
-      df_error <- i[i$Parameter == "Residuals", "df"]
-      cbind(
-        data.frame(Group = unique(i$Group), stringsAsFactors = FALSE),
-        .anova_es_model_params(i, f, df_num, df_error, type, ci, alternative)
+  approximate <- FALSE
+  if ("Sum_Squares" %in% colnames(model) && "Residuals" %in% model[["Parameter"]]) {
+    if ("Group" %in% colnames(model)) {
+      DVs <- unlist(insight::find_predictors(.get_object(model)))
+      out <- .es_aov_strata(
+        model, DV_names = DVs,
+        type = type, partial = partial, generalized = generalized,
+        ci = ci, alternative = alternative,
+        verbose = verbose, ...
       )
-    }))
-  } else {
-    if ("t" %in% colnames(model)) {
-      f <- model[["t"]]^2
-    }
-    if ("F" %in% colnames(model)) {
-      f <- model[["F"]]
-    }
-    if ("z" %in% colnames(model)) {
-      stop("Cannot compute effect size from models with no proper residual variance. Consider the estimates themselves as indices of effect size.")
-    }
-
-    df_col <- colnames(model)[colnames(model) %in% c("df", "Df", "NumDF")]
-    if (length(df_col)) {
-      df_num <- model[[df_col]][!is.na(f)]
     } else {
-      df_num <- 1
-    }
-
-    if ("df_error" %in% colnames(model)) {
-      df_error <- model$df_error
-    } else if ("Residuals" %in% model$Parameter) {
-      df_error <- model[model$Parameter == "Residuals", df_col]
-    } else {
-      stop("Cannot extract degrees of freedom for the error term. Try passing the model object directly to 'eta_squared()'.")
+      out <- .es_aov_simple(
+        model,
+        type = type, partial = partial, generalized = generalized,
+        ci = ci, alternative = alternative,
+        verbose = verbose, ...
+      )
     }
-    out <- .anova_es_model_params(model, f, df_num, df_error, type, ci, alternative)
+  } else {
+    out <- .es_aov_table(
+      model,
+      type = type, partial = partial, generalized = generalized,
+      ci = ci, alternative = alternative,
+      verbose = verbose, ...
+    )
+    approximate <- TRUE
   }
-
-  attr(out, "partial") <- partial
-  attr(out, "generalized") <- generalized
-  attr(out, "ci") <- ci
-  attr(out, "alternative") <- if (is.numeric(ci)) alternative
+  attr(out, "anova_type") <- attr(model, "anova_type")
+  attr(out, "approximate") <- approximate
   out
 }
 
@@ -1003,7 +1062,18 @@ cohens_f_squared <- function(model, partial = TRUE, ci = 0.95, alternative = "gr
 
   model <- lmerTest::as_lmerModLmerTest(model)
   model <- stats::anova(model)
-  .anova_es.anova(model, type = type, partial = partial, generalized = generalized, ci = ci, alternative = alternative, ...)
+  out <-
+    .anova_es.anova(
+      model,
+      type = type,
+      partial = partial,
+      generalized = generalized,
+      ci = ci,
+      alternative = alternative,
+      ...
+    )
+  attr(out, "approximate") <- TRUE
+  out
 }
 
 #' @keywords internal
@@ -1015,31 +1085,6 @@ cohens_f_squared <- function(model, partial = TRUE, ci = 0.95, alternative = "gr
                           ci = 0.95, alternative = "greater",
                           verbose = TRUE,
                           ...) {
-  type <- match.arg(type)
-  es_fun <- switch(type,
-    eta = F_to_eta2,
-    omega = F_to_omega2,
-    epsilon = F_to_epsilon2
-  )
-
-  if (verbose && !isTRUE(partial)) {
-    warning(
-      "Currently only supports partial ",
-      type,
-      " squared for repeated-measures / multi-variate ANOVAs",
-      call. = FALSE
-    )
-    partial <- TRUE
-  }
-
-  if (verbose && (isTRUE(generalized) || is.character(generalized))) {
-    warning(
-      "generalized ", type, " squared ",
-      "is not supported for this class of object."
-    )
-    generalized <- FALSE
-  }
-
   model <- stats::anova(model)
 
   p.table <- as.data.frame(model$pTerms.table)
@@ -1063,6 +1108,7 @@ cohens_f_squared <- function(model, partial = TRUE, ci = 0.95, alternative = "gr
     )
 
   attr(out, "anova_type") <- 3
+  attr(out, "approximate") <- TRUE
   out
 }
 
@@ -1131,12 +1177,12 @@ cohens_f_squared <- function(model, partial = TRUE, ci = 0.95, alternative = "gr
     aov_tab$`F` <- ifelse(aov_tab$Parameter == "Residuals", NA, 1)
     aov_tab$Mean_Square <- aov_tab$Sum_Squares/aov_tab$df
 
-    IVs <- c(id, names(attr(model, "within")), names(attr(model, "between")))
+    DV_names <- c(id, names(attr(model, "within")), names(attr(model, "between")))
 
     out <-
-      .es_aovlist(
+      .es_aov_strata(
         aov_tab,
-        IVs = IVs,
+        DV_names = DV_names,
         type = type,
         partial = partial,
         generalized = generalized,
@@ -1153,11 +1199,38 @@ cohens_f_squared <- function(model, partial = TRUE, ci = 0.95, alternative = "gr
     orig_terms <- c("(Intercept)", orig_terms)
   }
   out <- out[match(out$Parameter, orig_terms),]
+
   attr(out, "anova_type") <- attr(model, "type", exact = TRUE)
   attr(out, "approximate") <- FALSE
   out
 }
 
+#' @keywords internal
+.anova_es.mixed <- function(model,
+                            verbose = TRUE,
+                            include_intercept = FALSE,
+                            ...) {
+  aov_tab <- as.data.frame(model[["anova_table"]])
+
+  if (!"F" %in% colnames(aov_tab)) {
+    stop("Cannot estimate approx effect size for `mixed` type model - no F-statistic found.", call. = FALSE)
+  }
+
+  if (verbose && include_intercept) {
+    warning("Cannot estimate (Intercept) effect size for `mixed` model.", call. = FALSE)
+  }
+
+  aov_tab$Parameter <- rownames(aov_tab)
+  aov_tab$df <- aov_tab[["num Df"]]
+  aov_tab$df_error <- aov_tab[["den Df"]]
+  aov_tab <- aov_tab[, c("Parameter", "df", "df_error", "F")]
+
+  out <- .es_aov_table(aov_tab, verbose = verbose, ...)
+
+  attr(out, "anova_type") <- attr(model, "type")
+  attr(out, "approximate") <- TRUE
+  out
+}
 
 
 #' @keywords internal
@@ -1172,18 +1245,25 @@ cohens_f_squared <- function(model, partial = TRUE, ci = 0.95, alternative = "gr
   if (!inherits(model, "anova.rms")) {
     model <- stats::anova(model, test = "F")
   }
+
   i <- rownames(model)
   model <- as.data.frame(model)
-
-  colnames(model) <- gsub("F", "F value", colnames(model), fixed = TRUE)
-  colnames(model) <- gsub("d.f.", "NumDF", colnames(model), fixed = TRUE)
-  model$DenDF <- model$NumDF[rownames(model) == "ERROR"]
-
+  model$Parameter <- i
+  colnames(model) <- gsub("d.f.", "df", colnames(model), fixed = TRUE)
+  model$df_error <- model$df[rownames(model) == "ERROR"]
   model <- model[rownames(model) != "ERROR", ]
 
-  out <- .anova_es.anova(model, type = type, partial = partial, generalized = generalized, ci = ci, alternative = alternative, ...)
-  out$Parameter <- i[match(make.names(i), out$Parameter, nomatch = 0)]
+  out <- .es_aov_table(
+    model,
+    type = type,
+    partial = partial,
+    generalized = generalized,
+    ci = ci,
+    alternative = alternative,
+    ...
+  )
   attr(out, "anova_type") <- 2
+  attr(out, "approximate") <- FALSE
   out
 }
 
@@ -1207,14 +1287,4 @@ cohens_f_squared <- function(model, partial = TRUE, ci = 0.95, alternative = "gr
     verbose = verbose,
     ...
   )
-}
-
-# Utils -------------------------------------------------------------------
-
-#' @keywords internal
-.anova_es_model_params <- function(model, f, df_num, df_error, type, ci, alternative) {
-  # used by .anova_es.parameters_model
-  out <- .F_to_pve(stats::na.omit(f), df = df_num, df_error = df_error, ci = ci, alternative = alternative, es = paste0(type, "2"))
-  out$Parameter <- model$Parameter[!is.na(f)]
-  out[c(ncol(out), 1:(ncol(out) - 1))]
 }
\ No newline at end of file
diff --git a/R/print.effectsize_table.R b/R/print.effectsize_table.R
index 9ab03cbf5..23f1bbbb4 100644
--- a/R/print.effectsize_table.R
+++ b/R/print.effectsize_table.R
@@ -10,6 +10,8 @@
 print.effectsize_table <- function(x, digits = 2, ...) {
   x_orig <- x
 
+  footer <- attr(x, "table_footer")
+
   if (!is.null(alt <- attr(x, "alternative")) && alt != "two.sided") {
     ci_footer <- sprintf(
       "\n- One-sided CIs: %s bound fixed at (%s).",
@@ -17,10 +19,15 @@ print.effectsize_table <- function(x, digits = 2, ...) {
       as.character(if (alt == "less") x$CI_low[1] else x$CI_high[1])
     )
 
-    attr(x, "table_footer") <-
-      c(attr(x, "table_footer"), list(c(ci_footer, "cyan")))
+    footer <- c(footer, list(c(ci_footer, "cyan")))
   }
 
+  # if (isTRUE(attr(x, "approximate"))) {
+  #   footer <- c(footer, list(c("\n- Effect size is approximated.", "cyan")))
+  # }
+
+  attr(x, "table_footer") <- footer
+
   x <- format(x, digits = digits)
 
   cat(insight::export_table(x, digits = digits, ...))
@@ -192,9 +199,9 @@ print.effectsize_anova <- function(x, digits = 2, ...) {
   obs <- attr(x, "generalized")
   if (is.character(obs) || isTRUE(obs)) {
     if (isTRUE(obs)) {
-      footer <- "\n- Observed variabels: All"
+      footer <- "\n- Observed variables: All"
     } else {
-      footer <- paste0("\n- Observed variabels: ", paste0(obs, collapse = ", "))
+      footer <- paste0("\n- Observed variables: ", paste0(obs, collapse = ", "))
     }
     footer <- list(c(footer, "cyan"))
   }
diff --git a/R/standardize.models.R b/R/standardize.models.R
index 4d5a0103f..af81c96e2 100644
--- a/R/standardize.models.R
+++ b/R/standardize.models.R
@@ -14,12 +14,15 @@
 #'
 #' @param x A statistical model.
 #' @param weights If `TRUE` (default), a weighted-standardization is carried out.
-#' @param include_response For a model, if `TRUE` (default), the response value
-#'   will also be standardized. If `FALSE`, only the predictors will be
-#'   standardized. Note that for certain models (logistic regression, count
-#'   models, ...), the response value will never be standardized, to make
-#'   re-fitting the model work. (For `mediate` models, only applies to the y
-#'   model; m model's response will always be standardized.)
+#' @param include_response If `TRUE` (default), the response value will also be
+#'   standardized. If `FALSE`, only the predictors will be standardized.
+#'   - Note that for GLMs and models with non-linear link functions, the
+#'   response value will not be standardized, to make re-fitting the model work.
+#'   - If the model contains an [stats::offset()], the offset variable(s) will
+#'   be standardized only if the response is standardized. If `two_sd = TRUE`,
+#'   offsets are standardized by one-sd (similar to the response).
+#'   - (For `mediate` models, the `include_response` refers to the outcome in
+#'   the y model; m model's response will always be standardized when possible).
 #' @inheritParams datawizard::standardize
 #'
 #' @return A statistical model fitted on standardized data
@@ -50,6 +53,7 @@
 #' values, such as `log()` and `sqrt()`, the releven variables are shifted (post
 #' standardization) by `Z - min(Z) + 1` or `Z - min(Z)` (respectively).
 #'
+#'
 #' @family standardize
 #' @examples
 #' model <- lm(Infant.Mortality ~ Education * Fertility, data = swiss)
@@ -60,7 +64,8 @@
 #' @importFrom datawizard standardize
 #' @importFrom utils capture.output
 #' @export
-#' @aliases standardize-models
+#' @aliases standardize_models
+#' @aliases standardize.models
 standardize.default <- function(x,
                                 robust = FALSE,
                                 two_sd = FALSE,
@@ -68,7 +73,10 @@ standardize.default <- function(x,
                                 verbose = TRUE,
                                 include_response = TRUE,
                                 ...) {
-  m_info <- insight::model_info(x)
+  if (is.null(m_info <- list(...)[["m_info"]])){
+    m_info <- insight::model_info(x)
+  }
+
   data <- insight::get_data(x)
 
   if (insight::is_multivariate(x) && inherits(x, "brmsfit")) {
@@ -76,26 +84,50 @@ standardize.default <- function(x,
       "\nAs an alternative: you may standardize your data (and adjust your priors), and re-fit the model.",
       call. = FALSE
     )
-  } else if (m_info$is_bayesian) {
+  }
+
+  if (m_info$is_bayesian) {
     warning("Standardizing variables without adjusting priors may lead to bogus results unless priors are auto-scaled.",
       call. = FALSE, immediate. = TRUE
     )
   }
 
-  # for models with specific scale of the response value (e.g. count models
-  # with positive integers, or beta with ratio between 0 and 1), we need to
-  # make sure that the original response value will be restored after
-  # standardizing, as these models also require a non-standardized response.
-  if (!include_response || .no_response_standardize(m_info)) {
+  #  =-=-=-= Z the RESPONSE? =-=-=-=
+  # Some models have special responses that should not be standardized. This
+  # includes:
+  # - generalized linear models (counts, binomial, etc...)
+  # - Survival models
+
+  include_response <- include_response && .safe_to_standardize_response(m_info)
+
+  resp <- NULL
+  if (!include_response) {
     resp <- unique(c(insight::find_response(x), insight::find_response(x, combine = FALSE)))
   } else if (include_response && two_sd) {
     resp <- unique(c(insight::find_response(x), insight::find_response(x, combine = FALSE)))
-  } else {
-    resp <- NULL
   }
 
-  # Do not standardize weighting-variable, because negative weights will
-  # cause errors in "update()"
+  # If there's an offset, don't standardize offset OR response
+  offsets <- insight::find_offset(x)
+  if (length(offsets)) {
+    if (include_response) {
+      if (verbose) {
+        warning("Offset detected and will be standardized.", call. = FALSE)
+      }
+
+      if (two_sd) {
+        # Treat offsets like responses - only standardize by 1 SD
+        resp <- c(resp, offsets)
+        offsets <- NULL
+      }
+    } else if (!include_response) {
+      # Don't standardize offsets if not standardizing the response
+      offsets <- NULL
+    }
+  }
+
+  #  =-=-=-= DO NOT Z WEIGHTING-VARIABLE =-=-=-=
+  # because negative weights will cause errors in "update()"
   weight_variable <- insight::find_weights(x)
 
   if (!is.null(weight_variable) && !weight_variable %in% colnames(data) && "(weights)" %in% colnames(data)) {
@@ -104,10 +136,11 @@ standardize.default <- function(x,
     weight_variable <- c(weight_variable, "(weights)")
   }
 
-  # don't standardize random effects
+  #  =-=-=-= DO NOT Z RANDOM-EFFECTS =-=-=-=
   random_group_factor <- insight::find_random(x, flatten = TRUE, split_nested = TRUE)
 
-  # standardize data
+
+  #  =-=-=-= WHICH YES, WHICH NO? SUMMARY =-=-=-=
   dont_standardize <- c(resp, weight_variable, random_group_factor)
   do_standardize <- setdiff(colnames(data), dont_standardize)
 
@@ -121,76 +154,75 @@ standardize.default <- function(x,
       call. = FALSE
     )
     do_standardize <- setdiff(do_standardize, doller_vars)
+    dont_standardize <- c(dont_standardize, doller_vars)
   }
 
 
-  if (length(do_standardize)) {
-    w <- insight::get_weights(x, na_rm = TRUE)
-
-    data_std <- datawizard::standardize(data[do_standardize],
-      robust = robust,
-      two_sd = two_sd,
-      weights = if (weights) w else NULL,
-      verbose = verbose
-    )
-
-    if (!.no_response_standardize(m_info) && include_response && two_sd) {
-      # if two_sd, it must not affect the response!
-      data_std[resp] <- datawizard::standardize(data[resp],
-        robust = robust,
-        two_sd = FALSE,
-        weights = if (weights) w else NULL,
-        verbose = verbose
-      )
-      dont_standardize <- setdiff(dont_standardize, resp)
-    }
-  } else {
+  if (!length(do_standardize)) {
     warning("No variables could be standardized.", call. = FALSE)
     return(x)
   }
 
+  #  =-=-=-= STANDARDIZE! =-=-=-=
+  w <- insight::get_weights(x, na_rm = TRUE)
+
+  data_std <- datawizard::standardize(data[do_standardize],
+                                      robust = robust,
+                                      two_sd = two_sd,
+                                      weights = if (weights) w,
+                                      verbose = verbose)
+
+  # if two_sd, it must not affect the response!
+  if (include_response && two_sd) {
+    data_std[resp] <- datawizard::standardize(data[resp],
+                                              robust = robust,
+                                              two_sd = FALSE,
+                                              weights = if (weights) w,
+                                              verbose = verbose)
+
+    dont_standardize <- setdiff(dont_standardize, resp)
+  }
 
+
+  #  =-=-=-= FIX LOG-SQRT VARIABLES! =-=-=-=
   # if we standardize log-terms, standardization will fail (because log of
   # negative value is NaN). Do some back-transformation here
 
   log_terms <- .log_terms(x, data_std)
   if (length(log_terms) > 0) {
-    data_std[log_terms] <- lapply(data_std[log_terms], function(i) {
-      i - min(i, na.rm = TRUE) + 1
-    })
+    data_std[log_terms] <- lapply(data_std[log_terms],
+                                  function(i) i - min(i, na.rm = TRUE) + 1)
   }
 
   # same for sqrt
   sqrt_terms <- .sqrt_terms(x, data_std)
   if (length(sqrt_terms) > 0) {
-    data_std[sqrt_terms] <- lapply(data_std[sqrt_terms], function(i) {
-      i - min(i, na.rm = TRUE)
-    })
+    data_std[sqrt_terms] <- lapply(data_std[sqrt_terms],
+                                   function(i) i - min(i, na.rm = TRUE))
   }
 
-  if (verbose && (length(log_terms) > 0 || length(sqrt_terms) > 0)) {
+  if (verbose && length(c(log_terms, sqrt_terms))) {
     message("Formula contains log- or sqrt-terms. See help(\"standardize\") for how such terms are standardized.")
   }
 
 
-  # restore data that should not be standardized
-
+  #  =-=-=-= ADD BACK VARIABLES THAT WHERE NOT Z =-=-=-=
   if (length(dont_standardize)) {
     remaining_columns <- intersect(colnames(data), dont_standardize)
     data_std <- cbind(data[, remaining_columns, drop = FALSE], data_std)
   }
 
-  # update model with standardized data
-
+  #  =-=-=-= UPDATE MODEL WITH STANDARDIZED DATA =-=-=-=
   tryCatch({
     if (inherits(x, "brmsfit")) {
       text <- utils::capture.output(model_std <- stats::update(x, newdata = data_std))
     } else if (inherits(x, "biglm")) {
       text <- utils::capture.output(model_std <- stats::update(x, moredata = data_std))
+    } else if (inherits(x, "mixor")) {
+      data_std <- data_std[order(data_std[, random_group_factor, drop = FALSE]), ]
+      text <- utils::capture.output(model_std <- stats::update(x, data = data_std))
     } else {
-      if (inherits(x, "mixor")) {
-        data_std <- data_std[order(data_std[, random_group_factor, drop = FALSE]), ]
-      }
+      # DEFAULT METHOD
       text <- utils::capture.output(model_std <- stats::update(x, data = data_std))
     }
   }, error = function(er) {
@@ -208,74 +240,6 @@ standardize.default <- function(x,
 
 # exceptions, models that cannot use the default-method --------------------
 
-
-#' @export
-standardize.mlm <- function(x,
-                            robust = FALSE,
-                            two_sd = FALSE,
-                            weights = TRUE,
-                            verbose = TRUE,
-                            ...) {
-  standardize.default(x,
-    robust = robust, two_sd = two_sd, weights = weights, verbose = verbose,
-    include_response = FALSE, ...
-  )
-}
-
-
-#' @export
-standardize.coxph <- function(x,
-                              robust = FALSE,
-                              two_sd = FALSE,
-                              weights = TRUE,
-                              verbose = TRUE,
-                              ...) {
-
-
-  # for some models, the DV cannot be standardized when using
-  # "update()", so we only standardize model predictors
-  #
-  # survival models have some strange format for the response variable,
-  # so we don't use the default standardize function here, but
-  # use a different approach that only retrieves predictors that should
-  # be standardized.
-
-  pred <- insight::find_predictors(x, flatten = TRUE)
-  data <- insight::get_data(x)
-
-  # if we standardize log-terms, standardization will fail (because log of
-  # negative value is NaN)
-
-  log_terms <- .log_terms(x)
-  if (length(log_terms)) pred <- setdiff(pred, log_terms)
-
-  weight_variable <- insight::find_weights(x)
-  if (length(weight_variable)) pred <- setdiff(pred, weight_variable)
-
-  # standardize data, if we have anything left to standardize
-
-  if (length(pred)) {
-    w <- insight::get_weights(x, na_rm = TRUE)
-
-    data_std <- datawizard::standardize(data[pred],
-      robust = robust,
-      two_sd = two_sd,
-      weights = if (weights) w else NULL,
-      verbose = verbose
-    )
-    data[pred] <- data_std
-  }
-
-  text <- utils::capture.output(model_std <- stats::update(x, data = data))
-
-  model_std
-}
-
-
-#' @export
-standardize.coxme <- standardize.coxph
-
-
 #' @export
 #' @importFrom utils capture.output
 #' @importFrom insight get_data
@@ -360,23 +324,6 @@ standardize.mediate <- function(x,
   model_std
 }
 
-#' @keywords internal
-.rescale_fixed_values <- function(val, cov_nm,
-                                  y_data, m_data, y_data_std, m_data_std) {
-  if (cov_nm %in% colnames(y_data)) {
-    temp_data <- y_data
-    temp_data_std <- y_data_std
-  } else {
-    temp_data <- m_data
-    temp_data_std <- m_data_std
-  }
-
-  datawizard::change_scale(val,
-    to = range(temp_data_std[[cov_nm]]),
-    range = range(temp_data[[cov_nm]])
-  )
-}
-
 
 # Cannot ------------------------------------------------------------------
 
@@ -431,20 +378,47 @@ standardize.wbgee <- standardize.wbm
 
 
 #' @keywords internal
-.no_response_standardize <- function(info, verbose = TRUE) {
+.safe_to_standardize_response <- function(info, verbose = TRUE) {
   if (is.null(info)) {
     if (verbose) {
-      warning("Unable to varify if response should not be standardized.\nResponse will be standardized.",
-        immediate. = TRUE, call. = FALSE
-      )
+      warning("Unable to varify if response should not be standardized.",
+              "\nResponse will be standardized.",
+              immediate. = TRUE, call. = FALSE)
     }
-    return(FALSE)
+    return(TRUE)
   }
+
   # check if model has a response variable that should not be standardized.
-  !info$is_linear || info$is_censored || info$family == "inverse.gaussian"
+  info$is_linear &&
+    !info$family == "inverse.gaussian" &&
+    !info$is_survival &&
+    !info$is_censored
+
+  # # alternative would be to keep something like:
+  # !info$is_count &&
+  #   !info$is_ordinal &&
+  #   !info$is_multinomial &&
+  #   !info$is_beta &&
+  #   !info$is_censored &&
+  #   !info$is_binomial &&
+  #   !info$is_survival
+  # # And then treating response for "Gamma()" or "inverse.gaussian" similar to
+  # # log-terms...
+}
 
-  ## TODO alternative would be to keep the below line for checking if no std possible
-  ##      and then treat response for "Gamma()" or "inverse.gaussian" similar to log-terms
+#' @keywords internal
+.rescale_fixed_values <- function(val, cov_nm,
+                                  y_data, m_data, y_data_std, m_data_std) {
+  if (cov_nm %in% colnames(y_data)) {
+    temp_data <- y_data
+    temp_data_std <- y_data_std
+  } else {
+    temp_data <- m_data
+    temp_data_std <- m_data_std
+  }
 
-  # info$is_count | info$is_ordinal | info$is_multinomial | info$is_beta | info$is_censored | info$is_binomial | info$is_survival
-}
+  datawizard::change_scale(val,
+                           to = range(temp_data_std[[cov_nm]]),
+                           range = range(temp_data[[cov_nm]])
+  )
+}
\ No newline at end of file
diff --git a/R/standardize_info.R b/R/standardize_info.R
index acd3230dd..364878b28 100644
--- a/R/standardize_info.R
+++ b/R/standardize_info.R
@@ -25,6 +25,11 @@
 #' @importFrom parameters parameters_type
 #' @export
 standardize_info <- function(model, robust = FALSE, two_sd = FALSE, include_pseudo = FALSE, ...) {
+  UseMethod("standardize_info")
+}
+
+#' @export
+standardize_info.default <- function(model, robust = FALSE, two_sd = FALSE, include_pseudo = FALSE, ...) {
   params <- if (inherits(model, c("glmmTMB", "MixMod"))) {
     insight::find_parameters(model, effects = "fixed", component = "conditional", flatten = TRUE, ...)
   } else {
diff --git a/R/standardize_parameters.R b/R/standardize_parameters.R
index a497c6242..f12561ace 100644
--- a/R/standardize_parameters.R
+++ b/R/standardize_parameters.R
@@ -166,12 +166,13 @@ standardize_parameters.default <- function(model, method = "refit", ci = 0.95, r
   method <- match.arg(method, c("refit", "posthoc", "smart", "basic", "classic", "pseudo"))
 
   m_info <- insight::model_info(model)
-  include_response <- include_response && !.no_response_standardize(m_info, verbose = verbose)
+  include_response <- include_response && .safe_to_standardize_response(m_info, verbose = verbose)
 
   if (method == "refit") {
-    model <- standardize(model, robust = robust, two_sd = two_sd, verbose = verbose, include_response = include_response)
+    model <- standardize(model,
+                         robust = robust, two_sd = two_sd, include_response = include_response,
+                         verbose = verbose, m_info = m_info)
   }
-  mi <- insight::model_info(model)
 
   # need model_parameters to return the parameters, not the terms
   if (inherits(model, "aov")) class(model) <- class(model)[class(model) != "aov"]
@@ -183,6 +184,11 @@ standardize_parameters.default <- function(model, method = "refit", ci = 0.95, r
   coefficient_name <- attr(pars, "coefficient_name")
 
   if (method %in% c("posthoc", "smart", "basic", "classic", "pseudo")) {
+    if (m_info$is_multivariate) {
+      # TODO FIX
+      stop('Cannot post-hoc standardize multivariate models. Try using method "refit" instead.')
+    }
+
     pars <- .standardize_parameters_posthoc(pars, method, model, robust, two_sd, exponentiate, include_response, verbose)
 
     method <- attr(pars, "std_method")
@@ -194,17 +200,11 @@ standardize_parameters.default <- function(model, method = "refit", ci = 0.95, r
   colnm <- c("Component", "Response", "Group", "Parameter", head(.col_2_scale, -2), "CI", "CI_low", "CI_high")
   pars <- pars[, colnm[colnm %in% colnames(pars)]]
 
-  if (!is.null(coefficient_name) && coefficient_name == "Odds Ratio") {
-    colnames(pars)[colnames(pars) == "Coefficient"] <- "Odds_ratio"
-  }
-  if (!is.null(coefficient_name) && coefficient_name == "Risk Ratio") {
-    colnames(pars)[colnames(pars) == "Coefficient"] <- "Risk_ratio"
-  }
-  if (!is.null(coefficient_name) && coefficient_name == "IRR") {
-    colnames(pars)[colnames(pars) == "Coefficient"] <- "IRR"
+  if (!is.null(coefficient_name) && coefficient_name %in% c("Odds Ratio", "Risk Ratio", "IRR")) {
+    colnames(pars)[colnames(pars) == "Coefficient"] <- gsub(" ", "_", coefficient_name)
   }
 
-  i <- colnames(pars) %in% c("Coefficient", "Median", "Mean", "MAP", "Odds_ratio", "IRR")
+  i <- colnames(pars) %in% c("Coefficient", "Median", "Mean", "MAP", c("Odds_Ratio", "Risk_Ratio", "IRR"))
   colnames(pars)[i] <- paste0("Std_", colnames(pars)[i])
 
   ## SE attribute?
@@ -248,7 +248,7 @@ standardize_parameters.parameters_model <- function(model, method = "refit", ci
   if (is.null(model)) model <- attr(pars, "object")
 
   m_info <- insight::model_info(model)
-  include_response <- include_response && !.no_response_standardize(m_info, verbose = verbose)
+  include_response <- include_response && .safe_to_standardize_response(m_info, verbose = verbose)
 
   if (is.null(exponentiate <- attr(pars, "exponentiate"))) exponentiate <- FALSE
   pars <- .standardize_parameters_posthoc(pars, method, model, robust, two_sd, exponentiate, include_response, verbose)
@@ -295,16 +295,15 @@ standardize_parameters.bootstrap_model <-
     model <- attr(pars, "original_model")
 
     m_info <- insight::model_info(model)
-    include_response <- include_response && !.no_response_standardize(m_info, verbose = verbose)
+    include_response <- include_response && .safe_to_standardize_response(m_info, verbose = verbose)
 
     if (method == "refit") {
       stop("The 'refit' method is not supported for bootstrapped models.")
       ## But it would look something like this:
-      # model <- standardize(model, robust = robust, two_sd = two_sd, verbose = verbose)
+      # model <- standardize(model, robust = robust, two_sd = two_sd, verbose = verbose, m_info = m_info)
       # model <- parameters::bootstrap_model(model, iterations = 1000, verbose = verbose)
       # return(model)
     }
-    mi <- insight::model_info(model)
 
     # need model_parameters to return the parameters, not the terms
     if (inherits(model, "aov")) class(model) <- class(model)[class(model) != "aov"]
@@ -438,8 +437,8 @@ standardize_parameters.model_fit <-
     }
   )
 
-  if (length(i_missing) ||
-    any(to_complete <- apply(pars[, colnames(pars) %in% .col_2_scale], 1, anyNA))) {
+  to_complete <- apply(pars[, colnames(pars) %in% .col_2_scale], 1, anyNA)
+  if (length(i_missing) || any(to_complete)) {
     i_missing <- union(i_missing, which(to_complete))
 
     pars[i_missing, colnames(pars) %in% .col_2_scale] <-
@@ -467,11 +466,11 @@ standardize_posteriors <- function(model, method = "refit", robust = FALSE, two_
   object_name <- deparse(substitute(model), width.cutoff = 500)
 
   m_info <- insight::model_info(model)
-  include_response <- include_response && !.no_response_standardize(m_info, verbose = verbose)
+  include_response <- include_response && .safe_to_standardize_response(m_info, verbose = verbose)
 
   if (method == "refit") {
     model <- standardize(model, robust = robust, two_sd = two_sd, include_response = include_response,
-                         verbose = verbose)
+                         verbose = verbose, m_info = m_info)
   }
 
   pars <- insight::get_parameters(model)
diff --git a/R/utils_values_aov.R b/R/utils_values_aov.R
index e7e209255..9dd045766 100644
--- a/R/utils_values_aov.R
+++ b/R/utils_values_aov.R
@@ -1,9 +1,8 @@
 #' @keywords internal
-.values_aov <- function(params) {
+.values_aov <- function(params, group = FALSE) {
 
   # number of observations
-  # if ("Group" %in% names(params) && ("Within" %in% params$Group)) {
-  if ("Group" %in% names(params)) {
+  if (isTRUE(group)) {
     lapply(split(params, params$Group), function(.i) {
       N <- sum(.i$df) + 1
       .prepare_values_aov(.i, N)
@@ -17,16 +16,17 @@
 
 #' @keywords internal
 .prepare_values_aov <- function(params, N) {
+  iResid <- params$Parameter == "Residuals"
   # get mean squared of residuals
-  Mean_Square_residuals <- sum(params[params$Parameter == "Residuals", ]$Mean_Square)
+  Mean_Square_residuals <- sum(params[iResid, "Mean_Square"])
   # get sum of squares of residuals
-  Sum_Squares_residuals <- sum(params[params$Parameter == "Residuals", ]$Sum_Squares)
+  Sum_Squares_residuals <- sum(params[iResid, "Sum_Squares"])
   # get total sum of squares
   Sum_Squares_total <- sum(params$Sum_Squares)
   # number of terms in model
   N_terms <- nrow(params) - 1
   # df residuals
-  df_residuals <- sum(params[params$Parameter == "Residuals", ]$df)
+  df_residuals <- sum(params[iResid, "df"])
 
   list(
     "Mean_Square_residuals" = Mean_Square_residuals,
diff --git a/_pkgdown.yml b/_pkgdown.yml
index a0c98fb3a..b7ee2e142 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -73,12 +73,15 @@ reference:
   contents:
   - sd_pooled
   - is_effectsize_name
-  - effectsize_CIs
   - format_standardize
   - print.effectsize_table
 - subtitle: "Datasets"
   contents:
   - hardlyworking
+- subtitle: "Extra Docs"
+  contents:
+  - effectsize_CIs
+  - effectsize_API
 - subtitle: "Deprecated functions"
   contents:
   - effectsize_deprecated
@@ -110,6 +113,8 @@ navbar:
         href: reference/effectsize_CIs.html
       - text: "For Bayesian Models"
         href: articles/bayesian_models.html
+  - text: "Extending effectsize"
+    href: articles/effectsize_API.html
   - text: "Conversion"
     menu:
       - text: "Between d, r, OR"
diff --git a/man/effectsize_API.Rd b/man/effectsize_API.Rd
new file mode 100644
index 000000000..7d0526f17
--- /dev/null
+++ b/man/effectsize_API.Rd
@@ -0,0 +1,58 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/docs_extra.R, R/eta_squared.R
+\name{effectsize_API}
+\alias{effectsize_API}
+\alias{.es_aov_simple}
+\alias{.es_aov_strata}
+\alias{.es_aov_table}
+\title{\code{effectsize} API}
+\usage{
+.es_aov_simple(
+  aov_table,
+  type = c("eta", "omega", "epsilon"),
+  partial = TRUE,
+  generalized = FALSE,
+  ci = 0.95,
+  alternative = "greater",
+  verbose = TRUE,
+  include_intercept = FALSE
+)
+
+.es_aov_strata(
+  aov_table,
+  DV_names,
+  type = c("eta", "omega", "epsilon"),
+  partial = TRUE,
+  generalized = FALSE,
+  ci = 0.95,
+  alternative = "greater",
+  verbose = TRUE,
+  include_intercept = FALSE
+)
+
+.es_aov_table(
+  aov_table,
+  type = c("eta", "omega", "epsilon"),
+  partial = TRUE,
+  generalized = FALSE,
+  ci = 0.95,
+  alternative = "greater",
+  verbose = TRUE,
+  include_intercept = FALSE
+)
+}
+\arguments{
+\item{aov_table}{Input data frame}
+
+\item{type}{Which effect size to compute?}
+
+\item{partial, generalized, ci, alternative, verbose}{See \code{\link[=eta_squared]{eta_squared()}}.}
+
+\item{include_intercept}{Should the intercept (\code{(Intercept)}) be included?}
+
+\item{DV_names}{A character vector with the names of all the predictors,
+including the grouping variable (e.g., \code{"Subject"}).}
+}
+\description{
+Read the \href{https://easystats.github.io/effectsize/articles/effectsize_API.html}{\emph{Support functions for model extensions}} vignette.
+}
diff --git a/man/standardize.default.Rd b/man/standardize.default.Rd
index 488c6beb8..7faa82a15 100644
--- a/man/standardize.default.Rd
+++ b/man/standardize.default.Rd
@@ -2,7 +2,8 @@
 % Please edit documentation in R/standardize.models.R
 \name{standardize.default}
 \alias{standardize.default}
-\alias{standardize-models}
+\alias{standardize_models}
+\alias{standardize.models}
 \title{Re-fit a model with standardized data}
 \usage{
 \method{standardize}{default}(
@@ -33,12 +34,17 @@ outcome) (Gelman, 2008).}
 
 \item{verbose}{Toggle warnings and messages on or off.}
 
-\item{include_response}{For a model, if \code{TRUE} (default), the response value
-will also be standardized. If \code{FALSE}, only the predictors will be
-standardized. Note that for certain models (logistic regression, count
-models, ...), the response value will never be standardized, to make
-re-fitting the model work. (For \code{mediate} models, only applies to the y
-model; m model's response will always be standardized.)}
+\item{include_response}{If \code{TRUE} (default), the response value will also be
+standardized. If \code{FALSE}, only the predictors will be standardized.
+\itemize{
+\item Note that for GLMs and models with non-linear link functions, the
+response value will not be standardized, to make re-fitting the model work.
+\item If the model contains an \code{\link[stats:offset]{stats::offset()}}, the offset variable(s) will
+be standardized only if the response is standardized. If \code{two_sd = TRUE},
+offsets are standardized by one-sd (similar to the response).
+\item (For \code{mediate} models, the \code{include_response} refers to the outcome in
+the y model; m model's response will always be standardized when possible).
+}}
 
 \item{...}{Arguments passed to or from other methods.}
 }
diff --git a/tests/testthat/test-eta_squared.R b/tests/testthat/test-eta_squared.R
index eb8d54a45..e3743e772 100644
--- a/tests/testthat/test-eta_squared.R
+++ b/tests/testthat/test-eta_squared.R
@@ -2,14 +2,12 @@ if (require("testthat") && require("effectsize")) {
 
   # anova() -----------------------------------------------------------------
   test_that("anova()", {
-    m <- matrix(c(3, 1, 1),
-      nrow = 1,
-      dimnames = list(
-        "Term",
-        c("F value", "NumDF", "DenDF")
-      )
+    m <- structure(
+      c(3, 1, 1),
+      .Dim = c(1L, 3L),
+      .Dimnames = list("Term", c("F value", "NumDF", "DenDF")),
+      class = "anova"
     )
-    class(m) <- "anova"
 
     expect_error(eta_squared(m), regexp = NA)
     expect_equal(
diff --git a/tests/testthat/test-standardize_models.R b/tests/testthat/test-standardize_models.R
index 26b5e3156..61d5df75a 100644
--- a/tests/testthat/test-standardize_models.R
+++ b/tests/testthat/test-standardize_models.R
@@ -10,6 +10,14 @@ if (require("testthat") && require("effectsize")) {
     expect_equal(coef(m0), coef(model))
   })
 
+  test_that("standardize, mlm", {
+    m <- lm(cbind(mpg, hp) ~ cyl + am, data = mtcars)
+    m2 <- lm(scale(cbind(mpg, hp)) ~ scale(cyl) + scale(am), data = mtcars)
+
+    mz <- standardize(m)
+    expect_equal(coef(mz), coef(m2), ignore_attr = TRUE)
+  })
+
   test_that("standardize | errors", {
     my_lm_external_formula <- function(.dat, predicted, predictor){
       my_formula <- as.formula(paste0(predicted, "~", predictor))
@@ -177,15 +185,16 @@ if (require("testthat") && require("effectsize")) {
   test_that("variables evaluated in the environment", {
     m <- lm(mtcars$mpg ~ mtcars$cyl + am, data = mtcars)
     w <- capture_warnings(standardize(m))
-    expect_match(w[1], "mtcars$mpg", fixed = TRUE)
+    expect_true(any(grepl("mtcars$mpg", w, fixed = TRUE)))
+
 
     skip_if(packageVersion("base") == package_version(3.4))
     ## Note:
     # No idea why this is suddenly not giving a warning on older R versions.
     m <- lm(mtcars$mpg ~ mtcars$cyl + mtcars$am, data = mtcars)
     warns <- capture_warnings(standardize(m))
-    expect_true(grepl("mtcars$mpg", warns[1], fixed = TRUE))
-    expect_true(grepl("No variables", warns[2], fixed = TRUE))
+    expect_true(any(grepl("mtcars$mpg", warns, fixed = TRUE)))
+    expect_true(any(grepl("No variables", warns, fixed = TRUE)))
   })
 
 
@@ -221,4 +230,22 @@ if (require("testthat") && require("effectsize")) {
       tolerance = 0.1
     )
   })
+
+  # Offsets -----------------------------------------------------------------
+
+  test_that("offsets", {
+    m <- lm(mpg ~ hp + offset(wt), data = mtcars)
+
+    expect_warning(mz1 <- standardize(m))
+    expect_warning(mz2 <- standardize(m, two_sd = TRUE))
+    expect_equal(c(1, 2) * coef(mz1), coef(mz2))
+
+
+    m <- glm(cyl ~ hp + offset(wt), family = poisson(), data = mtcars)
+    expect_warning(mz <- standardize(m), regexp = NA)
+
+    par1 <- parameters::model_parameters(mz)
+    par2 <- standardize_parameters(m, method = "basic")
+    expect_equal(par2[2,2], par1[2,2], tolerance = 0.05)
+  })
 }
diff --git a/tests/testthat/test-standardize_parameters.R b/tests/testthat/test-standardize_parameters.R
index f0f5c8527..34ce7fc60 100644
--- a/tests/testthat/test-standardize_parameters.R
+++ b/tests/testthat/test-standardize_parameters.R
@@ -306,7 +306,7 @@ if (require("testthat") && require("effectsize")) {
       ID = sort(rep(letters, length.out = 1000))
     )
     dat <- transform(dat, Y = X + Z + rnorm(1000))
-    dat <- cbind(dat, parameters::demean(dat, c("X", "Z"), "ID"))
+    dat <- cbind(dat, datawizard::demean(dat, c("X", "Z"), "ID"))
 
 
     m <- lme4::lmer(Y ~ scale(X_within) * X_between + C + (scale(X_within) | ID),
diff --git a/vignettes/effectsize_API.Rmd b/vignettes/effectsize_API.Rmd
new file mode 100644
index 000000000..a74c726c0
--- /dev/null
+++ b/vignettes/effectsize_API.Rmd
@@ -0,0 +1,271 @@
+---
+title: "Support functions for model extensions"
+output: 
+  rmarkdown::html_vignette:
+    toc: true
+    fig_width: 10.08
+    fig_height: 6
+tags: [r, effect size, ANOVA, standardization, standardized coefficients]
+vignette: >
+  \usepackage[utf8]{inputenc}
+  %\VignetteIndexEntry{Support functions for model extensions}
+  %\VignetteEngine{knitr::rmarkdown}
+editor_options: 
+  chunk_output_type: console
+bibliography: bibliography.bib
+---
+
+```{r message=FALSE, warning=FALSE, include=FALSE}
+library(knitr)
+knitr::opts_chunk$set(comment = ">",
+                      warning = FALSE,
+                      message = FALSE)
+options(digits = 2)
+options(knitr.kable.NA = '')
+
+
+pkgs <- c("effectsize")
+if (!all(sapply(pkgs, requireNamespace, quietly = TRUE))) {
+  knitr::opts_chunk$set(eval = FALSE)
+}
+
+set.seed(333)
+```
+
+```{r}
+library(effectsize)
+```
+
+# Supporting ANOVA Effect Sizes
+
+To add support for you model, create a new `.anova_es()` method function. This functions should generally do 3 things:
+
+1. Build a data frame with all the required information.
+2. Pass the data frame to one of the 3 functions.
+3. Set some attributes to the output.
+
+## Simple ANOVA tables
+
+The input data frame must have these columns:
+- `Parameter` (char) - The name of the parameter or, more often, the term.
+- `Sum_Squares` (num) - The sum of squares.
+- `df` (num) - The degrees of freedom associated with the `Sum_Squares`.
+- `Mean_Square_residuals` (num; *optional*) - if *not* present, is calculated as `Sum_Squares / df`.
+(Any other column is ignored.)
+
+And exactly *1* row Where `Parameter` is `Residual`.
+
+Optionally, one of the rows can have a `(Intercept)` value for `Parameter`.
+
+An example of a minimally valid data frame:
+
+```{r}
+min_aov <- data.frame(
+  Parameter = c("(Intercept)", "A", "B", "Residuals"),
+  Sum_Squares = c(30, 40, 10, 100),
+  df = c(1, 1, 2, 50)
+)
+```
+
+Pass the data frame to `.es_aov_simple()`:
+
+```{r}
+.es_aov_simple(
+  min_aov,
+  type = "eta", partial = TRUE, generalized = FALSE,
+  include_intercept = FALSE,
+  ci = 0.95, alternative = "greater",
+  verbose = TRUE
+)
+```
+
+The output is a data frame with the columns: `Parameter`, the effect size, and (optionally) `CI` + `CI_low` + `CI_high`,
+
+And with the following attributes: `partial`, `generalized`, `ci`, `alternative`, `anova_type` (`NA` or `NULL`), `approximate`.
+
+You can then set the `anova_type` attribute to {1, 2, 3, or `NA`} and return the output.
+
+## ANOVA Tables with Multiple Error Strata
+
+(e.g., `aovlist` models.)
+
+The input data frame must have these columns:
+
+- `Group` (char) - The strata
+- `Parameter` (char)
+- `Sum_Squares` (num)
+- `df` (num)
+- `Mean_Square_residuals` (num; *optional*)
+
+And exactly *1* row ***per `Group`*** Where `Parameter` is `Residual`.
+
+Optionally, one of the rows can have a `(Intercept)` value for `Parameter`.
+
+An example of a minimally valid data frame:
+
+```{r}
+min_aovlist <- data.frame(
+  Group = c("S", "S", "S:A", "S:A"),
+  Parameter = c("(Intercept)", "Residuals", "A", "Residuals"),
+  Sum_Squares = c(34, 21, 34, 400),
+  df = c(1, 12, 4, 30)
+)
+```
+
+Pass the data frame to `.es_aov_strata()`, along with a list of predictors (including the stratifying variables) to the `DV_names` argument:
+
+```{r}
+.es_aov_strata(
+  min_aovlist, DV_names = c("S", "A"),
+  type = "omega", partial = TRUE, generalized = FALSE,
+  ci = 0.95, alternative = "greater",
+  verbose = TRUE,
+  include_intercept = TRUE
+)
+```
+
+The output is a data frame with the columns: `Group`, `Parameter`, the effect size, and (optionally) `CI` + `CI_low` + `CI_high`,
+
+And with the following attributes: `partial`, `generalized`, `ci`, `alternative`, `approximate`.
+
+You can then set the `anova_type` attribute to {1, 2, 3, or `NA`} and return the output.
+
+
+## Approximate Effect sizes
+
+When *sums of squares* cannot be extracted, we can still get *approximate* effect sizes based on the `F_to_eta2()` family of functions.
+
+The input data frame must have these columns:
+
+- `Parameter` (char)
+- `F` (num) - The *F* test statistic.
+- `df` (num) - effect degrees of freedom.
+- (Can also have a `t` col instead, in which case `df` is set to 1, and `F` is `t^2`).
+- `df_error` (num) - error degrees of freedom.
+
+Optionally, one of the rows can have `(Intercept)` as the `Parameter`.
+
+An example of a minimally valid data frame:
+
+```{r}
+min_anova <- data.frame(
+  Parameter = c("(Intercept)", "A", "B"),
+  F = c(4, 7, 0.7),
+  df = c(1, 1, 2),
+  df_error = 34
+)
+```
+
+Pass the table to `.es_aov_table()`:
+
+```{r}
+.es_aov_table(
+  min_anova,
+  type = "eta", partial = TRUE, generalized = FALSE,
+  include_intercept = FALSE,
+  ci = 0.95, alternative = "greater",
+  verbose = TRUE
+)
+```
+
+The output is a data frame with the columns: `Parameter`, the effect size, and (optionally) `CI` + `CI_low` + `CI_high`,
+
+And with the following attributes: `partial`, `generalized`, `ci`, `alternative`, `approximate`.
+
+You can then set the `anova_type` attribute to {1, 2, 3, or `NA`} and return the output, and optionally the `approximate` attribute, and return the output.
+
+## *Example*
+
+Let's fit a simple linear model and change its class:
+
+```{r}
+mod <- lm(mpg ~ factor(cyl) + am, mtcars)
+
+class(mod) <- "superMODEL"
+```
+
+We now need a new `.anova_es.superMODEL` function:
+
+```{r}
+.anova_es.superMODEL <- function(model, ...) {
+  # Get ANOVA table
+  anov <- suppressWarnings(stats:::anova.lm(model))
+  anov <- as.data.frame(anov)
+  
+  # Clean up
+  anov[["Parameter"]] <- rownames(anov)
+  colnames(anov)[2:1] <- c("Sum_Squares", "df")
+  
+  # Pass
+  out <- .es_aov_simple(anov, ...)
+  
+  # Set attribute
+  attr(out, "anova_type") <- 1
+  
+  out
+}
+```
+
+And... that's it! Our new `superMODEL` class of models is fully supported!
+
+```{r}
+eta_squared(mod)
+
+eta_squared(mod, partial = FALSE)
+
+omega_squared(mod)
+
+# Etc...
+```
+
+
+# Supporting Model Re-Fitting with Standardized Data
+
+`effectsize::standardize.default()` should support your model if you have methods for:
+
+1. `{insight}` functions.
+2. An `update()` method that can take the model and a data frame via the `data = ` argument.
+
+Or you can make your own `standardize.my_class()` function, DIY-style (possibly using `datawizard::standardize.data.frame()` or `datawizard::standardize.numeric()`). This function should return a fiffed model of the same class as the input model.
+
+# Supporting Standardized Parameters
+
+`standardize_parameters.default()` offers a few methods of parameter standardization:
+
+- For `method = "refit"` all you need is to have `effectsize::standardize()` support (see above) as well as `parameters::model_parameters()`.  
+- ***API for post-hoc methods coming soon...***  
+
+<!-- `standardize_parameters.default()` should support your model if it is already supported by `{parameters}` and `{insight}`. -->
+
+<!-- - For `method = "refit"`, to have `effectsize::standardize()` support (see above). -->
+<!-- - For the post-hoc methods, you will need to have a method for `standardize_info()` (or use the default method). See next section. -->
+
+<!-- Or you can make your own `standardize_parameters.my_class()` and/or `standardize_info.my_class()` functions. -->
+
+<!-- ## Extracting Post-Hoc Standardization Information (`standardize_info`) -->
+
+<!-- The `standardize_info()` function computes the standardized units needed for standardization; In order to standardize some slope $b_{xi}$, we need to multiply it by a scaling factor: -->
+
+<!-- $$b^*_{xi} = \frac{\text{Deviation}_{xi}}{\text{Deviation}_{y}}\times b_{xi}$$ -->
+
+<!-- These "deviations" are univariate scaling factors of the response and the specific parameter (usin its corresponding feature in the design matrix). Most often these are a single standard deviation (*SD*), but depending on the `robust` and `two_sd` arguments, these can be also be two *MAD*s, etc. -->
+
+<!-- Let's look at an example: -->
+
+<!-- ```{r} -->
+<!-- m <- lm(mpg ~ factor(cyl) * am, data = mtcars) -->
+
+<!-- standardize_info(m) -->
+<!-- ``` -->
+
+<!-- - The first 4 columns (`Parameter`, `Type`, `Link`, `Secondary_Parameter`) are taken from `parameters::parameters_type()`.   -->
+<!-- - The `EffectSize_Type` column is not used here, but is used in the the `{report}` package.   -->
+<!-- - `Deviation_Response_Basic` and `Deviation_Response_Smart` correspond to the $\text{Deviation}_{y}$ scalar using two different methods of post-hoc standardization (see `standardize_parameters()` docs for more details).   -->
+<!--     - Note then when the response is not standardized (either due to `standardize_parameters(include_response = FALSE)` or because the model uses a non-continuous response), both methods are fixed at **1** (i.e., no standardization with respect to the outcome).   -->
+<!-- - `Deviation_Basic` and `Deviation_Smart` correspond to the $\text{Deviation}_{xi}$ scaler using two different methods of post-hoc standardization. -->
+
+<!-- This information is then used by the `standardize_parameters()` to standardize the parameters. -->
+    
+
+# References
+
diff --git a/vignettes/standardize_parameters.Rmd b/vignettes/standardize_parameters.Rmd
index 9554fc7d4..029151eb3 100644
--- a/vignettes/standardize_parameters.Rmd
+++ b/vignettes/standardize_parameters.Rmd
@@ -338,7 +338,7 @@ are also called *pseudo*-standardized coefficients.[^Note that like method
 `"basic"`, these are based on the model matrix.]
 
 ```{r, eval=knitr::opts_chunk$get("eval") && require(lme4) && require(lmerTest), warning=FALSE}
-m <- lme4::lmer(mpg ~ cyl + am + vs + (1|cyl), mtcars)
+m <- lme4::lmer(Reaction ~ Days + (Days|Subject), data = lme4::sleepstudy)
 
 standardize_parameters(m, method = "pseudo", ci_method = "satterthwaite")