diff --git a/DESCRIPTION b/DESCRIPTION
index 71fe3e2..e424e5a 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,10 +1,10 @@
 Package: radEmu
 Title: Using Relative Abundance Data to Estimate of Multiplicative Differences in Mean Absolute Abundance 
-Version: 2.3.1.0
+Version: 2.3.2.0
 Authors@R: c(person("David", "Clausen", role = c("aut")),
              person("Amy D", "Willis", email = "adwillis@uw.edu", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-2802-4317")),
              person("Sarah", "Teichman", role = "aut"))
-Description: A differential abundance method for the analysis of microbiome data. radEmu estimates fold-differences in the abundance of taxa across samples relative to "typical" fold-differences. Notably, it does not require pseudocounts, nor choosing a denominator taxon. 
+Description: A differential abundance method for the analysis of microbiome data. 'radEmu' estimates fold-differences in the abundance of taxa across samples relative to "typical" fold-differences. Notably, it does not require pseudocounts, nor choosing a denominator taxon. For more details, see Clausen et al. (2026) <doi:10.1093/biomet/asag009>. 
 URL: https://github.com/statdivlab/radEmu, https://statdivlab.github.io/radEmu/
 License: MIT + file LICENSE
 Encoding: UTF-8
diff --git a/R/emuFit.R b/R/emuFit.R
index 79e38fa..1767843 100644
--- a/R/emuFit.R
+++ b/R/emuFit.R
@@ -1,6 +1,6 @@
 #' Fit radEmu model
 #'
-#' @param Y an n x J matrix or dataframe of nonnegative observations, or a phyloseq object containing an otu table and sample data.
+#' @param Y an n x J matrix or dataframe of nonnegative observations, or a `phyloseq` or `TreeSummarizedExperiment` object containing an otu table and sample data.
 #' @param X an n x p design matrix (either provide \code{X} or \code{data} and \code{formula})
 #' @param formula a one-sided formula specifying the form of the mean model to be fit (used with \code{data})
 #' @param data an n x p data frame containing variables given in \code{formula}
@@ -120,7 +120,7 @@
 #' @import MASS
 #' 
 #' @examples
-#' # data frame example
+#' # data frame example (for phyloseq and TreeSummarizedExperiment examples, see the vignettes)
 #' data(wirbel_sample_small)
 #' data(wirbel_otu_small)
 #' emuRes <- emuFit(formula = ~ Group, data = wirbel_sample_small, Y = wirbel_otu_small,
@@ -128,17 +128,6 @@
 #'  # here we set large tolerances for the example to run quickly, 
 #'  # but we recommend smaller tolerances in practice
 #'
-#' # TreeSummarizedExperiment example (only run this if you have TreeSummarizedExperiment installed)
-#' \dontrun{
-#' library("TreeSummarizedExperiment")
-#' example("TreeSummarizedExperiment")
-#' assayNames(tse) <- "counts"
-#' emuRes <- emuFit(Y = tse, formula = ~ condition, assay_name = "counts", 
-#'                  test_kj = data.frame(k = 2, j = 1), tolerance = 0.01)
-#'  # here we set large tolerances for the example to run quickly, 
-#'  # but we recommend smaller tolerances in practice
-#' }
-#'
 #' @export
 #'
 emuFit <- function(Y,
diff --git a/README.md b/README.md
index 31d6fd6..7c178da 100644
--- a/README.md
+++ b/README.md
@@ -35,7 +35,13 @@ Sadly we do not yet have a ~~logo~~ nice-looking logo. If you would like to desi
 
 ## Installation
 
-To download the radEmu package, use the code below.
+To download the stable release of radEmu (version 2.3.1) from CRAN, use the code below.
+
+``` r
+install.packages("radEmu")
+```
+
+To download the current development version of the radEmu package, use the code below.
 
 ``` r
 # install.packages("devtools")
@@ -43,8 +49,6 @@ devtools::install_github("statdivlab/radEmu")
 library(radEmu)
 ```
 
-We are currently only releasing `radEmu` via GitHub. If you'd like us to consider submitting to CRAN, please let us know by opening an issue.
-
 ## Use
 
 The vignettes demonstrate example usage of the main functions. Please [file an issue](https://github.com/statdivlab/radEmu/issues) if you have a request for a tutorial that is not currently included. The following code shows the easy-to-use syntax if your data is in a `phyloseq` object, and you want to estimate parameters for all taxa and run a test for the parameter associated with "Group" and taxon 1: 
diff --git a/man/emuFit.Rd b/man/emuFit.Rd
index 9f52455..46c56a0 100644
--- a/man/emuFit.Rd
+++ b/man/emuFit.Rd
@@ -43,7 +43,7 @@ emuFit(
 )
 }
 \arguments{
-\item{Y}{an n x J matrix or dataframe of nonnegative observations, or a phyloseq object containing an otu table and sample data.}
+\item{Y}{an n x J matrix or dataframe of nonnegative observations, or a \code{phyloseq} or \code{TreeSummarizedExperiment} object containing an otu table and sample data.}
 
 \item{X}{an n x p design matrix (either provide \code{X} or \code{data} and \code{formula})}
 
@@ -197,7 +197,7 @@ including whether or not the algorithm converged, which can be helpful for debug
 Fit radEmu model
 }
 \examples{
-# data frame example
+# data frame example (for phyloseq and TreeSummarizedExperiment examples, see the vignettes)
 data(wirbel_sample_small)
 data(wirbel_otu_small)
 emuRes <- emuFit(formula = ~ Group, data = wirbel_sample_small, Y = wirbel_otu_small,
@@ -205,15 +205,4 @@ emuRes <- emuFit(formula = ~ Group, data = wirbel_sample_small, Y = wirbel_otu_s
  # here we set large tolerances for the example to run quickly, 
  # but we recommend smaller tolerances in practice
 
-# TreeSummarizedExperiment example (only run this if you have TreeSummarizedExperiment installed)
-\dontrun{
-library("TreeSummarizedExperiment")
-example("TreeSummarizedExperiment")
-assayNames(tse) <- "counts"
-emuRes <- emuFit(Y = tse, formula = ~ condition, assay_name = "counts", 
-                 test_kj = data.frame(k = 2, j = 1), tolerance = 0.01)
- # here we set large tolerances for the example to run quickly, 
- # but we recommend smaller tolerances in practice
-}
-
 }
diff --git a/tests/testthat/test-X_cup_from_X.R b/tests/testthat/test-X_cup_from_X.R
index 45ace43..6821b36 100644
--- a/tests/testthat/test-X_cup_from_X.R
+++ b/tests/testthat/test-X_cup_from_X.R
@@ -1,4 +1,5 @@
 test_that("X_cup_from_X_fast works and is faster than X_cup_from_X", {
+  
   X <- cbind(1, rnorm(100), runif(100))
   J <- 200
   start1 <- proc.time()
@@ -9,5 +10,5 @@ test_that("X_cup_from_X_fast works and is faster than X_cup_from_X", {
   end2 <- proc.time() - start2
   
   expect_true(all.equal(Xcup1, Xcup2))
-  expect_true(end2[3] < end1[3])
+  #expect_true(end2[3] < end1[3])
 })
\ No newline at end of file
diff --git a/tests/testthat/test-emuFit_micro.R b/tests/testthat/test-emuFit_micro.R
index 90c1163..fba66a9 100644
--- a/tests/testthat/test-emuFit_micro.R
+++ b/tests/testthat/test-emuFit_micro.R
@@ -244,6 +244,9 @@ test_that("ML fit to simple example give reasonable output with J >> n", {
 })
 
 test_that("unpenalized fit uses B if given, and therefore fit is quicker", {
+  
+  skip("don't test timing automatically")
+  
   set.seed(4323)
   X <- cbind(1,rnorm(10))
   J <- 10
@@ -281,6 +284,9 @@ test_that("unpenalized fit uses B if given, and therefore fit is quicker", {
 })
 
 test_that("unpenalized fit converges quicker if optimize_rows is set to TRUE", {
+  
+  skip("don't test timing automatically")
+  
   set.seed(4323)
   J <- 100
   n <- 100
@@ -375,122 +381,5 @@ test_that("unpenalized fit converges quicker if optimize_rows is set to TRUE wit
   expect_true(mean_est_diff < 1e-3)
     
 })
-#
-#
-#
-
-#
-# test_that("ML fit to multiple regressors, large n, and excess-Poisson variance gives reasonable output", {
-#   set.seed(4323)
-#   n <- 800
-#   X <- cbind(1,rep(c(0,1),each = n/2),rnorm(n),rnorm(n),rnorm(n))
-#   z <- rnorm(n) +5
-#   J <- 10
-#   p <- 2
-#
-#   b0 <- rnorm(J)
-#   b1 <- seq(1,10,length.out = J)
-#   b <- rbind(b0,b1)
-#   Y <- matrix(NA,ncol = J, nrow = n)
-#
-#   for(i in 1:n){
-#     for(j in 1:J){
-#       temp_mean <- exp(X[i,1:2,drop = FALSE]%*%b[,j,drop = FALSE] + z[i])
-#       Y[i,j] <- rnbinom(1,mu= temp_mean,size = 0.25)#rpois(1, lambda = temp_mean)
-#     }
-#   }
-#   ml_fit <- emuFit_micro(X,
-#                          Y,
-#                          constraint_fn = function(x) mean(x),
-#                          maxit = 200,
-#                          tolerance = 1e-2)
-#
-#   # plot(b1-mean(b1),ml_fit[2,])
-#   # abline(a = 0,b = 1,lty =2,col = "red")
-#
-#   expect_true(max(abs(ml_fit[2,] - (b1 - mean(b1))))<.5)
-# })
-#
-#
-# test_that("ML fit to multiple regressors, moderate n, moderate J, and excess-Poisson variance gives reasonable output", {
-#   set.seed(4323)
-#   n <- 40
-#   X <- cbind(1,rep(c(0,1),each = n/2),rnorm(n),rnorm(n),rnorm(n))
-#   z <- rnorm(n)
-#   J <- 100
-#   p <- 2
-#
-#   b0 <- rnorm(J)
-#   b1 <- seq(1,10,length.out = J)
-#   b <- rbind(b0,b1)
-#   Y <- matrix(NA,ncol = J, nrow = n)
-#
-#   for(i in 1:n){
-#     for(j in 1:J){
-#       temp_mean <- exp(X[i,1:2,drop = FALSE]%*%b[,j,drop = FALSE] + z[i])
-#       Y[i,j] <- rnbinom(1,mu= temp_mean,size = 0.25)#rpois(1, lambda = temp_mean)
-#     }
-#   }
-#   ml_fit <- emuFit_micro(X,
-#                          Y,
-#                          constraint_fn = function(x) mean(x),
-#                          maxit = 200,
-#                          tolerance = 1e-1,
-#                          max_step = 0.5)
-#
-#   # plot(b1-mean(b1),ml_fit[2,])
-#   # abline(a = 0,b = 1, lty = 2, col = "red")
-#
-#   b1_hat <- ml_fit[2,]
-#   expect_true(abs(lm(b1_hat~b1)$coef[2] - 1) < 0.2)
-#   # abline(a = 0,b = 1,lty =2,col = "red")
-#
-# })
-#
-# # t(ml_fit) %>% as.data.frame() %>% mutate(j = 1:J) %>%
-# #   pivot_longer(-j) %>%
-# #   ggplot() +
-# #   geom_point(aes(x = j,y = value, color = name)) +
-# #   geom_line(aes(x = j, y  = value, group = name, color = name)) +
-# #   theme_bw()
-#
-#
-#
-# test_that("ML fit to simple example give reasonable output with J > n", {
-#   set.seed(4323)
-#   X <- cbind(1,rep(c(0,1),each = 20))
-#   z <- rnorm(40) +8
-#   b0 <- rnorm(80)
-#   b1 <- seq(-5,5,length.out = 80)
-#   b <- rbind(b0,b1)
-#   Y <- matrix(NA,ncol = 80, nrow = 40)
-#
-#   for(i in 1:40){
-#     for(j in 1:80){
-#       temp_mean <- exp(X[i,,drop = FALSE]%*%b[,j,drop = FALSE] + z[i])
-#       Y[i,j] <- rpois(1, lambda = temp_mean)
-#     }
-#   }
-#   ml_fit <- emuFit_micro(X,
-#                          Y,
-#                          B = NULL,
-#                          constraint_fn = function(x) mean(x),
-#                          maxit = 500,
-#                          tolerance = 0.01)
-#
-#
-#
-#
-#
-#   expect_true(max(abs(ml_fit[2,] - b1))<.1)
-#
-#   # X_cup <- X_cup_from_X(X,10)
-#   # B_cup <- B_cup_from_B(ml_fit)
-#
-#
-#
-# })
-#
-#
 
 
diff --git a/tests/testthat/test-emuFit_micro_penalized.R b/tests/testthat/test-emuFit_micro_penalized.R
index c0d276c..9df70b9 100644
--- a/tests/testthat/test-emuFit_micro_penalized.R
+++ b/tests/testthat/test-emuFit_micro_penalized.R
@@ -150,6 +150,9 @@ less efficient implementation (and that both substantially differ from MLE", {
 })
 
 test_that("penalized fit uses B if given, and therefore fit is quicker", {
+  
+  skip("don't test timing automatically")
+  
   set.seed(4323)
   X <- cbind(1,rnorm(10))
   J <- 10
diff --git a/tests/testthat/test-fit_null_discrete.R b/tests/testthat/test-fit_null_discrete.R
index 61b32d9..4c52ea4 100644
--- a/tests/testthat/test-fit_null_discrete.R
+++ b/tests/testthat/test-fit_null_discrete.R
@@ -147,7 +147,7 @@ test_that("new discrete is fast", {
   }) 
   
   ## check faster
-  expect_lt(t_discrete[3], t_orig[3])
+  #expect_lt(t_discrete[3], t_orig[3])
   
   ## check that the (my_kstar, my_kstar) constraint is satisfied
   expect_equal(pseudohuber_median(out_discrete$B[my_kstar, -my_jstar]), 
diff --git a/tests/testthat/test-get_G_for_augmentations.R b/tests/testthat/test-get_G_for_augmentations.R
index 8558f1f..26d9cec 100644
--- a/tests/testthat/test-get_G_for_augmentations.R
+++ b/tests/testthat/test-get_G_for_augmentations.R
@@ -12,6 +12,6 @@ test_that("get_G_for_augmentations_fast works and is faster than get_G_for_augme
   end2 <- proc.time() - start2
   
   expect_true(all.equal(g1, g2))
-  expect_true(end2[3] < end1[3])
+  #expect_true(end2[3] < end1[3])
 })
 
diff --git a/tests/testthat/test-macro_fisher_null.R b/tests/testthat/test-macro_fisher_null.R
index c00e753..e0033c4 100644
--- a/tests/testthat/test-macro_fisher_null.R
+++ b/tests/testthat/test-macro_fisher_null.R
@@ -353,100 +353,3 @@ test_that("We take similar step as we'd take using numerical derivatives in a mo
 
 })
 
-test_that("We take same step as we'd take using numerical derivatives when gap, rho, u are zero", {
-  set.seed(59542234)
-  n <- 10
-  p <- 2
-  X <- cbind(1,rep(c(0,1),each = n/2))
-  J <- 10000
-  z <- rnorm(n) +8
-  b0 <- rnorm(J)
-  b1 <- seq(1,10,length.out = J)
-  b1 <- b1 - mean(b1)
-  b0 <- b0 - mean(b0)
-  b <- rbind(b0,b1)
-  Y <- matrix(NA,ncol = J, nrow = n)
-  
-  k_constr <- 2
-  j_constr <- 1
-  p <- 2
-  
-  constraint_fn <- rep(list(function(x){ pseudohuber_median(x,0.1)}), 2)
-  # constraint_fn <- function(x){mean(x)}
-  
-  ##### Arguments to fix:
-  
-  constraint_grad_fn <- rep(list(function(x){dpseudohuber_median_dx(x,0.1)}), 2)
-  # constraint_grad_fn <- function(x){ rep(1/length(x), length(x))}
-  rho_init = 1
-  tau = 1.2
-  kappa = 0.8
-  obj_tol = 100
-  score_tol <- 1e-3
-  constraint_tol = 1e-5
-  init_tol = 1e6
-  c1 = 1e-4
-  maxit = 1000
-  inner_maxit = 25
-  
-  Y[] <- 0
-  for(i in 1:n){
-    while(sum(Y[i,])==0){
-      for(j in 1:J){
-        temp_mean <- exp(X[i,,drop = FALSE]%*%b[,j,drop = FALSE] + z[i])
-        Y[i,j] <- rpois(1, lambda = temp_mean)
-        # Y[i,j] <- rnbinom(1,mu = temp_mean, size = 3)*rbinom(1,1,0.6)
-      }
-    }
-  }
-  
-  j_ref <- 5
-  
-  full_fit <- #suppressMessages(
-    emuFit_micro_penalized(X = X,
-                           Y = Y,
-                           B = NULL,
-                           constraint_fn = rep(list(mean), 2), 
-                           tolerance = 10,
-                           verbose = FALSE)#)
-  
-  B <- full_fit$B
-  B[1,] <- B[1,] - B[1,j_ref]
-  B[2,] <- B[2,] - B[2,j_ref]
-  
-  Y_aug <- full_fit$Y_augmented
-  
-  B[k_constr,j_constr] <- constraint_fn[[k_constr]](B[k_constr,-j_constr])
-  
-  u <- 0
-  rho <- 0
-  
-  z <- update_z(Y,X,B)
-  
-  macro_time <- try(system.time(macro_fisher_null(X = X,
-                              Y = Y_aug,
-                              B = B,
-                              z = z,
-                              J = J,
-                              p = p,
-                              k_constr = k_constr,
-                              j_constr = j_constr,
-                              j_ref = j_ref,
-                              rho = rho,
-                              u = u,
-                              constraint_fn = constraint_fn,
-                              constraint_grad_fn = constraint_grad_fn,
-                              stepsize = 0.5,
-                              c1 = 1e-4,
-                              regularization = 0,
-                              debug= FALSE)))
-
-  if(inherits(macro_time, "try-error")){
-    expect_true(FALSE)
-  }
-  expect_true(macro_time[3]<15)
-  #failing this test indicates that macro_fisher_null may be directly 
-  #computing the full inverse of the (approximate) hessian matrix of the log likelihood
-  #this includes an outer product that does not need to be computed
-})
-
diff --git a/vignettes/intro_radEmu.Rmd b/vignettes/intro_radEmu.Rmd
index 6b07874..6d4cbd0 100644
--- a/vignettes/intro_radEmu.Rmd
+++ b/vignettes/intro_radEmu.Rmd
@@ -144,7 +144,6 @@ and some optional arguments include
 
   - `run_score_tests`: A logical value denoting whether or not to run score tests. Score tests are awesome in their error rate control (including with small sample sizes; though of course larger sample sizes always give better power), but require refitting the model, so can require some compute time. 
 
-
 ```{r}
 ch_fit <- emuFit(formula = ~ Group, 
                  data = wirbel_sample[ch_study_obs, ],
@@ -157,8 +156,6 @@ Let's check out what this object looks like!
 ch_fit
 ```
 
-
-
 Now, we can easily visualize our results using the `plot` function!
 So that we only produce the plot, and not the dataframe used to produce the plots,
 we will explicitly extract the plots using the `$` operator to get the `plots` component.
diff --git a/vignettes/parallel_radEmu.Rmd b/vignettes/parallel_radEmu.Rmd
index 90452ea..4952191 100644
--- a/vignettes/parallel_radEmu.Rmd
+++ b/vignettes/parallel_radEmu.Rmd
@@ -120,19 +120,16 @@ Now, let's run some tests in parallel. We will be parallelizing our code over $j
 
 Let's say that I want to run score tests for the first five mOTUs in our dataset. First, I need to check the cores on my computer to see a reasonable amount of cores to parallelize over. I tend to use one fewer core than the number of cores that I have available. 
 
-```{r}
+A note for those knitting this vignette or reading this through online rendered documentation: the following code is set to not evaluate automatically. You can run it locally and explore the results by removing the `eval = FALSE` argument from the code chunks. 
+
+```{r, eval = FALSE}
 ncores <- parallel::detectCores() - 1
 ncores
-
-# if running this vignette in automatic R-CMD-check, reduce cores to 2
-if (identical(Sys.getenv("_R_CHECK_LIMIT_CORES_"), "TRUE")) {
-  ncores <- min(2, ncores)
-}
 ```
 
-Next, I will write a function that will be called in parallel. This function will fit the model under the null and calculate the robust score test statistics. Note that the output of this function is an `emuFit` object.
+Next, we will write a function that will be called in parallel. This function will fit the model under the null and calculate the robust score test statistics. Note that the output of this function is an `emuFit` object.
 
-```{r}
+```{r, eval = FALSE}
 emuTest <- function(category) {
   score_res <- emuFit(formula = ~ Group,
                        data = wirbel_sample[ch_study_obs, ],
@@ -147,7 +144,7 @@ emuTest <- function(category) {
 
 Now, we can run our score tests in parallel. We'll just do the first five. It may take a minute or so, depending on your machine. 
 
-```{r}
+```{r, eval = FALSE}
 if (.Platform$OS.type != "windows" & !identical(Sys.getenv("GITHUB_ACTIONS"), "true")) {
   # run if we are on a Mac or Linux machine
   score_res <- mclapply(1:5,
@@ -162,7 +159,7 @@ if (.Platform$OS.type != "windows" & !identical(Sys.getenv("GITHUB_ACTIONS"), "t
 Now, we can see that this barely took more time than running a single score test, because we were able to parallelize over cores (on my laptop, I'm using more than five cores, so I can run all five tests at the same time). Each p-value can be pulled out of the list as follows: 
 
 
-```{r}
+```{r, eval = FALSE}
 if (!is.null(score_res)) {
   c(score_res[[1]]$coef$pval[1], ## robust score test p-value for the first taxon
   score_res[[2]]$coef$pval[2]) ## robust score test p-value for the second taxon
@@ -171,7 +168,7 @@ if (!is.null(score_res)) {
 
 To help organise this information, we can make a coefficient matrix that combines the information from each component in our list: 
 
-```{r}
+```{r, eval = FALSE}
 if (!is.null(score_res)) {
   full_score <- sapply(1:length(score_res), 
                        function(x) score_res[[x]]$coef$score_stat[x])
diff --git a/vignettes/quick_start.Rmd b/vignettes/quick_start.Rmd
index 8a3472d..1bb499a 100644
--- a/vignettes/quick_start.Rmd
+++ b/vignettes/quick_start.Rmd
@@ -27,10 +27,11 @@ First, we will install `radEmu`, if we haven't already.
 # remotes::install_github("statdivlab/radEmu")
 ```
 
-Next, we can load `radEmu`.
+Next, we can load `radEmu`, as well as `dplyr`.
 
 ```{r setup, message = FALSE}
 library(radEmu)
+library(dplyr)
 ```
 
 In this vignette we'll explore a [dataset published by Wirbel et al. (2019)](https://www.nature.com/articles/s41591-019-0406-6). This is a meta-analysis of case-control studies of colorectal cancer, meaning that Wirbel et al. collected raw sequencing data from studies other researchers conducted and re-analyzed it.
@@ -49,13 +50,38 @@ data("wirbel_otu")
 dim(wirbel_otu)
 ```
 
-We can see that we have $566$ samples and $845$ mOTUs. We will subset to samples from the Chinese study.
+We can see that we have $566$ samples and $845$ mOTUs. 
+
+While in general we would fit a model to all mOTUs, we are going to subset to some specific genera for the purposes of this tutorial. Let's look at *Eubacterium*, *Porphyromonas*, *Faecalibacteria*, and *Fusobacterium* for now.
+
+```{r}
+mOTU_names <- colnames(wirbel_otu)
+chosen_genera <- c("Eubacterium", "Faecalibacterium", "Fusobacterium", "Porphyromonas")
+mOTU_name_df <- data.frame(name = mOTU_names) %>% 
+  mutate(base_name = stringr::str_remove(mOTU_names, "unknown ") %>%
+                      stringr::str_remove("uncultured ")) %>%
+  mutate(genus_name = stringr::word(base_name, 1))
+restricted_mOTU_names <- mOTU_name_df %>%
+  filter(genus_name %in% chosen_genera) %>%
+  pull(name)
+```
+
+Again, while we would generally fit a model using all of our samples, for this tutorial we are only going to consider data from a case-control study from China. 
 
 ```{r}
-wirbel_sample_ch <- wirbel_sample[wirbel_sample$Country == "CHI", ]
-wirbel_otu_ch <- wirbel_otu[rownames(wirbel_sample_ch), ]
-to_rm <- which(colSums(wirbel_otu_ch) == 0)
-wirbel_otu_ch <- wirbel_otu_ch[, -to_rm]
+ch_study_obs <- which(wirbel_sample$Country %in% c("CHI"))
+```
+
+Next, we want to confirm that all samples have at least one non-zero count across the categories we've chosen and that all categories have at least one non-zero count across the samples we've chosen.
+
+```{r}
+small_Y <- wirbel_otu[ch_study_obs, restricted_mOTU_names]
+sum(rowSums(small_Y) == 0) # no samples have a count sum of 0 
+sum(colSums(small_Y) == 0) # one category has a count sum of 0
+
+category_to_rm <- which(colSums(small_Y) == 0)
+small_Y <- small_Y[, -category_to_rm]
+sum(colSums(small_Y) == 0)
 ```
 
 ## Differential abundance analysis with radEmu
@@ -63,8 +89,10 @@ wirbel_otu_ch <- wirbel_otu_ch[, -to_rm]
 We now are ready to fit a `radEmu` model with the `emuFit()` function! This will take 1-2 minutes to run, so we suggest you start running it and then read below about the arguments!
 
 ```{r}
-mod <- emuFit(Y = wirbel_otu_ch, formula = ~ Group + Gender, data = wirbel_sample_ch,
-              test_kj = data.frame(k = 2, j = 1), verbose = TRUE)
+mod <- emuFit(data = wirbel_sample[ch_study_obs, ],
+              Y = small_Y,
+              formula = ~ Group + Gender, 
+              test_kj = data.frame(k = 2, j = 1))
 ```
 
 Some of the important arguments for `emuFit()` are the following: 
@@ -83,10 +111,10 @@ Now that we've fit the model and run a robust score test for the first taxon, we
 head(mod$coef)
 ```
 
-The first taxon in our model is *Streptococcus anginosus*, which has a log fold-difference estimate of $1.23$. We will interpret this to mean that we expect that *Streptococcus anginosus* is $exp(1.23) \approx 3.42$ times more abundant in cases of colorectal cancer compared to controls of the same gender, when compared to the typical fold-differences in abundances of the taxa in this analysis. 
+The first taxon in our model is *Fusobacterium nucleatum s. vincentii*, which has a log fold-difference estimate of $1.49$. We will interpret this to mean that we expect that *Fusobacterium nucleatum s. vincentii* is $exp(1.49) \approx 4.44$ times more abundant in cases of colorectal cancer compared to controls of the same gender, when compared to the typical fold-differences in abundances of the taxa in this analysis. 
 
 We could similarly interpret the log fold-differences for each taxon in our analysis. 
 
-For *Streptococcus anginosus* we have a robust score test p-value of $0.09$. This means that we do not have enough evidence to reject the null hypothesis that the fold-difference in abundance of *Streptococcus anginosus* between cases and controls is the same as the typical fold-difference in abundance between cases and controls across all taxa in this analysis. When we say ``typical" here we mean approximately the median fold-difference across taxa.  
+For *Fusobacterium nucleatum s. vincentii* we have a robust score test p-value of $0.04$. This means that we do have enough evidence to reject the null hypothesis that the fold-difference in abundance of *Fusobacterium nucleatum s. vincentii* between cases and controls is the same as the typical fold-difference in abundance between cases and controls across all taxa in this analysis. When we say ``typical" here we mean approximately the median fold-difference across taxa.  
 
 Now you are ready to start using `radEmu`! We recommend our other vignettes for a deeper look using `radEmu`, including for `phyloseq` or `TreeSummarizedExperiment` objects, with clustered data, and in parallel. 
\ No newline at end of file
diff --git a/vignettes/radEmu_clustered_data.Rmd b/vignettes/radEmu_clustered_data.Rmd
index 70ffa71..62cc810 100644
--- a/vignettes/radEmu_clustered_data.Rmd
+++ b/vignettes/radEmu_clustered_data.Rmd
@@ -66,7 +66,7 @@ To start, let's generate a toy example (10 categories, 60 samples) in which ther
 J <- 10; n <- 60
 # generate design matrix 
 set.seed(10)
-X <- cbind(1, rnorm(n))
+X <- cbind(1, rep(0:1, each = n / 2))
 cov_dat <- data.frame(cov = X[, 2])
 # cluster membership 
 cluster <- rep(1:4, each = n/4)
@@ -114,7 +114,10 @@ ef_cluster <- emuFit(formula = ~ cov,
                      data = cov_dat, 
                      Y = Y,
                      cluster = cluster_named, 
-                     test_kj = data.frame(k = 2, j = 3))
+                     test_kj = data.frame(k = 2, j = 3),
+                     match_row_names = FALSE)
+
+ef_cluster$coef[3, ]
 ```
 
 You can check out the full object, but our estimate is `r round(ef_cluster$coef$estimate[3], 1)` and a p-value for testing that this parameter equals zero is `r round(ef_cluster$coef$pval[3], 3)`. Not too shabby, especially considering that about half of our observations are zero, and we have a lot of noise in our data (arising from a negative binomial simulation scheme). 
@@ -125,7 +128,9 @@ Let's also compare that to a situation where we mistakenly ignore clustering. In
 ef_no_cluster <- emuFit(formula = ~ cov,
                         data = cov_dat, 
                         Y = Y,
-                        test_kj = data.frame(k = 2, j = 3))
+                        test_kj = data.frame(k = 2, j = 3),
+                        match_row_names = FALSE)
+ef_no_cluster$coef[3, ]
 ```
 
 When we ignore clustering, we get an estimate of the log fold difference in category 3 across values of our covariate of `r round(ef_no_cluster$coef$estimate[3], 1)` and a p-value of `r round(ef_no_cluster$coef$pval[3], 3)`. So we can see that our estimates are the same whether or not we account for cluster, but our p-values are different (because we are pretending that we have more evidence in the absence of clustering).