niceday version 1.0.0

.Rbuildignore

@@ -3,3 +3,5 @@
 ^LICENSE\.md$
 ^README\.Rmd$
 ^\.github$
+^doc$
+^Meta$

.gitignore

@@ -3,3 +3,5 @@
 .RData
 .Ruserdata
 inst/doc
+/doc/
+/Meta/

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

NAMESPACE

@@ -1,4 +1,5 @@
 # Generated by roxygen2: do not edit by hand
 
+S3method(print,ndFit)
 export(est_nuis)
 export(ndFit)

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"

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"

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,

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,15 +204,15 @@ 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,
                                           upper_marg = upper_log_noadj_marg,
                                           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"))
 }

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.")
+}