diff --git a/.Rbuildignore b/.Rbuildignore
index 5f3d199..612dedc 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -3,3 +3,5 @@
 ^LICENSE\.md$
 ^README\.Rmd$
 ^\.github$
+^doc$
+^Meta$
diff --git a/.gitignore b/.gitignore
index c833a2c..7a7155e 100644
--- a/.gitignore
+++ b/.gitignore
@@ -3,3 +3,5 @@
 .RData
 .Ruserdata
 inst/doc
+/doc/
+/Meta/
diff --git a/DESCRIPTION b/DESCRIPTION
index 756ed99..05172c7 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: niceday
 Title: Nonparametric estimates and confidence intervals for fold-difference parameters
-Version: 0.1.0
+Version: 1.0.0
 Authors@R: c(person("Grant", "Hopkins", 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"))
@@ -10,17 +10,20 @@ Encoding: UTF-8
 LazyData: true
 RoxygenNote: 7.3.2
 Depends: 
-  gam,
-  MASS,
-  glmnet,
-  hal9001,
-  ranger,
-  xgboost,
-  withr,
-  SuperLearner
+    gam,
+    MASS,
+    glmnet,
+    hal9001,
+    ranger,
+    xgboost,
+    withr,
+    SuperLearner,
+    R (>= 3.5)
 Suggests:
     knitr,
     rmarkdown,
+    dplyr,
+    ggplot2,
     testthat (>= 3.0.0)
 Config/testthat/edition: 3
 VignetteBuilder: knitr
diff --git a/NAMESPACE b/NAMESPACE
index 544c099..a1e40f6 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,4 +1,5 @@
 # Generated by roxygen2: do not edit by hand
 
+S3method(print,ndFit)
 export(est_nuis)
 export(ndFit)
diff --git a/R/EcoZUR_count.R b/R/EcoZUR_count.R
new file mode 100644
index 0000000..0b7e436
--- /dev/null
+++ b/R/EcoZUR_count.R
@@ -0,0 +1,7 @@
+#' Count data from EcoZUR study, example data used in examples and vignettes
+#'
+#' A dataset containing the count matrix from the EcoZur study.
+#'
+#' @format A data frame with with 358 rows and 229 columns. Each row is an observation and each column is a microbial taxon.
+#' @source Jesser et al. (2023). "Why are so many enteric pathogen infections asymptomatic? Pathogen and gut microbiome characteristics associated with diarrhea symptoms and carriage of diarrheagenic E. coli in northern Ecuador." Gut Microbes 15(2).
+"EcoZUR_count"
diff --git a/R/EcoZUR_meta.R b/R/EcoZUR_meta.R
new file mode 100644
index 0000000..d03b636
--- /dev/null
+++ b/R/EcoZUR_meta.R
@@ -0,0 +1,14 @@
+#' Metadata from EcoZUR study, example data used in examples and vignettes
+#'
+#' A dataset containing metadata from the EcoZur study.
+#'
+#' @format A data frame with 358 rows and 5 variables:
+#' \describe{
+#'   \item{Diarrhea}{Indicator of whether sample is from a participant with diarrhea (Case) or a participant without diarrhea (control) (character)}
+#'   \item{sex}{Sex of participant (character)}
+#'   \item{age_months}{Age in months of participant (integer)}
+#'   \item{urban}{How urban of a population the sample is from, from 0 to 5 (integer)}
+#'   \item{site}{Study site (character)}
+#' }
+#' @source Jesser et al. (2023). "Why are so many enteric pathogen infections asymptomatic? Pathogen and gut microbiome characteristics associated with diarrhea symptoms and carriage of diarrheagenic E. coli in northern Ecuador." Gut Microbes 15(2).
+"EcoZUR_meta"
diff --git a/R/est_nuis.R b/R/est_nuis.R
index abb6dff..d7e24dc 100644
--- a/R/est_nuis.R
+++ b/R/est_nuis.R
@@ -5,19 +5,26 @@
 #' @param X Matrix of covariates to adjust for.
 #' @param num_crossval_folds Number of folds for cross-validation. Default is `10`.
 #' @param num_crossfit_folds Number of folds for cross-fitting. Default is `10`.
-#' @param gtrunc ?? Default is `min(0.05, 5 / sqrt(NROW(W)) / log(NROW(W)))`.
-#' @param sl.lib.pi Libraries used to estimate ???. Default is `c("SL.mean", "SL.glm.binom", "SL.glmnet.binom", "SL.gam.binom", "SL.ranger.binom", "SL.hal9001.binom")`.
-#' @param sl.lib.m Libraries used to estimate ???. Default is `c("SL.mean", "SL.glm.pois", "SL.glmnet.pois", "SL.ranger.pois")`.
-#' @param sl.lib.q Libraries used to estimate ???. Default is the input to `sl.lib.pi`.
+#' @param gtrunc Truncation parameter, bounding the estimated propensity scores away from `0` and `1` by `gtrunc`. Default is `min(0.05, 5 / sqrt(NROW(W)) / log(NROW(W)))`.
+#' @param sl.lib.pi Libraries used to estimate the propensity score nuisance function. Default is `c("SL.mean", "SL.glm.binom", "SL.glmnet.binom", "SL.gam.binom", "SL.ranger.binom", "SL.hal9001.binom")`.
+#' @param sl.lib.m Libraries used to estimate conditional mean mu_j nuisance functions. Default is `c("SL.mean", "SL.glm.pois", "SL.glmnet.pois", "SL.ranger.pois")`.
+#' @param sl.lib.q Libraries used to estimate conditional probabilities on nonzero observations. Default is the input to `sl.lib.pi`.
 #' @param allow_warnings Allow warnings? Default is `TRUE`.
 #' @param enforce_pos_reg Should estimates of \eqn{E[W_j|A=a,X]} to forced to be strictly positive? Default is `FALSE`.
 #' @param verbose Do you want to receive updates as this function runs? Default is `TRUE`.
 #'
-#' @return A list containing elements `noadjust`, `adjust`, `nuis`, `cf_nuis`, and `variance`. `noadjust` gives unadjusted parameter estimates.
-#' `adjust` gives covariate adjusted parameter estimates. `nuis` gives estimates of nuisance parameters. Add a description of the rest here!
+#' @return A list of estimated nuisance functions and additional diagnostic information.
 #'
 #' @examples
-#' # add example here!
+#' data(EcoZUR_meta)
+#' data(EcoZUR_count)
+#' nuis <- est_nuis(W = EcoZUR_count[, 1:50], # consider only the first 50 taxa to run quickly
+#'                  A = (EcoZUR_meta[, "Diarrhea"] == "Case") * 1,
+#'                  X = model.matrix(~ sex + age_months, EcoZUR_meta),
+#'                  num_crossval_folds = 2,
+#'                  num_crossfit_folds = 2,
+#'                  sl.lib.pi = c("SL.mean"), # choosing single learner for the example to run quickly,
+#'                  sl.lib.m = c("SL.mean"))  # in practice would use other options as well
 #'
 #' @export
 est_nuis <- function(W,
diff --git a/R/ndFit.R b/R/ndFit.R
index baa57e2..9d14ca1 100644
--- a/R/ndFit.R
+++ b/R/ndFit.R
@@ -1,34 +1,46 @@
 #' Fit niceday model
 #'
 #' @param W Matrix of multivariate outcome variables.
-#' @param A Binary covariate of interest.
-#' @param X Matrix of covariates to adjust for.
+#' @param A Formula indicating binary covariate of interest.
+#' @param X Formula indicating covariates to adjust for. Default is no adjustment, `~ 1`
+#' @param data Data frame containing any covariates in formulas for `A` and `X`.
 #' @param alpha Desired asymptotic type I error rate. Default is `0.05`.
-#' @param uniform_CI ?? Default is `TRUE`.
+#' @param uniform_CI Generate simultaneous confidence intervals across al categories? Default is `TRUE`.
 #' @param d Parameter for smoothed median centering. Default is `0.1`.
-#' @param gtrunc ?? Default is `0.01`.
+#' @param gtrunc Truncation parameter passed to `est_nuis`, bounding the estimated propensity scores away from `0` and `1` by `gtrunc`. Default is `0.01`.
 #' @param bs_rep Number of replicates for bootstrap sampling of uniform confidence intervals. Default is `1e5`.
 #' @param num_crossval_folds Number of folds for cross-validation. Default is `10`.
 #' @param num_crossfit_folds Number of folds for cross-fitting. Default is `10`.
 #' @param enforce_pos_reg Should estimates of \eqn{E[W_j|A=a,X]} to forced to be strictly positive? Default is `FALSE`.
-#' @param sl.lib.pi Libraries used to estimate ???. Default is `c("SL.mean", "SL.lm", "SL.glm.binom", "SL.xgboost.binom")`.
-#' @param sl.lib.m Libraries used to estimate ???. Default is `c("SL.mean", "SL.lm", "SL.glm.qpois", "SL.xgboost.pois")`.
-#' @param sl.lib.q Libraries used to estimate ???. Default is the input to `sl.lib.pi`.
-#' @param adjust_covariates Should covariates be adjusted for in the estimator? Default is `TRUE`.
-#' @param nuis ??. Default is `NULL`.
+#' @param sl.lib.pi Libraries used to estimate the propensity score nuisance function. Default is `c("SL.mean", "SL.lm", "SL.glm.binom", "SL.xgboost.binom")`.
+#' @param sl.lib.m Libraries used to estimate conditional mean mu_j nuisance functions. Default is `c("SL.mean", "SL.lm", "SL.glm.qpois", "SL.xgboost.pois")`.
+#' @param sl.lib.q Libraries used to estimate conditional probabilities on nonzero observations. Default is the input to `sl.lib.pi`.
+#' @param nuis Optional list of estimated nuisance functions (from function `est_nuis()`). Default is `NULL`.
 #' @param verbose Do you want to receive updates as this function runs? Default is `FALSE`.
 #' @param cross_fit Should cross-fitting be run? Default is `TRUE`.
 #'
-#' @return A list containing elements `noadjust`, `adjust`, `nuis`, `cf_nuis`, and `variance`. `noadjust` gives unadjusted parameter estimates.
-#' `adjust` gives covariate adjusted parameter estimates. `nuis` gives estimates of nuisance parameters. Add a description of the rest here!
+#' @return A list containing elements `coef`, `nuisances`, `crossfit_nuisances`, `variance`, and `call`. `coef` contains estimates,
+#' standard errors, and confidence intervals for all categories. `nuisances` and `crossfit_nuisances` contain estimated nuisance
+#' parameters and nuisance parameters from the crossfitting procedure (included when adjustment covariates are included).
+#' `variance` contains the estimated covariance matrix. `call` contains the call to `ndFit()`.
 #'
 #' @examples
-#' # add example here!
+#' data(EcoZUR_meta)
+#' data(EcoZUR_count)
+#' ndFit(W = EcoZUR_count[, 1:50], # consider only the first 50 taxa to run quickly
+#'       data = EcoZUR_meta,
+#'       A = ~ Diarrhea,
+#'       X = ~ sex + age_months,
+#'       num_crossval_folds = 2, # use more folds in practice
+#'       num_crossfit_folds = 2, # for cross validation and cross fitting
+#'       sl.lib.pi = c("SL.mean"), # choosing single learner for the example to run quickly,
+#'       sl.lib.m = c("SL.mean"))  # in practice would use other options as well
 #'
 #' @export
 ndFit <- function(W,
                   A,
-                  X,
+                  X = ~ 1,
+                  data,
                   alpha = 0.05,
                   uniform_CI = TRUE,
                   d = 0.1,
@@ -46,22 +58,33 @@ ndFit <- function(W,
                                "SL.glm.qpois",
                                "SL.xgboost.pois"),
                   sl.lib.q = sl.lib.pi,
-                  adjust_covariates = TRUE,
                   nuis = NULL,
                   verbose = FALSE,
                   cross_fit = TRUE) {
 
+  call <- match.call(expand.dots = FALSE)
+
   # perform checks to ensure valid data input
   if (!is.data.frame(W)) {
     stop("W must be an nxJ data.frame")
   }
 
-  if (!is.vector(A)) {
-    stop("A must be a length n vector")
+  if (!inherits(A, "formula")) {
+    stop("A must be a formula")
+  }
+
+  # transform A from formula into a vector
+  A <- stats::model.matrix(A, data)
+  if (ncol(A) > 2) {
+    stop("The `A` formula must only include a single covariate")
+  }
+  A <- A[, 2]
+  if (any(!(A %in% c(0, 1)))) {
+    stop("The covariate A must be binary")
   }
 
-  if (length(A) != nrow(W)) {
-    stop("The length of A does not match the number of rows of W")
+  if (nrow(data) != nrow(W)) {
+    stop("The number of rows in `data` does not match the number of rows of `W`")
   }
 
   if (any(colMeans(W[A == 1, ]) == 0) | any(colMeans(W[A == 0, ]) == 0)) {
@@ -73,11 +96,12 @@ ndFit <- function(W,
     stop("One of the subgroups of interest (0 or 1) is entirely absent.")
   }
 
-  if (adjust_covariates) {
-    if (!(is.data.frame(X) | is.matrix(X))) {
-      stop("User requested covariate adjustment. The input X must be ",
-           "provided as a matrix or a data.frame object.")
-    }
+  # transform X formula into design matrix
+  X <- stats::model.matrix(X, data)
+  if (ncol(X) == 1) {
+    adjust_covariates = FALSE
+  } else {
+    adjust_covariates = TRUE
   }
 
   if (!is.null(nuis)) {
@@ -180,7 +204,7 @@ ndFit <- function(W,
     upper_g_log_noadj_sim <- NA
   }
 
-  results_simple <- list(ABC = data.frame(taxon = 1:J,
+  results_simple <- list(Psi1 = data.frame(category = 1:J,
                                           est = psi_hat_ABC_simp,
                                           se = se_hat_psi_hat_ABC_simp,
                                           lower_marg = lower_log_noadj_marg,
@@ -188,7 +212,7 @@ ndFit <- function(W,
                                           lower_sim = lower_log_noadj_sim,
                                           upper_sim = upper_log_noadj_sim),
 
-                         ABC_g = data.frame(taxon = 1:J,
+                         Psi1g = data.frame(category = 1:J,
                                             est = psi_hat_ABC_g_simp,
                                             se = se_hat_psi_hat_ABC_g_simp,
                                             lower_marg = lower_g_log_noadj_marg,
@@ -227,7 +251,7 @@ ndFit <- function(W,
                              uniform_CI = uniform_CI,
                              verbose = verbose)
 
-    results_adjust <- list(ABCD = data.frame(taxon = 1:J,
+    results_adjust <- list(Psi2 = data.frame(category = 1:J,
                                              est = psi_ABCD$res$est,
                                              se = psi_ABCD$res$se,
                                              lower_marg = psi_ABCD$res$lower_marg,
@@ -236,7 +260,7 @@ ndFit <- function(W,
                                              upper_sim = psi_ABCD$res$upper_sim,
                                              type = psi_ABCD$res$type),
 
-                           ABCD_g = data.frame(taxon = 1:J,
+                           Psi2g = data.frame(category = 1:J,
                                                est = psi_ABCD$res_g$est,
                                                se = psi_ABCD$res_g$se,
                                                lower_marg = psi_ABCD$res_g$lower_marg,
@@ -252,14 +276,17 @@ ndFit <- function(W,
   }
 
   # now return results
-  results <- list(noadjust = results_simple,
-                  adjust = results_adjust,
-                  nuis = nuis,
-                  cf_nuis = cf_nuis,
-                  variance = list(noadjust = list("Sigmahat" = Sigmahat,
-                                                  "Sigmahat_g" = Sigmahat_g),
-                                  adjust = list("Sigmahat" = psi_ABCD$Sigmahat$Sigmahat_log_AIPW,
-                                                "Sigmahat_g" = psi_ABCD$Sigmahat$Sigmahat_g_log_AIPW)))
-
-  return(results)
+  if (adjust_covariates) {
+    results <- list(coef = results_adjust,
+                    nuisances = nuis,
+                    crossfit_nuisances = cf_nuis,
+                    variance = list("Sigmahat" = psi_ABCD$Sigmahat$Sigmahat_log_AIPW,
+                                    "Sigmahat_g" = psi_ABCD$Sigmahat$Sigmahat_g_log_AIPW))
+  } else {
+    results <- list(coef = results_simple,
+                    variance = list("Sigmahat" = Sigmahat, "Sigmahat_g" = Sigmahat_g))
+  }
+  results$call <- call
+
+  return(structure(results, class = "ndFit"))
 }
diff --git a/R/print.ndFit.R b/R/print.ndFit.R
new file mode 100644
index 0000000..6a4e2ee
--- /dev/null
+++ b/R/print.ndFit.R
@@ -0,0 +1,44 @@
+#' Print function
+#'
+#' @param x Object of class \code{ndFit}
+#' @param n The number of coefficient estimates to be printed (ordered by largest absolute value to smallest)
+#' @param ... No optional arguments are accepted at this time.
+#'
+#' @return \code{NULL}. Displays printed model summary.
+#'
+#' @method print ndFit
+#'
+#' @export
+print.ndFit <- function(x,
+                        n = 20,
+                        ...) {
+
+  cat("\nCall:\n",
+      paste(deparse(x$call), sep = "\n", collapse = "\n"), "\n\n", sep = "")
+
+  cat("\nCoefficient estimates with the largest magnitudes (Psi1/Psi2):\n")
+  n_or_nrow <- min(n, nrow(x$coef[[1]]))
+  sorted_ind <- order(abs(x$coef[[1]]$est), decreasing = TRUE)[1:n_or_nrow]
+  coefs_tab <- x$coef[[1]][sorted_ind,, drop = FALSE]
+  cols_NA <- which(colMeans(is.na(x$coef[[1]])) == 1)
+  if (length(cols_NA) > 0) {
+    coefs_tab <- coefs_tab[, -(cols_NA)]
+  }
+  print(data.frame(coefs_tab))
+
+  cat("\n")
+
+  cat("\nCoefficient estimates with the largest magnitudes (Psi1g/Psi2g):\n")
+  n_or_nrow <- min(n, nrow(x$coef[[2]]))
+  sorted_ind <- order(abs(x$coef[[2]]$est), decreasing = TRUE)[1:n_or_nrow]
+  coefs_tab <- x$coef[[2]][sorted_ind,, drop = FALSE]
+  cols_NA <- which(colMeans(is.na(x$coef[[2]])) == 1)
+  if (length(cols_NA) > 0) {
+    coefs_tab <- coefs_tab[, -(cols_NA)]
+  }
+  print(data.frame(coefs_tab))
+
+  cat("\n")
+
+  message("To obtain the entire coefficient table, use the command `ndFit_object[[1]]$coef` for the Psi1 or Psi2 parameter and `ndFit_object[[2]]$coef` for the Psi1g or Psi2g parameter.")
+}
diff --git a/README.Rmd b/README.Rmd
index cd30df1..63d6194 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -16,42 +16,53 @@ knitr::opts_chunk$set(
 # niceday
 
 <!-- badges: start -->
+[![R-CMD-check](https://github.com/statdivlab/niceday/workflows/R-CMD-check/badge.svg)](https://github.com/statdivlab/niceday/actions)
 <!-- badges: end -->
 
-`niceday` is an `R` package for estimating fold-differences in multivariate outcomes across binary treatments and in the presence of other covariates. It uses a nonparametric approach, which means that under a given set of identifiability assumptions, no assumptions about the distribution of outcomes or the relationship between covariates and outcomes need to be made.
+`niceday` is an `R` package for non-parametrically estimating fold-differences in multivariate outcomes. `niceday` specifically estimates log-fold differences in covariate-weighted conditional means of an outcome *even when the desired outcome may not be directly observed*. Instead, the observed outcome reflects the desired outcome but with category- and sample-specific multiplicative distortions. One of the key advantages of niceday is that no assumptions about the distribution of outcomes or the structural relationship between the outcome and adjustment covariates (e.g., linearity) need to be made.
+
+We expect `niceday` to be especially useful for **microbiome researchers**, because 
+
+1. High-throughput sequencing of microbiomes displays variation in sequencing depth, and [unequal detection of taxa](https://elifesciences.org/articles/46923) (e.g., due to extraction bias, amongst other things). Therefore, `niceday` allows estimation of fold-differences on the true abundance scale, but just the observed sequencing scale. 
+2. Researchers often want to adjust for additional covariates (either known confounders or simply to make more reasonable comparisons between groups), but don't know what form the adjustment should take -- e.g., linear vs nonlinear, etc. 
+
+We welcome your feedback and questions, and hope you find this package useful!
 
 ## Installation
 
-To install and load `niceday`, use the following code.
+To install and load `niceday`, use the following code. You may get some messages about the loaded packages, but these aren't a problem. 
 
-``` r
+```{r, message=FALSE}
 # install.packages("remotes")
-remotes::install_github("statdivlab/niceday")
+# remotes::install_github("statdivlab/niceday")
 library(niceday)
 ```
 
 ## Use
 
-Consider a dataset with `20` samples split into a treatment and control group and `30` outcomes. `W` is the multivariate outcome matrix, `A` is the treatment vector, and `X` is the matrix of covariates that we would like to adjust for. 
+Here's an example of how to run `niceday`'s main fitting function, `ndFit`. Check out the vignettes for more information! Here, we have an observed data matrix `W`, metadata `data`, contrast of interest `A`, and adjustment covariates `X`.
 
 ```{r example}
 library(niceday)
-
-W <- data.frame(matrix(rpois(20 * 30, 100), nrow = 20, ncol = 30))
-A <- rep(0:1, each = 10)
-X <- data.frame(rnorm(20, 0, 1), sample(c(1:3), size = 20, replace = TRUE))
-
-# fit niceday model
-niceday_res <- ndFit(W = W, A = A, X = X)
-# look at estimates across categories
-plot(niceday_res$adjust$ABCD_g$taxon, niceday_res$adjust$ABCD_g$est)
+data(EcoZUR_meta)
+data(EcoZUR_count)
+my_ndfit <- ndFit(W = EcoZUR_count[, 1:50], # consider only the first 50 taxa to run quickly
+      data = EcoZUR_meta,
+      A = ~ Diarrhea,
+      X = ~ sex + age_months,
+      num_crossval_folds = 2, # use more folds in practice
+      num_crossfit_folds = 2, # for cross validation and cross fitting
+      sl.lib.pi = c("SL.mean"), # choosing single learner for the example to run quickly,
+      sl.lib.m = c("SL.mean"))  # in practice would use other options as well
 ```
 
+If you want to go deep down the rabbit hole, you can look at how we ran it for the data analysis in our paper [here](https://github.com/statdivlab/niceday_supplementary/tree/main/data_application/R_scripts). 
+
 ## Citation
 
 If you use `niceday` for your analysis, please cite our manuscript. 
 
-Grant Hopkins, Sarah Teichman, Ellen Graham, and Amy Willis. "Nonparametric Identification and Estimation of Ratios of Multi-Category Means under Preferential Sampling." 
+Grant Hopkins, Sarah Teichman, Ellen Graham, and Amy Willis. "Nonparametric Identification and Estimation of Ratios of Multi-Category Means under Preferential Sampling." https://arxiv.org/abs/2510.23920 
 
 ## Bug reports and feature requests 
 
diff --git a/README.md b/README.md
index b508d3d..01e0625 100644
--- a/README.md
+++ b/README.md
@@ -5,56 +5,71 @@
 
 <!-- badges: start -->
 
+[![R-CMD-check](https://github.com/statdivlab/niceday/workflows/R-CMD-check/badge.svg)](https://github.com/statdivlab/niceday/actions)
 <!-- badges: end -->
 
-`niceday` is an `R` package for estimating fold-differences in
-multivariate outcomes across binary treatments and in the presence of
-other covariates. It uses a nonparametric approach, which means that
-under a given set of identifiability assumptions, no assumptions about
-the distribution of outcomes or the relationship between covariates and
-outcomes need to be made.
+`niceday` is an `R` package for non-parametrically estimating
+fold-differences in multivariate outcomes. `niceday` specifically
+estimates log-fold differences in covariate-weighted conditional means
+of an outcome *even when the desired outcome may not be directly
+observed*. Instead, the observed outcome reflects the desired outcome
+but with category- and sample-specific multiplicative distortions. One
+of the key advantages of niceday is that no assumptions about the
+distribution of outcomes or the structural relationship between the
+outcome and adjustment covariates (e.g., linearity) need to be made.
+
+We expect `niceday` to be especially useful for **microbiome
+researchers**, because
+
+1.  High-throughput sequencing of microbiomes displays variation in
+    sequencing depth, and [unequal detection of
+    taxa](https://elifesciences.org/articles/46923) (e.g., due to
+    extraction bias, amongst other things). Therefore, `niceday` allows
+    estimation of fold-differences on the true abundance scale, but just
+    the observed sequencing scale.
+2.  Researchers often want to adjust for additional covariates (either
+    known confounders or simply to make more reasonable comparisons
+    between groups), but don’t know what form the adjustment should take
+    – e.g., linear vs nonlinear, etc.
+
+We welcome your feedback and questions, and hope you find this package
+useful!
 
 ## Installation
 
-To install and load `niceday`, use the following code.
+To install and load `niceday`, use the following code. You may get some
+messages about the loaded packages, but these aren’t a problem.
 
 ``` r
 # install.packages("remotes")
-remotes::install_github("statdivlab/niceday")
+# remotes::install_github("statdivlab/niceday")
 library(niceday)
 ```
 
 ## Use
 
-Consider a dataset with `20` samples split into a treatment and control
-group and `30` outcomes. `W` is the multivariate outcome matrix, `A` is
-the treatment vector, and `X` is the matrix of covariates that we would
-like to adjust for.
+Here’s an example of how to run `niceday`’s main fitting function,
+`ndFit`. Check out the vignettes for more information! Here, we have an
+observed data matrix `W`, metadata `data`, contrast of interest `A`, and
+adjustment covariates `X`.
 
 ``` r
 library(niceday)
-
-W <- data.frame(matrix(rpois(20 * 30, 100), nrow = 20, ncol = 30))
-A <- rep(0:1, each = 10)
-X <- data.frame(rnorm(20, 0, 1), sample(c(1:3), size = 20, replace = TRUE))
-
-# fit niceday model
-niceday_res <- ndFit(W = W, A = A, X = X)
-#> Loading required package: SuperLearner
-#> Loading required package: nnls
-#> Loading required package: gam
-#> Loading required package: splines
-#> Loading required package: foreach
-#> Loaded gam 1.22-5
-#> Super Learner
-#> Version: 2.0-29
-#> Package created on 2024-02-06
-#> Loading required package: xgboost
-# look at estimates across categories
-plot(niceday_res$adjust$ABCD_g$taxon, niceday_res$adjust$ABCD_g$est)
+data(EcoZUR_meta)
+data(EcoZUR_count)
+my_ndfit <- ndFit(W = EcoZUR_count[, 1:50], # consider only the first 50 taxa to run quickly
+      data = EcoZUR_meta,
+      A = ~ Diarrhea,
+      X = ~ sex + age_months,
+      num_crossval_folds = 2, # use more folds in practice
+      num_crossfit_folds = 2, # for cross validation and cross fitting
+      sl.lib.pi = c("SL.mean"), # choosing single learner for the example to run quickly,
+      sl.lib.m = c("SL.mean"))  # in practice would use other options as well
 ```
 
-<img src="man/figures/README-example-1.png" width="100%" />
+If you want to go deep down the rabbit hole, you can look at how we ran
+it for the data analysis in our paper
+[here](https://github.com/statdivlab/niceday_supplementary/tree/main/data_application/R_scripts).
 
 ## Citation
 
@@ -62,7 +77,7 @@ If you use `niceday` for your analysis, please cite our manuscript.
 
 Grant Hopkins, Sarah Teichman, Ellen Graham, and Amy Willis.
 “Nonparametric Identification and Estimation of Ratios of Multi-Category
-Means under Preferential Sampling.”
+Means under Preferential Sampling.” <https://arxiv.org/abs/2510.23920>
 
 ## Bug reports and feature requests
 
diff --git a/data/EcoZUR_count.rda b/data/EcoZUR_count.rda
new file mode 100644
index 0000000..81e217f
Binary files /dev/null and b/data/EcoZUR_count.rda differ
diff --git a/data/EcoZUR_meta.rda b/data/EcoZUR_meta.rda
new file mode 100644
index 0000000..950deec
Binary files /dev/null and b/data/EcoZUR_meta.rda differ
diff --git a/man/EcoZUR_count.Rd b/man/EcoZUR_count.Rd
new file mode 100644
index 0000000..5b24742
--- /dev/null
+++ b/man/EcoZUR_count.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/EcoZUR_count.R
+\docType{data}
+\name{EcoZUR_count}
+\alias{EcoZUR_count}
+\title{Count data from EcoZUR study, example data used in examples and vignettes}
+\format{
+A data frame with with 358 rows and 229 columns. Each row is an observation and each column is a microbial taxon.
+}
+\source{
+Jesser et al. (2023). "Why are so many enteric pathogen infections asymptomatic? Pathogen and gut microbiome characteristics associated with diarrhea symptoms and carriage of diarrheagenic E. coli in northern Ecuador." Gut Microbes 15(2).
+}
+\usage{
+EcoZUR_count
+}
+\description{
+A dataset containing the count matrix from the EcoZur study.
+}
+\keyword{datasets}
diff --git a/man/EcoZUR_meta.Rd b/man/EcoZUR_meta.Rd
new file mode 100644
index 0000000..7d26a87
--- /dev/null
+++ b/man/EcoZUR_meta.Rd
@@ -0,0 +1,26 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/EcoZUR_meta.R
+\docType{data}
+\name{EcoZUR_meta}
+\alias{EcoZUR_meta}
+\title{Metadata from EcoZUR study, example data used in examples and vignettes}
+\format{
+A data frame with 358 rows and 5 variables:
+\describe{
+  \item{Diarrhea}{Indicator of whether sample is from a participant with diarrhea (Case) or a participant without diarrhea (control) (character)}
+  \item{sex}{Sex of participant (character)}
+  \item{age_months}{Age in months of participant (integer)}
+  \item{urban}{How urban of a population the sample is from, from 0 to 5 (integer)}
+  \item{site}{Study site (character)}
+}
+}
+\source{
+Jesser et al. (2023). "Why are so many enteric pathogen infections asymptomatic? Pathogen and gut microbiome characteristics associated with diarrhea symptoms and carriage of diarrheagenic E. coli in northern Ecuador." Gut Microbes 15(2).
+}
+\usage{
+EcoZUR_meta
+}
+\description{
+A dataset containing metadata from the EcoZur study.
+}
+\keyword{datasets}
diff --git a/man/est_nuis.Rd b/man/est_nuis.Rd
index 728b934..fa8bb47 100644
--- a/man/est_nuis.Rd
+++ b/man/est_nuis.Rd
@@ -31,13 +31,13 @@ est_nuis(
 
 \item{num_crossval_folds}{Number of folds for cross-validation. Default is `10`.}
 
-\item{gtrunc}{?? Default is `min(0.05, 5 / sqrt(NROW(W)) / log(NROW(W)))`.}
+\item{gtrunc}{Truncation parameter, bounding the estimated propensity scores away from `0` and `1` by `gtrunc`. Default is `min(0.05, 5 / sqrt(NROW(W)) / log(NROW(W)))`.}
 
-\item{sl.lib.pi}{Libraries used to estimate ???. Default is `c("SL.mean", "SL.glm.binom", "SL.glmnet.binom", "SL.gam.binom", "SL.ranger.binom", "SL.hal9001.binom")`.}
+\item{sl.lib.pi}{Libraries used to estimate the propensity score nuisance function. Default is `c("SL.mean", "SL.glm.binom", "SL.glmnet.binom", "SL.gam.binom", "SL.ranger.binom", "SL.hal9001.binom")`.}
 
-\item{sl.lib.m}{Libraries used to estimate ???. Default is `c("SL.mean", "SL.glm.pois", "SL.glmnet.pois", "SL.ranger.pois")`.}
+\item{sl.lib.m}{Libraries used to estimate conditional mean mu_j nuisance functions. Default is `c("SL.mean", "SL.glm.pois", "SL.glmnet.pois", "SL.ranger.pois")`.}
 
-\item{sl.lib.q}{Libraries used to estimate ???. Default is the input to `sl.lib.pi`.}
+\item{sl.lib.q}{Libraries used to estimate conditional probabilities on nonzero observations. Default is the input to `sl.lib.pi`.}
 
 \item{allow_warnings}{Allow warnings? Default is `TRUE`.}
 
@@ -46,13 +46,20 @@ est_nuis(
 \item{verbose}{Do you want to receive updates as this function runs? Default is `TRUE`.}
 }
 \value{
-A list containing elements `noadjust`, `adjust`, `nuis`, `cf_nuis`, and `variance`. `noadjust` gives unadjusted parameter estimates.
-`adjust` gives covariate adjusted parameter estimates. `nuis` gives estimates of nuisance parameters. Add a description of the rest here!
+A list of estimated nuisance functions and additional diagnostic information.
 }
 \description{
 Estimate nuisance parameters for niceday model
 }
 \examples{
-# add example here!
+data(EcoZUR_meta)
+data(EcoZUR_count)
+nuis <- est_nuis(W = EcoZUR_count[, 1:50], # consider only the first 50 taxa to run quickly
+                 A = (EcoZUR_meta[, "Diarrhea"] == "Case") * 1,
+                 X = model.matrix(~ sex + age_months, EcoZUR_meta),
+                 num_crossval_folds = 2,
+                 num_crossfit_folds = 2,
+                 sl.lib.pi = c("SL.mean"), # choosing single learner for the example to run quickly,
+                 sl.lib.m = c("SL.mean"))  # in practice would use other options as well
 
 }
diff --git a/man/ndFit.Rd b/man/ndFit.Rd
index 029b7dd..131d33b 100644
--- a/man/ndFit.Rd
+++ b/man/ndFit.Rd
@@ -7,7 +7,8 @@
 ndFit(
   W,
   A,
-  X,
+  X = ~1,
+  data,
   alpha = 0.05,
   uniform_CI = TRUE,
   d = 0.1,
@@ -19,7 +20,6 @@ ndFit(
   sl.lib.pi = c("SL.mean", "SL.lm", "SL.glm.binom", "SL.xgboost.binom"),
   sl.lib.m = c("SL.mean", "SL.lm", "SL.glm.qpois", "SL.xgboost.pois"),
   sl.lib.q = sl.lib.pi,
-  adjust_covariates = TRUE,
   nuis = NULL,
   verbose = FALSE,
   cross_fit = TRUE
@@ -28,17 +28,19 @@ ndFit(
 \arguments{
 \item{W}{Matrix of multivariate outcome variables.}
 
-\item{A}{Binary covariate of interest.}
+\item{A}{Formula indicating binary covariate of interest.}
 
-\item{X}{Matrix of covariates to adjust for.}
+\item{X}{Formula indicating covariates to adjust for. Default is no adjustment, `~ 1`}
+
+\item{data}{Data frame containing any covariates in formulas for `A` and `X`.}
 
 \item{alpha}{Desired asymptotic type I error rate. Default is `0.05`.}
 
-\item{uniform_CI}{?? Default is `TRUE`.}
+\item{uniform_CI}{Generate simultaneous confidence intervals across al categories? Default is `TRUE`.}
 
 \item{d}{Parameter for smoothed median centering. Default is `0.1`.}
 
-\item{gtrunc}{?? Default is `0.01`.}
+\item{gtrunc}{Truncation parameter passed to `est_nuis`, bounding the estimated propensity scores away from `0` and `1` by `gtrunc`. Default is `0.01`.}
 
 \item{bs_rep}{Number of replicates for bootstrap sampling of uniform confidence intervals. Default is `1e5`.}
 
@@ -48,28 +50,37 @@ ndFit(
 
 \item{enforce_pos_reg}{Should estimates of \eqn{E[W_j|A=a,X]} to forced to be strictly positive? Default is `FALSE`.}
 
-\item{sl.lib.pi}{Libraries used to estimate ???. Default is `c("SL.mean", "SL.lm", "SL.glm.binom", "SL.xgboost.binom")`.}
-
-\item{sl.lib.m}{Libraries used to estimate ???. Default is `c("SL.mean", "SL.lm", "SL.glm.qpois", "SL.xgboost.pois")`.}
+\item{sl.lib.pi}{Libraries used to estimate the propensity score nuisance function. Default is `c("SL.mean", "SL.lm", "SL.glm.binom", "SL.xgboost.binom")`.}
 
-\item{sl.lib.q}{Libraries used to estimate ???. Default is the input to `sl.lib.pi`.}
+\item{sl.lib.m}{Libraries used to estimate conditional mean mu_j nuisance functions. Default is `c("SL.mean", "SL.lm", "SL.glm.qpois", "SL.xgboost.pois")`.}
 
-\item{adjust_covariates}{Should covariates be adjusted for in the estimator? Default is `TRUE`.}
+\item{sl.lib.q}{Libraries used to estimate conditional probabilities on nonzero observations. Default is the input to `sl.lib.pi`.}
 
-\item{nuis}{??. Default is `NULL`.}
+\item{nuis}{Optional list of estimated nuisance functions (from function `est_nuis()`). Default is `NULL`.}
 
 \item{verbose}{Do you want to receive updates as this function runs? Default is `FALSE`.}
 
 \item{cross_fit}{Should cross-fitting be run? Default is `TRUE`.}
 }
 \value{
-A list containing elements `noadjust`, `adjust`, `nuis`, `cf_nuis`, and `variance`. `noadjust` gives unadjusted parameter estimates.
-`adjust` gives covariate adjusted parameter estimates. `nuis` gives estimates of nuisance parameters. Add a description of the rest here!
+A list containing elements `coef`, `nuisances`, `crossfit_nuisances`, `variance`, and `call`. `coef` contains estimates,
+standard errors, and confidence intervals for all categories. `nuisances` and `crossfit_nuisances` contain estimated nuisance
+parameters and nuisance parameters from the crossfitting procedure (included when adjustment covariates are included).
+`variance` contains the estimated covariance matrix. `call` contains the call to `ndFit()`.
 }
 \description{
 Fit niceday model
 }
 \examples{
-# add example here!
+data(EcoZUR_meta)
+data(EcoZUR_count)
+ndFit(W = EcoZUR_count[, 1:50], # consider only the first 50 taxa to run quickly
+      data = EcoZUR_meta,
+      A = ~ Diarrhea,
+      X = ~ sex + age_months,
+      num_crossval_folds = 2, # use more folds in practice
+      num_crossfit_folds = 2, # for cross validation and cross fitting
+      sl.lib.pi = c("SL.mean"), # choosing single learner for the example to run quickly,
+      sl.lib.m = c("SL.mean"))  # in practice would use other options as well
 
 }
diff --git a/man/print.ndFit.Rd b/man/print.ndFit.Rd
new file mode 100644
index 0000000..40e3f95
--- /dev/null
+++ b/man/print.ndFit.Rd
@@ -0,0 +1,21 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/print.ndFit.R
+\name{print.ndFit}
+\alias{print.ndFit}
+\title{Print function}
+\usage{
+\method{print}{ndFit}(x, n = 20, ...)
+}
+\arguments{
+\item{x}{Object of class \code{ndFit}}
+
+\item{n}{The number of coefficient estimates to be printed (ordered by largest absolute value to smallest)}
+
+\item{...}{No optional arguments are accepted at this time.}
+}
+\value{
+\code{NULL}. Displays printed model summary.
+}
+\description{
+Print function
+}
diff --git a/tests/testthat/test-ndFit.R b/tests/testthat/test-ndFit.R
index bedb9fb..de9fd0f 100644
--- a/tests/testthat/test-ndFit.R
+++ b/tests/testthat/test-ndFit.R
@@ -1,14 +1,17 @@
 test_that("ndFit works as expected", {
 
-  # fix this to test something more meaningful!
-  W <- data.frame(matrix(rpois(20 * 30, 100), nrow = 20, ncol = 30))
-  A <- rep(0:1, each = 10)
-  X <- data.frame(rnorm(20, 0, 1), sample(c(1:3), size = 20, replace = TRUE))
+  data(EcoZUR_meta)
+  data(EcoZUR_count)
+  ndRes <- ndFit(W = EcoZUR_count[, 1:50], # consider only the first 50 taxa to run quickly
+                 data = EcoZUR_meta,
+                 A = ~ Diarrhea,
+                 X = ~ sex + age_months,
+                 num_crossval_folds = 2,
+                 num_crossfit_folds = 2,
+                 sl.lib.pi = c("SL.mean"),
+                 sl.lib.m = c("SL.mean"))
 
-  # fit niceday model
-  niceday_res <- ndFit(W = W, A = A, X = X)
-
-  expect_type(niceday_res, "list")
+  expect_true(inherits(ndRes, "ndFit"))
 
   # also add tests that things fail when expected (if you give the wrong inputs), and that test different
   # argument values
diff --git a/vignettes/intro_niceday.Rmd b/vignettes/intro_niceday.Rmd
new file mode 100644
index 0000000..d54afd8
--- /dev/null
+++ b/vignettes/intro_niceday.Rmd
@@ -0,0 +1,155 @@
+---
+title: "Introduction to niceday"
+output: rmarkdown::html_vignette
+author: "Grant Hopkins, Sarah Teichman, and Amy Willis"
+date: "`r Sys.Date()`"
+vignette: >
+  %\VignetteIndexEntry{intro_niceday}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>"
+)
+```
+
+First we will install `niceday` if haven't already.
+
+```{r, eval = FALSE}
+# if (!require("remotes", quietly = TRUE))
+#     install.packages("remotes")
+# 
+# remotes::install_github("statdivlab/niceday")
+```
+
+Next we will load `niceday`, as well as `tidyverse` suite packages that we will use.
+
+```{r setup, message = FALSE}
+library(niceday)
+library(dplyr)
+library(ggplot2)
+```
+
+## Introduction 
+
+This vignette provides an introduction to `niceday`, a statistical package that implements a nonparametric approach to estimating ratios of multi-category means under preferential sampling. First, let's break this down:
+
+- "nonparametric": we identify the parameters that we target without imposing parametric assumptions on the data distribution or the relationship between the outcome and covariates
+- "ratios of multi-category means": this refers to the parameter that we target. Specifically, we estimate ratios of means (also referred to as fold-differences), i.e. the multiplicative differences in means of the same category across levels of a binary treatment
+- "preferential sampling": we assume that observed data are proportional to "absolute-scale" data that we wish we could observe, with potential unknown sample-specific and category-specific scaling factors 
+
+In `niceday`, we identify our parameters using assumptions about our observed data and parameters. We assume the following:
+
+1. Our expected observed level of category $j$ in sample $i$ is proportional to the true level of category $j$ in sample $i$, an unknown category-specific effect, and an unknown sample-specific effect. This is trivially fulfilled if there are no category- and/or sample-specific effects. 
+
+2. Conditional on the treatment that we care about and any additional measured covariates, the true category levels, sample-specific effects, and category-specific effects are independent of each other. 
+
+3. Category-specific effects are independent of the treatment that we care about and any additional measured covariates, for all categories. 
+
+All parameters in `niceday` rely on assumptions 1-3. There are two types of parameters targeted in `niceday`, and they differ due to which of the following assumptions is preferred: 
+
+4. Sample-specific effects are independent of the treatment that we care about and any additional measured covariates. 
+
+This is the stronger independence assumption. With this assumption we can estimate log fold-differences in true category levels, across the binary treatment (with or without accounting for additional measured covariates). We refer to these parameters as `Psi1` (no additional covariates) and `Psi2` (including additional covariates).
+
+5. Conditional on additional measured covariates, sample-specific effects are independent of the treatment. 
+
+This is the weaker independence assumption. With this assumption we can estimate log fold-differences in true category levels across binary treatment levels relative to the median fold-differences across all categories. We refer to these parameters as `Psi1g` (no additional covariates) and `Psi2g` (including additional covariates).
+
+## Data
+
+We will demonstrate use of `niceday` using a dataset from the EcoZUR study, [published](https://pmc.ncbi.nlm.nih.gov/articles/PMC10730187/) by Jesser et al. (2023). We will use metadata from the study as well as a count table, which provides observed levels of a set of microbial taxa from each sample. 
+
+```{r}
+data("EcoZUR_meta")
+data("EcoZUR_count")
+```
+
+Let's take a quick look at this data! 
+
+```{r}
+EcoZUR_meta %>% head()
+EcoZUR_count[1:5, 1:5]
+```
+
+## Fitting model
+
+Now that we have the data, we can use the `niceday` fit function `ndFit()` to estimate parameters and generate confidence intervals. 
+
+```{r}
+ndRes <- ndFit(W = EcoZUR_count, 
+               data = EcoZUR_meta,
+               A = ~ Diarrhea,
+               X = ~ sex + age_months,
+               num_crossval_folds = 2, # use more folds in practice
+               num_crossfit_folds = 2, # for cross validation and cross fitting
+               sl.lib.pi = c("SL.mean", "SL.glm.binom"), # choosing single learner for the example to run quickly,
+               sl.lib.m = c("SL.mean", "SL.glm.pois", "SL.gam.pois"), # in practice would use other options as well
+               verbose = TRUE)  
+```
+
+Let's look closer at some of the arguments for `ndFit()`. 
+
+- `W`: this is our observed data table, with observations on the rows and categories on the columns. Here's its the number of times the ASV was observed in each sample. 
+- `data`: this is our metadata as a data frame (which includes all covariates in `A` and `X` formulas)
+- `A`: this is the binary covariate that we would like to estimate fold-differences across (given as a `formula`)
+- `X`: these are the covariates that we would like to adjust for in our model
+- `num_crossval_folds` and `num_crossfit_folds`: the numbers of folds used for cross-validation and cross-fitting. We recommend using the defaults of `10` for each of these, we just use `2` in this example to run more quickly.
+- `sl.lib.pi` and `sl.lib.m`: the types of statistical learners that we use for our nuisance parameters (these are additional parameters that we must estimate in order to estimate our target parameters). In this example we only use the simplest learner "SL.mean", but in practice we recommend a wider set of learners to account for more complex relationships between the outcome and covariates. For our analysis, we used `sl.lib.pi = c("SL.mean", "SL.glm.binom", "SL.glmnet.binom", "SL.gam.binom", "SL.ranger.binom", "SL.xgboost.binom", "SL.hal9001.binom")` and `sl.lib.m = c("SL.mean", "SL.glm.pois", "SL.glmnet.pois", "SL.gam.pois", "SL.ranger.pois", "SL.xgboost.pois")`. It takes longer, but gives us the performance guarantees that we'd like. 
+
+There are more optional arguments for `ndFit()` that we won't cover here. Feel free to post an issue if you have questions about them! 
+
+## Investigating results 
+
+Now, let's look at our results. The results for our target parameters are in the `coef` component in our `ndFit` object.
+
+```{r}
+names(ndRes$coef)
+```
+
+We can see that the `ndFit$coef` object includes two components, `Psi2`, which includes parameter estimates and confidence intervals for log fold-difference parameters (relying on stronger independence assumption 4), and `Psi2g`, which includes parameter estimates and confidence intervals for log fold-difference relative to the median log fold-difference across categories parameters (relying on weaker independence assumption 5). 
+
+```{r}
+ndRes$coef$Psi2 %>% head()
+ndRes$coef$Psi2g %>% head()
+
+(ndRes$coef$Psi2$est - ndRes$coef$Psi2g$est) %>% head()
+```
+
+Here we can see that the set of parameters `Psi2` and `Psi2g` are the same up to a linear shift - this is the effect of centering all values by the same factor. 
+
+How can we interpret these results? To get the effect size on the fold scale (rather than the log-fold scale), exponentiate the estimate, and take the inverse ($1/\hat{\Psi_j}$) to flip the group comparison. Here's an example:
+
+- `Psi2` = `r signif(ndRes$coef$Psi2$est[1], 2)`: We estimate that the average abundance of ASV 1 is `r signif({psihat <- ndRes$coef$Psi2$est[1]; ifelse(psihat >0, exp(psihat), 1/exp(psihat))}, 3)` times greater in `r {psihat <- ndRes$coef$Psi2$est[1]; ifelse(psihat >0, "Diarrheal", "Non-Diarrheal")}` cohort.
+- `Psi2` = `r signif(ndRes$coef$Psi2$est[2], 2)`: We estimate that the average abundance of ASV 2 is `r signif({psihat <- ndRes$coef$Psi2$est[2]; ifelse(psihat >0, exp(psihat), 1/exp(psihat))}, 3)` times greater in `r {psihat <- ndRes$coef$Psi2$est[2]; ifelse(psihat >0, "Diarrheal", "Non-Diarrheal")}` cohort.
+- `Psi2g` = `r signif(ndRes$coef$Psi2g$est[1], 2)`: We estimate that the average abundance of ASV 1 is `r signif({psihat <- ndRes$coef$Psi2g$est[1]; ifelse(psihat >0, exp(psihat), 1/exp(psihat))}, 3)` times greater in `r {psihat <- ndRes$coef$Psi2g$est[1]; ifelse(psihat >0, "Diarrheal", "Non-Diarrheal")}` cohort, compared to typical fold-differences across groups.
+- `Psi2g` = `r signif(ndRes$coef$Psi2g$est[2], 2)`: We estimate that the average abundance of ASV 2 is `r signif({psihat <- ndRes$coef$Psi2g$est[2]; ifelse(psihat >0, exp(psihat), 1/exp(psihat))}, 3)` times greater in `r {psihat <- ndRes$coef$Psi2g$est[2]; ifelse(psihat >0, "Diarrheal", "Non-Diarrheal")}` cohort, compared to typical fold-differences across groups.
+
+You should state which assumptions you're using in the methods section of your work! Remember -- to use `Psi2`, you need to make a stronger assumption (independence of sampling effects and contrast of interest). 
+
+Finally, let's plot our estimates and 95% confidence intervals. We will move forward with our weaker assumption 5, with the parameter `Psi2g`. 
+
+```{r}
+ndRes$coef$Psi2g %>%
+  mutate(includes_zero = ifelse(lower_marg < 0 & upper_marg > 0, TRUE, FALSE)) %>% 
+  ggplot(aes(x = category, y = est, color = includes_zero,
+             ymin = upper_marg, ymax = lower_marg)) + 
+  geom_point() + 
+  geom_errorbar() + 
+  labs(x = "Category", y = "Log fold-difference estimate",
+       color = "Interval includes zero") + 
+  theme_bw() + 
+  theme(legend.position = "bottom") + 
+  ylim(c(-5, 5))
+```
+
+Here we can see that our 95% confidence intervals range from approximately $-4.5$ to $4$, and approximately 75% of the categories in our analysis have 95% confidence intervals that include $0$. 
+
+## Citation
+
+If you use `niceday` please cite the paper:
+
+Hopkins et al. (2025+) "Nonparametric Identification and Estimation of Ratios of Multi-Category Means under Preferential Sampling." https://arxiv.org/abs/2510.23920 
diff --git a/vignettes/niceday_vignette.Rmd b/vignettes/niceday_vignette.Rmd
deleted file mode 100644
index f571f18..0000000
--- a/vignettes/niceday_vignette.Rmd
+++ /dev/null
@@ -1,36 +0,0 @@
----
-title: "Introduction to niceday"
-author: "Amy Willis"
-date: "`r Sys.Date()`"
-output: rmarkdown::html_vignette
-vignette: >
-  %\VignetteIndexEntry{Introduction to niceday}
-  %\VignetteEngine{knitr::rmarkdown}
-  %\VignetteEncoding{UTF-8}
----
-
-```{r, include = FALSE}
-knitr::opts_chunk$set(
-  collapse = TRUE,
-  comment = "#>"
-)
-```
-
-First, install `niceday` if you haven't already.
-
-```{r, eval = FALSE}
-# if (!require("remotes", quietly = TRUE))
-#     install.packages("remotes")
-# 
-# remotes::install_github("statdivlab/niceday")
-```
-
-Now we can load `niceday`.
-
-```{r}
-library(niceday)
-```
-
-## Introduction 
-
-[use this vignette to demonstrate use of `niceday` on data analysis from the paper!]