diff --git a/README.Rmd b/README.Rmd
index 08de0dd..a3cf586 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -71,7 +71,7 @@ sim_data <- readRDS(url("https://zenodo.org/records/8433077/files/scLANE_sim_dat
 
 The PCA embeddings show us a pretty simple trajectory that's strongly correlated with the first principal component. 
 
-```{r plot-sims-pt, results='hold'}
+```{r plot-sims-pt, results='hold', message=FALSE}
 data.frame(sim_data@int_colData$reducedDims@listData$PCA[, 1:2]) %>% 
   mutate(pseudotime = sim_data$cell_time_normed) %>% 
   ggplot(aes(x = PC1, y = PC2, color = pseudotime)) + 
@@ -83,7 +83,7 @@ data.frame(sim_data@int_colData$reducedDims@listData$PCA[, 1:2]) %>%
 
 We also see that the data are not clustered by subject, which indicates that gene dynamics are mostly homogeneous across subjects. 
 
-```{r plot-sims-subj, results='hold'}
+```{r plot-sims-subj, results='hold', message=FALSE}
 data.frame(sim_data@int_colData$reducedDims@listData$PCA[, 1:2]) %>% 
   mutate(subject = sim_data$subject) %>% 
   ggplot(aes(x = PC1, y = PC2, color = subject)) + 
@@ -96,12 +96,8 @@ data.frame(sim_data@int_colData$reducedDims@listData$PCA[, 1:2]) %>%
 
 Since we have multi-subject data, we can use any of the three model modes to run our DE testing. We'll start with the simplest model, the GLM, then work our way through the other options in order of increasing complexity. We first prepare our inputs - a dataframe containing our cell ordering, a set of genes to build models for, and a vector of per-cell size factors to be used as offsets during estimation. In reality, it's usually unnecessary to fit a model for every single gene in a dataset, as trajectories are usually estimated using a subset of the entire set of genes (usually a few thousand most highly variable genes). For the purpose of demonstration, we'll select 50 genes each from the dynamic and non-dynamic populations. 
 
-{% note %}
-
 **Note:** In this case we're working with a single pseudotime lineage, though in real datasets several lineages often exist; in order to fit models for a subset of lineages simply remove the corresponding columns from the cell ordering dataframe passed as input to `testDynamic()`. 
 
-{% endnote %}
-
 ```{r prep-data}
 set.seed(312)
 gene_sample <- c(sample(rownames(sim_data)[rowData(sim_data)$geneStatus_overall == "Dynamic"], size = 50), 
@@ -175,12 +171,8 @@ scLANE_models_glmm <- testDynamic(sim_data,
                                   n.cores = 4)
 ```
 
-{% note %}
-
 **Note:** The GLMM mode is still under development, as we are working on further reducing runtime and increasing the odds of the underlying optimization process converging successfully. As such, updates will be frequent and functionality / results may shift slightly. 
 
-{% endnote %}
-
 Like the GLM mode, the GLMM mode uses a likelihood ratio test to compare the null & alternate models. We fit the two nested models using maximum likelihood estimation instead of [REML](https://en.wikipedia.org/wiki/Restricted_maximum_likelihood) in order to perform this test; the null model in this case is a negative binomial GLMM with a random intercept for each subject.  
 
 ```{r glmm-results}
diff --git a/tests/testthat/test_scLANE.R b/tests/testthat/test_scLANE.R
index 3cd0100..ca8e429 100644
--- a/tests/testthat/test_scLANE.R
+++ b/tests/testthat/test_scLANE.R
@@ -62,7 +62,7 @@ withr::with_output_sink(tempfile(), {
   # get results tables by interval
   glm_slope_test <- testSlope(glm_gene_stats)
   gee_slope_test <- testSlope(gee_gene_stats)
-  glmm_slope_test <- testSlope(glmm_gene_stats)
+  # glmm_slope_test <- testSlope(glmm_gene_stats)
   # run NB GAMs of varying structure
   gam_mod_bs <- nbGAM(expr = counts_test[, 1],
                       pt = pt_test,
@@ -294,13 +294,13 @@ test_that("getResultsDE() output", {
 test_that("testSlope() output", {
   expect_s3_class(glm_slope_test, "data.frame")
   expect_s3_class(gee_slope_test, "data.frame")
-  expect_s3_class(glmm_slope_test, "data.frame")
+  # expect_s3_class(glmm_slope_test, "data.frame")
   expect_gt(nrow(glm_slope_test), 0)
   expect_gt(nrow(gee_slope_test), 0)
-  expect_gt(nrow(glmm_slope_test), 0)
+  # expect_gt(nrow(glmm_slope_test), 0)
   expect_gt(sum(glm_slope_test$P_Val_Adj < 0.01, na.rm = TRUE), 0)
   expect_gt(sum(gee_slope_test$P_Val_Adj < 0.01, na.rm = TRUE), 0)
-  expect_gt(sum(glmm_slope_test$P_Val_Adj < 0.01, na.rm = TRUE), 0)
+  # expect_gt(sum(glmm_slope_test$P_Val_Adj < 0.01, na.rm = TRUE), 0)
 })
 
 test_that("nbGAM() output", {