diff --git a/R/add_discriminator_auc.R b/R/add_discriminator_auc.R
index cc8e213..8175333 100644
--- a/R/add_discriminator_auc.R
+++ b/R/add_discriminator_auc.R
@@ -4,18 +4,19 @@
 #' added using add_propensities())
 #' @param split A logical for if the metric should be calculated separately for 
 #' the training/testing split. Defaults to TRUE.
-#'
+#' @param group A list of variable names to group by
+#' 
 #' @return A discrimination object with propensities (likely added using 
 #' add_propensities()) with discriminator AUC
 #' 
 #' @export
 #'
-add_discriminator_auc <- function(discrimination, split = TRUE) {
+add_discriminator_auc <- function(discrimination, group = c(),split = TRUE) {
   
   if (split) {
     
     discriminator_auc <- discrimination$propensities %>%
-      dplyr::group_by(.data$.sample) %>%
+      dplyr::group_by(across(all_of(c(".sample", group))))%>%
       yardstick::roc_auc(".source_label", ".pred_synthetic") %>%
       dplyr::mutate(.sample = factor(.data$.sample, levels = c("training", "testing"))) %>%
       dplyr::arrange(.data$.sample)
@@ -23,6 +24,7 @@ add_discriminator_auc <- function(discrimination, split = TRUE) {
   } else {
     
     discriminator_auc <- discrimination$propensities %>%
+      dplyr::group_by(across(all_of(group)))%>%
       yardstick::roc_auc(".source_label", ".pred_synthetic") %>%
       dplyr::mutate(.sample = factor("overall", levels = "overall"))
     
diff --git a/R/add_pmse.R b/R/add_pmse.R
index 02f2f3d..144f4f0 100644
--- a/R/add_pmse.R
+++ b/R/add_pmse.R
@@ -1,9 +1,11 @@
-#' Add pMSE to discrimination object
+#' Add pMSE to discrimination object, with option to assess separately on groups
+#' indicated by a grouping variable
 #'
 #' @param discrimination A discrimination object with propensities (likely 
 #' added using add_propensities())
 #' @param split A logical for if the metric should be calculated separately for 
 #' the training/testing split. Defaults to TRUE.
+#' @param group A set of variables to group the pmse by
 #'
 #' @return A discrimination object with propensities (likely added using 
 #' add_propensities()) with a pMSE
@@ -12,27 +14,30 @@
 #' 
 #' @export
 #' 
-add_pmse <- function(discrimination, split = TRUE) {
+
+add_pmse <- function(discrimination, group = c(), split = TRUE) {
   
   calc_pmse <- function(propensities) {
     
     # calculate the expected propensity
     prop_synthetic <- propensities %>%
-      dplyr::summarize(
-        n_synthetic = sum(.data$.source_label == "synthetic"),
-        n_total = dplyr::n()
-      ) %>%
-      dplyr::mutate(prop_synthetic = .data$n_synthetic / .data$n_total) %>%
-      dplyr::pull("prop_synthetic")
+      dplyr::group_by(across(all_of(group))) %>%
+      dplyr::summarize("prop_synthetic" = list(c(sum(.data$.source_label == "synthetic")/dplyr::n()))) %>%
+      dplyr::pull(prop_synthetic)
     
     propensities_vec <- propensities %>%
-      dplyr::pull(".pred_synthetic")
+      dplyr::group_by(across(all_of(group))) %>%
+      dplyr::summarise(".pred_synthetic" = list(c(.pred_synthetic))) %>%
+      dplyr::pull(.pred_synthetic)
     
+    # function for pmse
+    pmse_func <- function(propensities_vec, prop_synthetic){
+      mean((propensities_vec - prop_synthetic)^2)
+    }
     # calculate the observed pMSE
-    pmse <- mean((propensities_vec - prop_synthetic) ^ 2)
+    pmse <- mapply(pmse_func, propensities_vec, prop_synthetic)
     
     return(pmse)
-    
   }
   
   if (split) {
@@ -45,21 +50,47 @@ add_pmse <- function(discrimination, split = TRUE) {
       dplyr::filter(.data$.sample == "testing") %>%
       calc_pmse()
     
-    pmse <- tibble::tibble(
-      .source = factor(c("training", "testing"), levels = c("training", "testing")),
-      .pmse = c(pmse_training, pmse_testing)
-    )
+    if (length(group)==0){ # original case
+      pmse <- tibble::tibble(
+        .source = factor(c("training", "testing"), levels = c("training", "testing")),
+        .pmse = c(pmse_training, pmse_testing)
+      )
+    }
+    else{ # we have passed a list of grouping variables
+      groups <- discrimination$propensities %>%
+        dplyr::group_by(across(all_of(group))) %>%
+        group_keys()
+      groups <- rbind(groups, groups) # make 2 copies for train/test
+      
+      pmse <- tibble::tibble(
+        groups,
+        .source = factor(c(rep("training", length(pmse_training)), rep("testing", length(pmse_testing))), levels = c("training", "testing")),
+        .pmse = c(pmse_training, pmse_testing)
+      )
+    }
     
   } else {
     
     pmse_overall <- discrimination$propensities %>%
       calc_pmse()
     
-    pmse <- tibble::tibble(
-      .source = factor("overall", levels = "overall"),
-      .pmse = pmse_overall
-    )
-    
+    if (length(group)==0){ # original case
+      pmse <- tibble::tibble(
+        .source = factor("overall", levels = "overall"),
+        .pmse = pmse_overall
+      )
+    }
+    else{ # we have passed a list of grouping variables
+      groups <- discrimination$propensities %>%
+        dplyr::group_by(across(all_of(group))) %>%
+        group_keys()
+      
+      pmse <- tibble::tibble(
+        groups,
+        .source = factor(c(rep("overall", length(pmse_overall))), levels = c("overall")),
+        .pmse = pmse_overall
+      )
+    }
   }
   
   discrimination$pmse <- pmse
diff --git a/R/add_pmse_ratio.R b/R/add_pmse_ratio.R
index 393c8bf..622b11c 100644
--- a/R/add_pmse_ratio.R
+++ b/R/add_pmse_ratio.R
@@ -6,6 +6,7 @@
 #' @param prop The proportion of data to be retained for modeling/analysis in 
 #' the training/testing split. The sampling is stratified by the original and
 #' synthetic data.
+#' @param group A set of variables to group the pmse and null pmse by
 #' @param times The number of bootstrap samples.
 #'
 #' @return A discrimination with pMSE
@@ -13,7 +14,7 @@
 #' @family Utility metrics
 #' 
 #' @export
-add_pmse_ratio <- function(discrimination, split = TRUE, prop = 3 / 4, times) {
+add_pmse_ratio <- function(discrimination, split = TRUE, prop = 3 / 4, group = c(), times) {
   
   if (is.null(discrimination$pmse)) {
     
@@ -25,48 +26,86 @@ add_pmse_ratio <- function(discrimination, split = TRUE, prop = 3 / 4, times) {
     
     # calculate the expected propensity
     prop_synthetic <- propensities %>%
-      dplyr::summarize(
-        n_synthetic = sum(.data$.source_label == "synthetic"),
-        n_total = dplyr::n()
-      ) %>%
-      dplyr::mutate(prop_synthetic = .data$n_synthetic / .data$n_total) %>%
-      dplyr::pull("prop_synthetic")
+      dplyr::group_by(across(all_of(group))) %>%
+      dplyr::summarize("prop_synthetic" = list(c(sum(.data$.source_label == "synthetic")/dplyr::n()))) %>%
+      dplyr::pull(prop_synthetic)
     
     propensities_vec <- propensities %>%
-      dplyr::pull(".pred_synthetic")
+      dplyr::group_by(across(all_of(group))) %>%
+      dplyr::summarise(".pred_synthetic" = list(c(.pred_synthetic))) %>%
+      dplyr::pull(.pred_synthetic)
     
+    # function for pmse
+    pmse_func <- function(propensities_vec, prop_synthetic){
+      mean((propensities_vec - prop_synthetic)^2)
+    }
     # calculate the observed pMSE
-    pmse <- mean((propensities_vec - prop_synthetic) ^ 2)
+    pmse <- mapply(pmse_func, propensities_vec, prop_synthetic)
     
     return(pmse)
-    
   }
-
-  pmse_null_overall <- vector(mode = "numeric", length = times)
-  pmse_null_training <- vector(mode = "numeric", length = times)
-  pmse_null_testing <- vector(mode = "numeric", length = times)
+  # function to sample 2x size of grouped data, by group, 
+  # for the dataset with both confidential and synthetic data
+  # group_resample <- function(data){
+  #   # split by grouping variables
+  #   split_data = dplyr::group_split(data %>% dplyr::ungroup() %>% dplyr::group_by(across(all_of(group))))
+  #   bootstrap_sample = list()
+  #   for (elem in split_data){ # iterate through grouped dataset
+  #     # in each group, sample twice as much data
+  #     bootstrap_sample <- append(bootstrap_sample, list(dplyr::bind_cols(
+  #       elem %>%
+  #         dplyr::filter(.data$.source_label == "original") %>%
+  #         dplyr::group_by(across(all_of(group))) %>%
+  #         dplyr::slice_sample(n = nrow(elem), replace = TRUE) %>%
+  #         dplyr::select(-".source_label"),
+  #       elem %>%
+  #         dplyr::select(".source_label"))
+  #     ))
+  #   }
+  #   bootstrap_sample = dplyr::bind_rows(bootstrap_sample)
+  # }
   
-  for (a in seq_along(pmse_null_overall)) {
-    
+  group_resample <- function(data){
+    # split by grouping variables
+    split_data = dplyr::group_split(data %>% dplyr::ungroup() %>% dplyr::group_by(across(all_of(group))))
+    bootstrap_sample = vector("list", length = length(split_data))
+    for (i in 1:length(split_data)){ # iterate through grouped dataset
+      # in each group, sample twice as much data
+      bootstrap_sample[[i]] <- list(dplyr::bind_cols(
+        split_data[[i]] %>%
+          dplyr::filter(.data$.source_label == "original") %>%
+          dplyr::group_by(across(all_of(group))) %>%
+          dplyr::slice_sample(n = nrow(split_data[[i]]), replace = TRUE) %>%
+          dplyr::select(-".source_label"),
+        split_data[[i]] %>%
+          dplyr::select(".source_label"))
+      )
+    }
+    bootstrap_sample = dplyr::bind_rows(bootstrap_sample)
+  }
+  
+  # matrix instead of vector, where each entry is a simulation, containing a vector with groups
+  #pmse_null_overall <- c()
+  #pmse_null_training <- c()
+  #pmse_null_testing <- c()
+  
+  pmse_null_overall <- rep(NA, times)
+  pmse_null_training <- rep(NA, times)
+  pmse_null_testing <- rep(NA, times)
+  
+  for (a in 1:times) {
     # bootstrap sample original observations to equal the size of the combined 
-    # data
+    # data, with a vector of one set of observations per grouping variable (?)
     # append the original labels so the proportions match
-    bootstrap_sample <- dplyr::bind_cols(
-      discrimination$combined_data %>%
-        dplyr::filter(.data$.source_label == "original") %>%
-        dplyr::slice_sample(n = nrow(discrimination$combined_data), replace = TRUE) %>%
-        dplyr::select(-".source_label"),
-      discrimination$combined_data %>%
-        dplyr::select(".source_label")
-    )
-
+    bootstrap_sample <- group_resample(discrimination$combined_data)
+    
     if (split) {
       
-      # make training/testing split
+      # make training/testing split 
       data_split <- rsample::initial_split(
         data = bootstrap_sample,
         prop = prop,
-        strata = ".source_label"
+        strata = ".source_label" # NOTE: Should this also be stratified by group?
       )
       
       # fit the model from the pMSE on the bootstrap sample
@@ -89,14 +128,20 @@ add_pmse_ratio <- function(discrimination, split = TRUE, prop = 3 / 4, times) {
         )
       
       # calculate the pmse for each bootstrap
+      # pmse_null_overall <- append(pmse_null_overall, calc_pmse(propensities_df))
+      # pmse_null_training <- append(pmse_null_training, propensities_df %>%
+      #                                dplyr::filter(.data$.sample == "training") %>%
+      #                                calc_pmse())
+      # pmse_null_testing <- append(pmse_null_testing, propensities_df %>%
+      #                               dplyr::filter(.data$.sample == "testing") %>%
+      #                               calc_pmse())
       pmse_null_overall[a] <- calc_pmse(propensities_df)
       pmse_null_training[a] <- propensities_df %>%
-        dplyr::filter(.data$.sample == "training") %>%
-        calc_pmse()
+                                     dplyr::filter(.data$.sample == "training") %>%
+                                     calc_pmse()
       pmse_null_testing[a] <- propensities_df %>%
-        dplyr::filter(.data$.sample == "testing") %>%
-        calc_pmse()
-      
+                                    dplyr::filter(.data$.sample == "testing") %>%
+                                    calc_pmse()
     } else {
       
       # fit the model from the pMSE on the bootstrap sample
@@ -113,15 +158,17 @@ add_pmse_ratio <- function(discrimination, split = TRUE, prop = 3 / 4, times) {
       
       # calculate the pmse for each bootstrap
       pmse_null_overall[a] <- calc_pmse(propensities_df)
-    
+      
     }
     
   }
   
   # find the mean of the bootstrapped pMSEs
-  mean_null_pmse_overall <- mean(pmse_null_overall)
-  mean_null_pmse_training <- mean(pmse_null_training)
-  mean_null_pmse_testing <- mean(pmse_null_testing)
+  mean_null_pmse_overall <- colMeans(t(matrix(pmse_null_overall, ncol = times))) # each row is a new sample
+  if (split){
+    mean_null_pmse_training <- colMeans(t(matrix(pmse_null_training, ncol = times)))
+    mean_null_pmse_testing <- colMeans(t(matrix(pmse_null_testing, ncol= times)))
+  }
   
   # calculate the ratio for the training/testing split or overall data
   if (all(c("training", "testing") %in% discrimination$pmse$.source)) {
diff --git a/R/add_specks.R b/R/add_specks.R
index cac92fd..04f5be9 100644
--- a/R/add_specks.R
+++ b/R/add_specks.R
@@ -3,31 +3,47 @@
 #' @param discrimination A discrimination with added propensities
 #' @param split A logical for if the metric should be calculated separately for 
 #' the training/testing split. Defaults to TRUE.
+#' @param group A list of variable names to group by
 #'
 #' @family Utility metrics
 #'
 #' @return A discrimination with SPECKS
 #' 
 #' @export
-add_specks <- function(discrimination, split = TRUE) {
+
+add_specks <- function(discrimination, group = c(), split = TRUE) {
   
   calc_specks <- function(propensities) {
     
     propensities_original <- propensities %>%
+      dplyr::group_by(across(all_of(group))) %>%
       dplyr::filter(.data$.source_label == "original") %>%
+      dplyr::summarize(".pred_synthetic" = list(c(.data$.pred_synthetic))) %>%
       dplyr::pull(".pred_synthetic")
     
     propensities_synthetic <- propensities %>%
+      dplyr::group_by(across(all_of(group))) %>%
       dplyr::filter(.data$.source_label == "synthetic") %>%
+      dplyr::summarize(".pred_synthetic" = list(c(.data$.pred_synthetic))) %>%
       dplyr::pull(".pred_synthetic")
     
+    specks_func <- function(original, synthetic){
+      suppressWarnings(
+        stats::ks.test(
+          original, 
+          synthetic, 
+          exact = FALSE)$statistic
+      )
+    }
     # Calculate KS Distance of the original and synthetic ECDFS
-    specks <- suppressWarnings(
-      stats::ks.test(
-        propensities_original, 
-        propensities_synthetic, 
-        exact = FALSE)$statistic
-    )
+    #specks <- suppressWarnings(
+    # stats::ks.test(
+    #    propensities_original, 
+    #    propensities_synthetic, 
+    #    exact = FALSE)$statistic
+    #)
+    # calculate the observed pMSE
+    specks <- mapply(specks_func, propensities_original, propensities_synthetic)
     
     specks <- unname(specks)
     
@@ -45,25 +61,52 @@ add_specks <- function(discrimination, split = TRUE) {
       dplyr::filter(.data$.sample == "testing") %>%
       calc_specks()
     
-    specks <- tibble::tibble(
-      .source = factor(c("training", "testing"), levels = c("training", "testing")),
-      .specks = c(specks_training, specks_testing)
-    )
+    if (length(group)==0){ # original case
+      specks <- tibble::tibble(
+        .source = factor(c("training", "testing"), levels = c("training", "testing")),
+        .specks = c(specks_training, specks_testing)
+      )
+    }
+    else{ # we have passed a list of grouping variables
+      groups <- discrimination$propensities %>%
+        dplyr::group_by(across(all_of(group))) %>%
+        group_keys()
+      groups <- rbind(groups, groups) # make 2 copies for train/test
+      
+      specks <- tibble::tibble(
+        groups,
+        .source = factor(c(rep("training", length(specks_training)), rep("testing", length(specks_testing))), levels = c("training", "testing")),
+        .specks = c(specks_training, specks_testing)
+      )
+    }
     
   } else {
     
     specks_overall <- discrimination$propensities %>%
       calc_specks()
     
-    specks <- tibble::tibble(
-      .source = factor("overall", levels = "overall"),
-      .specks = specks_overall
-    )
+    if (length(group)==0){ # original case
+      specks <- tibble::tibble(
+        .source = factor("overall", levels = "overall"),
+        .specks = specks_overall
+      )
+    }
+    else{ # we have passed a list of grouping variables
+      groups <- discrimination$propensities %>%
+        dplyr::group_by(across(all_of(group))) %>%
+        group_keys()
+      
+      specks <- tibble::tibble(
+        groups,
+        .source = factor(c(rep("overall", length(specks_overall))), levels = c("overall")),
+        .specks = specks_overall
+      )
+    }
     
   }
   
   discrimination$specks <- specks
   
   return(discrimination)
-
+  
 }
\ No newline at end of file
diff --git a/README.qmd b/README.qmd
index ad1d683..ea83a4b 100644
--- a/README.qmd
+++ b/README.qmd
@@ -267,7 +267,6 @@ disc1 %>%
   add_specks() %>%
   add_pmse() %>%
   add_pmse_ratio(times = 25)
-
 ```
 
 Finally, we can look at variable importance and the decision tree from our discriminator. 
@@ -355,6 +354,15 @@ util_moments(
 
 ```
 
+All four of discriminant-based metrics can also be evaluated by group
+```{r}
+disc1 %>%
+  add_discriminator_auc(group = c("sex", "species")) %>%
+  add_specks(group = c("sex", "species")) %>%
+  add_pmse(group = c("sex", "species")) %>%
+  add_pmse_ratio(times = 25, group = c("sex", "species"))
+```
+
 ### Weighting
 
 Many utility metrics include a `weight_var` argument to use weighted statistics during calculation. For example, this code weights the moments by the body weight of the penguins.
diff --git a/man/add_discriminator_auc.Rd b/man/add_discriminator_auc.Rd
index 84e2cbf..e8914e8 100644
--- a/man/add_discriminator_auc.Rd
+++ b/man/add_discriminator_auc.Rd
@@ -4,12 +4,14 @@
 \alias{add_discriminator_auc}
 \title{Add discriminator AUC to discrimination object}
 \usage{
-add_discriminator_auc(discrimination, split = TRUE)
+add_discriminator_auc(discrimination, group = c(), split = TRUE)
 }
 \arguments{
 \item{discrimination}{A discrimination object with propensities (likely
 added using add_propensities())}
 
+\item{group}{A list of variable names to group by}
+
 \item{split}{A logical for if the metric should be calculated separately for
 the training/testing split. Defaults to TRUE.}
 }
diff --git a/man/add_pmse.Rd b/man/add_pmse.Rd
index 2a12257..7b44943 100644
--- a/man/add_pmse.Rd
+++ b/man/add_pmse.Rd
@@ -2,14 +2,17 @@
 % Please edit documentation in R/add_pmse.R
 \name{add_pmse}
 \alias{add_pmse}
-\title{Add pMSE to discrimination object}
+\title{Add pMSE to discrimination object, with option to assess separately on groups
+indicated by a grouping variable}
 \usage{
-add_pmse(discrimination, split = TRUE)
+add_pmse(discrimination, group = c(), split = TRUE)
 }
 \arguments{
 \item{discrimination}{A discrimination object with propensities (likely
 added using add_propensities())}
 
+\item{group}{A set of variables to group the pmse by}
+
 \item{split}{A logical for if the metric should be calculated separately for
 the training/testing split. Defaults to TRUE.}
 }
@@ -18,7 +21,8 @@ A discrimination object with propensities (likely added using
 add_propensities()) with a pMSE
 }
 \description{
-Add pMSE to discrimination object
+Add pMSE to discrimination object, with option to assess separately on groups
+indicated by a grouping variable
 }
 \seealso{
 Other Utility metrics: 
diff --git a/man/add_pmse_ratio.Rd b/man/add_pmse_ratio.Rd
index fc404bb..9b9f49f 100644
--- a/man/add_pmse_ratio.Rd
+++ b/man/add_pmse_ratio.Rd
@@ -4,7 +4,7 @@
 \alias{add_pmse_ratio}
 \title{Add pMSE ratio to discrimination object}
 \usage{
-add_pmse_ratio(discrimination, split = TRUE, prop = 3/4, times)
+add_pmse_ratio(discrimination, split = TRUE, prop = 3/4, group = c(), times)
 }
 \arguments{
 \item{discrimination}{A discrimination with added propensities}
@@ -16,6 +16,8 @@ the training/testing split. Defaults to TRUE.}
 the training/testing split. The sampling is stratified by the original and
 synthetic data.}
 
+\item{group}{A set of variables to group the pmse and null pmse by}
+
 \item{times}{The number of bootstrap samples.}
 }
 \value{
diff --git a/man/add_specks.Rd b/man/add_specks.Rd
index 8383940..bb538ad 100644
--- a/man/add_specks.Rd
+++ b/man/add_specks.Rd
@@ -4,11 +4,13 @@
 \alias{add_specks}
 \title{Add SPECKS to discrimination object}
 \usage{
-add_specks(discrimination, split = TRUE)
+add_specks(discrimination, group = c(), split = TRUE)
 }
 \arguments{
 \item{discrimination}{A discrimination with added propensities}
 
+\item{group}{A list of variable names to group by}
+
 \item{split}{A logical for if the metric should be calculated separately for
 the training/testing split. Defaults to TRUE.}
 }