From 8eb43ed24f985fd4aaf709061ea200248faffb13 Mon Sep 17 00:00:00 2001 From: Sarah Teichman Date: Mon, 6 Oct 2025 15:34:26 -0700 Subject: [PATCH 1/5] updating inputs --- DESCRIPTION | 2 +- R/est_nuis.R | 6 +++--- R/ndFit.R | 49 ++++++++++++++++++++++++++++++------------------- 3 files changed, 34 insertions(+), 23 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 756ed99..23944c4 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")) diff --git a/R/est_nuis.R b/R/est_nuis.R index abb6dff..ca05c11 100644 --- a/R/est_nuis.R +++ b/R/est_nuis.R @@ -6,9 +6,9 @@ #' @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 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`. diff --git a/R/ndFit.R b/R/ndFit.R index baa57e2..6432e48 100644 --- a/R/ndFit.R +++ b/R/ndFit.R @@ -1,21 +1,21 @@ #' 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 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`. #' @@ -28,7 +28,8 @@ #' @export ndFit <- function(W, A, - X, + X = ~ 1, + data, alpha = 0.05, uniform_CI = TRUE, d = 0.1, @@ -46,7 +47,6 @@ 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) { @@ -56,12 +56,22 @@ ndFit <- function(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") } - if (length(A) != nrow(W)) { - stop("The length of A does not match the number of rows of W") + # transform A from formula into a vector + A <- 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 (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 +83,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 <- model.matrix(X, data) + if (ncol(X) == 1) { + adjust_covariates = TRUE + } else { + adjust_covariates = FALSE } if (!is.null(nuis)) { From 8db64c7476064d9d7249f0201bca40eea6093d00 Mon Sep 17 00:00:00 2001 From: Sarah Teichman Date: Tue, 7 Oct 2025 11:48:02 -0700 Subject: [PATCH 2/5] adding print function and doing more documentation --- NAMESPACE | 1 + R/est_nuis.R | 5 ++--- R/ndFit.R | 24 +++++++++++++----------- R/print.ndFit.R | 32 ++++++++++++++++++++++++++++++++ man/est_nuis.Rd | 9 ++++----- man/ndFit.Rd | 22 +++++++++++----------- man/print.ndFit.Rd | 21 +++++++++++++++++++++ 7 files changed, 84 insertions(+), 30 deletions(-) create mode 100644 R/print.ndFit.R create mode 100644 man/print.ndFit.Rd 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/est_nuis.R b/R/est_nuis.R index ca05c11..b044781 100644 --- a/R/est_nuis.R +++ b/R/est_nuis.R @@ -5,7 +5,7 @@ #' @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 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`. @@ -13,8 +13,7 @@ #' @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 some diagnostic information. #' #' @examples #' # add example here! diff --git a/R/ndFit.R b/R/ndFit.R index 6432e48..557f221 100644 --- a/R/ndFit.R +++ b/R/ndFit.R @@ -7,7 +7,7 @@ #' @param alpha Desired asymptotic type I error rate. Default is `0.05`. #' @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`. @@ -263,14 +263,16 @@ 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)) + } + + return(structure(results, class = "ndFit")) } diff --git a/R/print.ndFit.R b/R/print.ndFit.R new file mode 100644 index 0000000..12912cf --- /dev/null +++ b/R/print.ndFit.R @@ -0,0 +1,32 @@ +#' 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:\n") + n_or_nrow <- min(n, nrow(x$coef)) + sorted_ind <- order(abs(x$coef$est), decreasing = TRUE)[1:n_or_nrow] + coefs_tab <- x$coef[sorted_ind,, drop = FALSE] + cols_NA <- which(colMeans(is.na(x$coef)) == 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$coef`.") +} diff --git a/man/est_nuis.Rd b/man/est_nuis.Rd index 728b934..7eb24e7 100644 --- a/man/est_nuis.Rd +++ b/man/est_nuis.Rd @@ -33,11 +33,11 @@ est_nuis( \item{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,8 +46,7 @@ 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 some diagnostic information. } \description{ Estimate nuisance parameters for niceday model diff --git a/man/ndFit.Rd b/man/ndFit.Rd index 029b7dd..ed9e5ba 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,13 +28,15 @@ 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`.} @@ -48,15 +50,13 @@ 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`.} 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 +} From eac882451212c5f18ca45e5530ab1e773a09fdc4 Mon Sep 17 00:00:00 2001 From: Sarah Teichman Date: Tue, 7 Oct 2025 14:54:08 -0700 Subject: [PATCH 3/5] building vignette and readme --- .Rbuildignore | 2 + .gitignore | 2 + DESCRIPTION | 19 +++--- R/EcoZUR_count.R | 7 +++ R/EcoZUR_meta.R | 14 +++++ R/est_nuis.R | 12 +++- R/ndFit.R | 36 +++++++---- R/print.ndFit.R | 10 +-- README.Rmd | 21 ++++--- README.md | 98 +++++++++++++++++++++++------ data/EcoZUR_count.rda | Bin 0 -> 47861 bytes data/EcoZUR_meta.rda | Bin 0 -> 2033 bytes man/EcoZUR_count.Rd | 19 ++++++ man/EcoZUR_meta.Rd | 26 ++++++++ man/est_nuis.Rd | 14 ++++- man/ndFit.Rd | 19 ++++-- tests/testthat/test-ndFit.R | 19 +++--- vignettes/intro_niceday.Rmd | 120 ++++++++++++++++++++++++++++++++++++ 18 files changed, 369 insertions(+), 69 deletions(-) create mode 100644 R/EcoZUR_count.R create mode 100644 R/EcoZUR_meta.R create mode 100644 data/EcoZUR_count.rda create mode 100644 data/EcoZUR_meta.rda create mode 100644 man/EcoZUR_count.Rd create mode 100644 man/EcoZUR_meta.Rd create mode 100644 vignettes/intro_niceday.Rmd 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 23944c4..05172c7 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -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/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 b044781..d7e24dc 100644 --- a/R/est_nuis.R +++ b/R/est_nuis.R @@ -13,10 +13,18 @@ #' @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 of estimated nuisance functions and some diagnostic information. +#' @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 557f221..9d14ca1 100644 --- a/R/ndFit.R +++ b/R/ndFit.R @@ -19,11 +19,22 @@ #' @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, @@ -51,6 +62,8 @@ ndFit <- function(W, 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") @@ -61,7 +74,7 @@ ndFit <- function(W, } # transform A from formula into a vector - A <- model.matrix(A, data) + A <- stats::model.matrix(A, data) if (ncol(A) > 2) { stop("The `A` formula must only include a single covariate") } @@ -84,11 +97,11 @@ ndFit <- function(W, } # transform X formula into design matrix - X <- model.matrix(X, data) + X <- stats::model.matrix(X, data) if (ncol(X) == 1) { - adjust_covariates = TRUE - } else { adjust_covariates = FALSE + } else { + adjust_covariates = TRUE } if (!is.null(nuis)) { @@ -191,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, @@ -199,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, @@ -238,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, @@ -247,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, @@ -273,6 +286,7 @@ ndFit <- function(W, 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 index 12912cf..c946c95 100644 --- a/R/print.ndFit.R +++ b/R/print.ndFit.R @@ -17,10 +17,10 @@ print.ndFit <- function(x, paste(deparse(x$call), sep = "\n", collapse = "\n"), "\n\n", sep = "") cat("\nCoefficient estimates with the largest magnitudes:\n") - n_or_nrow <- min(n, nrow(x$coef)) - sorted_ind <- order(abs(x$coef$est), decreasing = TRUE)[1:n_or_nrow] - coefs_tab <- x$coef[sorted_ind,, drop = FALSE] - cols_NA <- which(colMeans(is.na(x$coef)) == 1) + 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)] } @@ -28,5 +28,5 @@ print.ndFit <- function(x, cat("\n") - message("To obtain the entire coefficient table, use the command `ndFit_object$coef`.") + 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..a072c47 100644 --- a/README.Rmd +++ b/README.Rmd @@ -16,6 +16,7 @@ knitr::opts_chunk$set( # niceday +[![R-CMD-check](https://github.com/statdivlab/niceday/workflows/R-CMD-check/badge.svg)](https://github.com/statdivlab/niceday/actions) `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. @@ -32,19 +33,21 @@ 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. +Consider a count matrix `W`, metadata `data`, the treatment `A`, and covariates to adjust for `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) +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 ``` ## Citation diff --git a/README.md b/README.md index b508d3d..6d5b07c 100644 --- a/README.md +++ b/README.md @@ -5,6 +5,7 @@ +[![R-CMD-check](https://github.com/statdivlab/niceday/workflows/R-CMD-check/badge.svg)](https://github.com/statdivlab/niceday/actions) `niceday` is an `R` package for estimating fold-differences in @@ -26,35 +27,94 @@ 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. +Consider a count matrix `W`, metadata `data`, the treatment `A`, and +covariates to adjust for `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 +#> Loaded gam 1.22-6 +#> Loading required package: MASS +#> Loading required package: glmnet +#> Loading required package: Matrix +#> Loaded glmnet 4.1-10 +#> Loading required package: hal9001 +#> Loading required package: Rcpp +#> hal9001 v0.4.6: The Scalable Highly Adaptive Lasso +#> note: fit_hal defaults have changed. See ?fit_hal for details +#> Loading required package: ranger +#> Loading required package: xgboost +#> Loading required package: withr +#> Loading required package: SuperLearner +#> Loading required package: nnls #> 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) +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 +#> +#> Call: +#> ndFit(W = EcoZUR_count[, 1:50], A = ~Diarrhea, X = ~sex + age_months, +#> data = EcoZUR_meta, num_crossval_folds = 2, num_crossfit_folds = 2, +#> sl.lib.pi = c("SL.mean"), sl.lib.m = c("SL.mean")) +#> +#> +#> Coefficient estimates with the largest magnitudes: +#> category est se lower_marg upper_marg lower_sim upper_sim +#> 13 13 -3.2741166 0.6118426 -4.47330614 -2.07492704 -5.2715008 -1.2767324 +#> 3 3 -2.6157685 0.7393277 -4.06482416 -1.16671283 -5.0293326 -0.2022044 +#> 22 22 2.1413861 0.9688809 0.24241444 4.04035780 -1.0215634 5.3043356 +#> 39 39 2.0849190 0.7705275 0.57471289 3.59512520 -0.4304980 4.6003361 +#> 30 30 -1.8174642 0.5879262 -2.96977845 -0.66515004 -3.7367723 0.1018438 +#> 14 14 1.7374620 0.4585679 0.83868546 2.63623844 0.2404493 3.2344746 +#> 21 21 1.6905260 0.5946739 0.52498671 2.85606538 -0.2508099 3.6318620 +#> 20 20 1.5157064 0.5558763 0.42620878 2.60520394 -0.2989736 3.3303863 +#> 33 33 -1.4345277 0.5776259 -2.56665373 -0.30240164 -3.3202101 0.4511547 +#> 32 32 -1.3580239 0.6794349 -2.68969173 -0.02635606 -3.5760655 0.8600177 +#> 15 15 -1.3004635 0.6041093 -2.48449595 -0.11643105 -3.2726018 0.6716748 +#> 31 31 1.2779538 0.6321374 0.03898736 2.51692033 -0.7856833 3.3415909 +#> 45 45 -1.2336750 0.4405577 -2.09715223 -0.37019774 -2.6718928 0.2045428 +#> 27 27 1.1777250 0.6579759 -0.11188403 2.46733410 -0.9702630 3.3257130 +#> 42 42 -1.0418668 0.5415186 -2.10322383 0.01949026 -2.8096755 0.7259420 +#> 12 12 1.0119062 0.8021781 -0.56033395 2.58414645 -1.6068355 3.6306480 +#> 25 25 1.0117772 0.3666362 0.29318338 1.73037094 -0.1851211 2.2086754 +#> 23 23 1.0073781 0.5866781 -0.14248979 2.15724602 -0.9078554 2.9226116 +#> 8 8 -0.8832539 0.4902861 -1.84419704 0.07768917 -2.4838120 0.7173042 +#> 4 4 -0.8691163 0.3683387 -1.59104691 -0.14718560 -2.0715724 0.3333399 +#> type +#> 13 TMLE +#> 3 TMLE +#> 22 TMLE +#> 39 TMLE +#> 30 TMLE +#> 14 TMLE +#> 21 TMLE +#> 20 TMLE +#> 33 TMLE +#> 32 TMLE +#> 15 TMLE +#> 31 TMLE +#> 45 TMLE +#> 27 TMLE +#> 42 TMLE +#> 12 TMLE +#> 25 TMLE +#> 23 TMLE +#> 8 TMLE +#> 4 TMLE +#> 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. +``` ## Citation diff --git a/data/EcoZUR_count.rda b/data/EcoZUR_count.rda new file mode 100644 index 0000000000000000000000000000000000000000..81e217f1747621a6f02d03188cc6fa0b9f56e90e GIT binary patch literal 47861 zcmV)5K*_&CT4*^jL0KkKStdOCZ~;Qt|NsC0|NsC0|NsC0|NsC0|NsC0|NsC0|NsC0 z|NsC0|Nr2wz7ItqRZ&SnK}6Q7rP%93yQb{s%#45o00000pcEt#AOwX91AysNCE1g6 zHpJ<`0007#NjFqgcJ%4bAOH#@?|MZSy#N$XaC@oX1mW4hPU^=8yypjZz2(zQ9RLB% zXz+aXK*ynTzVA>_z7)X+3yG}9n@BWe#& z&^e6(S*pU{t+~J3F>W328y1>ri0ZUqiUPdpQ)ImWT%O> z40=@4$kFO$9;d0Mo})~ksig5H)ICpAPek=S6Vwk$k*BCnzWHa+lx{>IdWP?4ZuyM| z3h6LIBj%!jf&pZN8ApkxL=Z{!bxa$GFp4smWD*GhFTh|TP1h6BRB{c(n~69=41aK- ze}&0hL!AtRg?N8X z2Qrx+6O@ch>y!(Vkz-z7!(B}%ecB9+ArH!RYu}sg7gye#G^E_>hdlC1(A^-re8*5a zSvps8dbE7k7&67tDsjwrDG$MpZVhpN(_MP3QV#;V69yW)OCh1cP$z=CP@s~ElHR;_ zP10Y#KvFx$`KO+B5lj%T6}QBl4&MXdVH=42la_MA{yb?2k!(L1;IcnHi3)OjtZns` z)VcP4ru71ZZZqswCIsBShxm$IfUgiLG7Nm17Ud2YTo(S7CgdFGtYH?XT)evCKhG^u_-?+koXY3Ci)!bw5=m_qb8@^+SUJIAO$i=G_I z2^}BMV)-Sg<3NxWGsE>xk4p%DiL1ZAUb*E zB!aJ6;WJGrs(QO{N=yy`;ji9*qSYY)eJl8+5aIURzaxuV*;;LU_u+>thbsT(jpr9ZW?mY?5C@a7)M z^-fv~qy!{Mx@%y>_)Ad(eLFjB_1e5A5J=4~I=jfjDmfz!k)#?P$7U_fLbgzN71=o; z{+G&Q^3Up-rE|22BSIIRPzV$v4#s((CMm(geb9|BBhBgaBv;pW`tqA{Rk7(305R>I z-2)-#9bcAvdiLPH@GSesllnZ8$}eNO->VdbFQ^Ns@q3!&$w?S4DsJ2D>$kwm0|E%V zbr4-5ET7=-k%88#+^DodZ8z>(`vRK=a9p+1DeXFwmhAb<6 zI%s4Xxp2|_=;lpX;+X>Ogb??kC>=Q)`ptmVM)=YLH^j{C8r2LJp_`@Dfsgj9^gKAF z<5?e_HLqGfoTml*+Ml>&Wv5A#h2ipdBJ_{rNtSsz8Sm!3$=19T3QwF9DUVh=K9%V9 z(v)Iie?-eB_Ca)QlFWJ`rs|>aSWpg(Uh0{6`S92DBjJp#r(fEA_wwMf9V2>e;btx% zw*vMcU#W7|l_Dr_YmfQP0i>rv+lZYIrV3gCep^gm{B`>fZ9)bMlwGI+-GANc_&vNE z&r9omK*{|-0ut-r6k_-5RAK-%W@xQ(-VE*#2(lh^ROr1We0Wr2;Fy9F|8gWPKQXEb ze_Ns|`#lNn4%W!>jhNyhzAe!~t@^Q#LsCA7iyb9y!{DUH zs$D0i5Z>_5+*|inz-~4DZymA9St6NozEzC1viEV%$ZAeQ9**sLY|#mSSK#|h!t*}< zS+2wYjfG(UEYEWc6+nhQschHmgNq0KKOcEp4A8M#R}^ak!GV#`z_)1E^xR+AZbeC? zm-QyX+q;#yMZDa`;coZ%#^X6^|9Nwc6z?_vUu;t!v0f(oQn+!~pxY98@y_0T+t`2S zTag&@H|?}T(HJ9%1)x25hOIEKpS)}%sDCn5y>ap1r6bqJPO;ehH^Hbsk-9Y4Pny8Y zXd)s6awM5VIzmDULmjeHOsW;ICgDq}BV=EUV=4BKyPll~O1{S+pyRZKDf*0BQ23Odtm$iF~A1v~pdfHorlzRO>&7^rLjR)4$VN{-c zEYE^+s-`dFz`&0$XVJQ0ka=-zW4E@|fm2hu^@eBLI7uK?&j;eZw<6*A6Ix8x&A}Nx z|A9#WurlT|{pgA>+Rd`nr2}%tDZAlt(%ln%slgO~mACJm(bp)Dyajc<{GOnOEj*w4 z|F=w^bSZ6L$u-`r3;Cm(m^S{$&8mssL{O%F2}E8&|4#<4BO+8F-O+L#FlG1DFQ^LNHM)n%9D{SK8R9N#yB?kT5F!ouk{3@RU*W0?b^veH_oMqWk=wtlw`|B9oSaGc+W zau9Wy6`n+E7_4yt@AHrDtKCzL58)n22e~?wy!hnO z854jpPjCN(fr0^HmBLyQhwK<*dVv$-c~v*H_Kpi1O|H!k`t;2+J#U5iYBMeLp%uVn zZF1B3<-y?Or@UZye$Lp77vPNk&%PKE?0=7pm9L_I-&}f1MA9D!!$WUgGHWb%pu!PK z9Z`>W*^DNBek`7QAby~;M9{9roLCJ}Xo599`;XL4lG|jl>aRR)N> zf&xSfX^*}gNP96j2-e?1Qq7S|t->1_TlF>crBff7M&Ui-gODIaH}-81C+pCeu&J zh|`IvDKiDr)Y17I-Th>zfSM2>0fYF0$Dm3w({fGT38xYq7t=D<)u(8k=Fk%4iN({ zd}!$TP0>rERk*^+Zm(9ec`CS=t1?ey=TfJ^h@Zqr6A=}oUSm?FqzC##XlbcD(Qs)$ z`bHZZsqOns+5Xb&O5KMXkO*~h3>zaJ#E_0e$yWOxXwXqp0Km2tV~r5ZBH)6n)}9o3 zkK3l7(LQ{Q_cqJZvJgYYR)x>xz`rV9@f1_C*v;Gj-@EctxxPUU%opQ5f;d9qEV3UU( zB9452LzFU`3uE^Z^C_du_h_@&$z851$8jZJcKKTpUtaccI`eW9;4hXX;&YgOj8#eY z*&I1*PbbFTm2YtEDjw+m`E5U8NsN2m?*->{5}(0bgz1w`oCjv5yA$iGa|}{HBQ*?O zCyvt+A}Y7PSHCX#ae>4J5j~j4>bmn!cRU!V`&fS}DXFnkC_QgRiEKJ-*lQuV)Ws@PMQ^9w!IxmxC8iZZBw| z#K>R$+Jh^k9*N{nT!G{_SddY171s?ODEni0+t$?dl!Tid>Dh)aY>~xcWI~=t-@K`X zK>yxA-ZlhE97y9adV_dXVy;Acw-^zAhvgBpLIYWiY{uOe?{mLGy$?yu#z&jC?bZiH zIP!F4@pAq|{hu~N*ECn}Kc%8k$?>WcFOvBcG-7y+sQ|RK#u8mB#Q41#2*}mv+bMrF z-j?4xCeD429vZ+buOxp@=4(esmM6SeyL~@W{I3+@XZ*hTrfTZ*<;YX|=Vn*x*@Pc| z@%m)7)bmN%`J>nFG~xRH^3cH@4{N!XzF*2?d`VqSN5x|7t^!hkUnfQ1_$Ydq)$x6A z!zVF%qWHVrE7E?G8drW$zRyw@Ih-j>#JwMbNq?akxmI`o62T7tncAD1-8I(jz1z`( zf5`Po_BtPCs4x0H>&FPmL+*^PC%J(d0E|>YcnAH+*6%r-a&jf(|1aqtf7U=eJxnUU zaYmjeAbo=Lh5rkH6ON7AA)Jv#_~HUV46fb_nrF-i`e{rb&Tskd8|wUKT^MG@pjK@# zJk*Hdvr~Zg$zl+54zoBqbRt>2jb3Frul@&=3jgBy8f9)&;N#x7;8>@dMA^v?)5xZz zrv~>DC1>h=c_Abukmx-=lp-R>QD-{X;aYA_P;F8{?A`= z3`EYoqq4^~nc}-QIupC~d=v7(oW2%OWD4tnSJVF^G3hq~ZY273{H-=3{I#uiNP4b% z2srL(h-p;wqhG*3_Bk<}S~9XJvxJbktuXECE(JU=Z_p|lu1TRcSz;Y?LZG|&H4A(O z9uu}j{JFC{hmul$+Q>@g;z;1tYT|?QTIA3SEX}yxxFHaQ5QO%=Jwg#TD?&Q$|7tF1 zp?V$Jq2a=Ut2b>8<`X3CUQlP~^c?2;wQ42P0y#ihAC=QS2vKU0hhutGo*`uh6No$Y z68DFJ6pn(5E_wpsVnNw!(r)AJ_LF?{5PB?J(9LV%&j0pZrfFsxA$Ng+x;1e;g%cNe=gaEC*+g=pE(Nu%D5M_-bfT--L}A?4nzXXPp6$E4cS$eY^b zoDrd=c|nTnizAit$)ynOlB70yZk%F;nR2_>9%KP*5FhlALD8I4Uzx_8@z9jb%h|sR z%?yO1@8WdfsK(y({8vke2BJOJPFXTw|5e)j>&IE_*gZuX!5H~5{5+Y>ka)Hy+!o!F z_Z)_uaq1p)xl?pF=nN!YVd)Qz)IBY4t(=-&$s2<-Kik-+tiF+x%I@7wlZ|RHaNuoH zc58lwK?xCqwsIbQg6&A>LmaA-3zJiAn>={YvXoO%obqX@5pLQ(h-8cb>s)iQJHT|h z;OmIi2cXEn&I;~`=;l^h-q>mSG}NbgE_{_l_p7c?FLE=un?pVI=B@fSK=MHlix#91 zJEh%9tJ~wthO<6={m!T4>wNS?7%4qlBf-0%J`lU2b2XS-1X#LV4!!59ZsLrI?i9|6 z{{!N|iP4*yxjSCxd|rMh>aN-wO09kJueeaJiO1g zt`0j}=7Qc!f3+`UUApq_?1loW314%%jIGK(Xh?Lf98{ zQyd!eVf!{+C_KI%oZ0*QWiAzRT<|kh8{0^>9-vxIC@lh>2@(@r!3XN=hV%_9KQ0IpQgp= z&<3lOfSfQEaZZ`H>w7_0r4$Po=_bk` zJXIzKP29IGd!q2>pM+aQeV=z`y`SW!#locs_lPS(Pff zVbozutJ<>>Sd+W~ork*LM0%1IJ>eGkwyt-fb@r>~@b_Q6yP< z`S2ejJ42UU8)ES2I(eq~6blO#4RNI+D#;(_7 z%}a*^-mS*IxoF5%x@>4ojH}5K9CvW}UH6Gr=Q@@4b9u+A4J{jMGhOhz@`P-?2qZr--Q3OCNn>v1>*bfNOWy_%kEhd^RybO3m$fF1GAfRF zI^6PapN}DFDr3HCCI6d!8U1xIEX}{4MSOC3l0$E%Pd)FOtCcL}vo;uvuqfPLAxNwH4iO`vl4oo|H# z3_b~$f1335_?q_-3x*GG=BH1awpu)>A*&nOsPHd<;f^Xtk@;wY-x$ripIFbzUON*pPER>BJ65<>xYa zQPJ4)p7ql%lZ(mSs<|j`^mqEOdJb<(%NCy9y&gH}vqwe=0dpPw7R^^e#o|pp>dNZz zZ*ZNR^s@NN7mL~6jaz8d2X;k|mnyv6hck(iG7Tlmu}Reki%3{JM;0pJPVyd>awlUu3JR1woG4n(a%Rz=dnW!H8q?z8%0{IcL75^gv&*3_PALuWV%4XiCTv(_Wh~KetS|XR#_}9PL_Mtt4Ywv#XsJ z9#BS;vxX&LJt&f=XwlGsoo*Ifi#lh#(+J1?)9l*5%-#DE=E~NX5-C8buZJ|Z_#05T z88=|+TA$Rka^t%2`p!)bTXl<*y`;WESL+>i=WnJb-4oolLD!2rUDDXppKSwopo;Dj zqpz-2X6;*(vfDTrEA=qU;5sCg&LisOBuUGB-M&(Z8cvQxk>y4qDaD}y_h4_MJD;D( zEeaFP%%(|SDm!M*w?H6Z$l~p(bsj}kE90}WRLUrIBv;`^y#1`{=l@cK$argDKCcAix*RHI3*deQ}6F5_YP!{Er{%rl@br2+@@AGZqgO?Ftu<1V<){KVk?-%1BI-EvZ3g{x zQHfzHj|c4zB1^IXyepGD{Xm;CdHsj!abE;IEa@07t#39OJGtErh;K~0 z0R)?@(U({%9lZ<6Mz+N{a8luHTf7M03RB=$ZUWFtuQQ_==J&?3yO0)H$F&{TwJIls zcLsL5t)s_ZmB{^`&s(NYrT*uu+hu%_xzW=9UF_uWY{=AIuWGeu+sUD`Uzfz5s#%jR zXHM7EyT^$(y|tjhjXq^G!!|z5?SAQup)Fb=A>y%+lnBqKXB3@244#Qfob0vYEu7HF zYkYfaTvKr@d`t4NWc-JtRLsPRE*Jkuk?F9}z1w zIODghTWaz@7R~y(SMwIEv}sNov8;2#$=sQtZw%G45sXueRdpkg^yBDI`PrOxb6}?$ zPsF`btIE5rnD=)JYG$oI0i0Sj57N@Q>Hr^bKmZ^RARz#NNFaxbEZ+9VBz!YPliBin z{^c!D06~Eu0Dxkl8dmwp+E;Zns1xATmie zBF{XBW_4nH5dVMAjel#)fn+(^*<>GM~vn)N;{qqBy8 zY4&^lSCvdML*aMS#(kExU4=547+X^#cXGh}U%JxS@lo!}>i$ky*1fH;_mHL#S27Vz zSnD@IywlrW?vO_SIpTn+z@)IM1i>EY_}!Hq#_xlahHer;A12Q2@P|<J4ztH)^l%6O34s&{^nLex8}RuG5Bn&IrnC)ld8IVXJowOb3bpY|CV5cUMNXDg zG?x$wS5$~0{iG)&uI^;VS)=3XSa(?lF9dR&;(6=DV82Tb>FlNc@W5UGFWFA(-&^5@ zi3h51J^RmPdpeLH8tv~-6knC9M~MSv2bgP(UsLau%#3Nxe*ZKFOI9x&)}A(O&-`3& zJ#+dTtk|89c>HptqJck{WryYJ=w5;RO~2pyaBPmT$(y{z63Il*o=b|M3MB@Sp@w@q z;)K@BG7S)oki#H_f5yH<80|z~AzXE4h8`(b24vvBMaue-(E9T8p|Ipxdvrj2y5Rv^IC&tm~REpf8C$(8hRO|Y1_j;0l@Pb5eyoH&d{t7 zT5+x0hDB;HC?6n|&Tu2l)j**8X&47t#T%J0{oQ|gKP#>`)XIe8h2^^FhH0Q3{R0W#(ESPd(hv-~5%! zvyZXAb~$vPV$|kn5!=YFY;g8#h}0@5Dy5*LY9=B4TBbcsFZgz+WwJVfiK;ZIwZ)Z(5Ty(gj>-M6JzDXc1fYn%@gAM*$D zL;mW7j%!j>1YESi*@_Ey&Ns=lrQ1G*B}bS6F*WfqKLPM3Qj-l_jyAA3U3o3^G4r$k zXT)SiC@q)>k59Z;3?@+lih{4=DWwVqK2CzPgxxH3dzPHqmiMjeqf~XyWGQ`iM!)&{Bq(;9_0u&f2Y=@$E8c!cMv9y6 zP6jAH0D>4#&w$|~%)uL*10Wt4);~wn8SxUmCsjPzp5;mz%A}&j z9?@U%RW{_f5lX^$MG57{3ePcsMKx8_7ydDQhkxfkYtxT54WGkjvzRfL=B);;#y~=% z7by8cpZMxF&9G1^Dd;zcRbhRx-4whQ} zab2K3wD&z7f3jNsqbvW3_!x8ZV?k0Lm%5jHqv)g*dzQ$Q@s3c@`NeR9_UjJ6uFhI9 z)tf>SKgWq5V_(U3BZsc*PZLM7ez)FQ{eO-7^_2WOE+R-e+6L;_SfOjNt|qgUxU`lODg$2KyO zCEh6Uu9&99mj`gB5A)fM{DAOIe?dRVgM>nX!J2)R4mJ3q8_KEqa{g27{p>^6W)bk0 zpXHx)h;Un*a#rvYD?tFrp_5~j1fYNRFHJ~tOQG2#(EFPNmBXOX$n+X0*lU2_j?};% z!GIaMCywbQ>aeC^@liBxN;IwbUoQ@K#t3(5`;+!KxA5xTt3HUvS6z=_WM9vks`QC; z)xVENT>!^T`!1Qh=(QT*sr0`Y^t;orNz?Jx=_nla$Tz>3aBvATi=uP_LQFp8LwRL1 zA|@{gMNJSKg!WZ-sDBA-=6xfxHz#0XjlzUJawx>7L7kbuuYN!#IA7 zE+gwTfweoK))APnms;j`FpnuX`itw#?0wq;NX`iL(CxjK6!E)=XxU3#I&0ezck4_S zSr5(&A)hm|JgrYnH}BsOlKY|0PO0>wu8e#08-1MwBNsO|eu#5`@`sEwa+fxTznML^ zr60~#ZRo4DwbIx@x56iPi*1%n!w9gyYI^=$eocqmVGjL$wkLtj&N&fnzv*)2$$@!@fb8Vx*&3VI!e&cu-{@!0 zv%uV=Ys8vW^KT|@n`^)MYx8TfCTRNhe0Y2D+0EvhU(>gSY|_&qIL01**T;r_xG}ph zCGBZu7g&_%li2 zyL_BOcsw(xVOkVze}dlJqdy=M>e zy%IpJ{RfXMHQeupPR^NQBKdB!LAyFf(`#q!dniDKj*+IF)2qNky-35KT^U^R{8M%F zs;IPOm!1VDrp}dE1R*^&&zzujjC9kL_q* zEXX;qBS9VgNJ2*Rx)*{B1hsC_n4?Nq&~kNQQELGFdS)+x)(iPOu|nHWef34)(243vAw+q@jPP#QGTr+0Xx;k} zMtouZcdzHq!M|uWHL{`u#}L+Lriw1 zzNKbu=r^X_+?kTzzYp~feT=MBB*uZ+w;T_P!H5*)u~2o1EEA$0zVAJ2y;}+3zReT& zdc_<>>EC2vgjDdGv3&GBd=DGF&#q~hPUfmzny16kufoTS7g|3Q6ysG4?7FPKkL9`f zG?^A*{Vr#-NUr?;e_2V<-Z|>TrIn)~SLi_Rk+zGvN1r182H{n+r4s$CR&^lpPO+Gw zVw6}&j-1_kP)z|_&GL7r)RQk~Jur%ih-pkW;QpB0@)ZE|Y6Cxw6=Cphe!TQPP z?$-+abl1DV&**nq;E$6{eu-BaD$6duqOg8jW)}2!;1!2Ed+F&N^ z#vXa_Kpox6k5l#(1B}5YPa*M)Kwo9&|B&a^QRqF}`Zx9XQPp?#QZyfmxh2_*0WDAZ zhWzw)dLi;R*7JP}K}+2#hPQ488!Taiw4y+vyH@zq6(Y!2T4J0}Dvx@_8O6Mfi>$&WZ*{u}Z{r_94{s`&j43JNg_y%KLp` z@Ypgj=r`~7rdDf+M$Pa6IJHtcQ_**e?m^p(uBkfWayOQ{At`;nBH9wu229$er0XE%0s$SjhMt*@thYysZO- znn~$}eLUU$9yEwDP*_1_F2y0}Y1GMm^?0E{pl00PhfeAyi9vlEz|xt+#@J`lf#p!a zGMWuint8H$DPY9We_Splg=f}T$E?E+9;LZL`|>DQMY5=$=W&#hl_d9 zOPC)0*{%4L7rJL8K|a$rFU(RHAh61Mp1%OLtZdubFX{dYf1MYJIW!`*RSDS1^fHAw z#XdM@4BGrF{uF-4lQE;$Z<_jQ?*7c??iwL{)WwT%z|-~f`ONr!tBP^vOkNz5%~Tr= z_2c33x9(&seYbfG$IIE1=#?RU<7ZjnI=iR%((7u{e`Wc${44$C{5)O;+@mdD-W-gG zC~fK0?DL*c6jZZ!*!njc`R~>AIWX!N$$R*5%(Gy&W>Ux5yq2`^?%9SvYKE^i8FBZv zXGThW8tm9{=*^*%oaktphd4;zREJPu3PQl5SO?eM+BiaQ+-M{je=5g_h4e^J!T-Dr zTNH^dik69)46X}(K7}h$y#U<*)lmizeGb78FMf%Hwf@JRf7Uqvsr&P&V&*NDqMn5L zQu;m?NxnaDwt3>d_AUAuO09+xO#u z3Fw?4k(UeX+kxfw2-v)HClt^bMt7#2`J;g@36}jrCI-G2zYtf939gHqOZ9#te&f)) z8wNqB9{{jD6b}nBi#LUabwu5hi<>@4Ft)Wp5)r(?m?>N zxuST3Eo1v`>Uo|QM728Aipq{f=^hcsvp8rwFiqtNS%a<(BMXY zwXpU@GU1kheXDaH0UXH5by-#Bo@mD82+R>jcN?X2US2r8Z!lI3FbbEXxibf8jdeGY zJ03+fXXpH-z4p2Uq(wI_c1;bOz)YQpMZ^(iZ7$|5p~vnzq&u!C)`n&qmh6_NVr&o( zENP=doTPd*>Gt7HnG^$&#-Eg8?l*k3S+M*imp0cMKws1-`W*Ly?op;jggb>Cz5(~ zD~+dYXW-cXKaM$AVC-tmP#xBIJjc$m1C^qpx#NxNHuKG!41D+`=K5mU{97t#@nG&r zI}T9QIAxM*X5WU)_VfnpcJgp^+u_dEKQUz+WbNY4W?{Ea6MCCsviNQ)jYp1HiCrBt zt{)>8bc(8^&b~P5i|xaA76XmcXoy~SF`YeCgaaW55_TN!{d<7s-kHar9@<%&P-Yh# z=H}e@nqs%Vk6e!P8aQl#SjQ0c_j$)10Bz?6az~mG5y(PnwiN49{TCvGS3X>_{nP#K znE0OikYkP@JomH}igHAbY#TSAw!E8L`FJK>n68ZXQ@fE#$r$lCoK@>pA52vaGeqy@ zbmya?(p!VsnmADhr)m~}o8`J75g75c4jGpwvI(alO)r6Y=WyENn(wstmSy!*&9lNK zGx!p5ht<=f=&RWuaqYQ8cx{@$2cJGVv-ARhNb)Vilejg@V5roG>f^my+FT#Uz7j{{DdI@$KST4jugKOWmPSJJqJ?isTv!6c1=^Evq^Nk>@U~J9BI` zrb@8$u){>rr(lJmx-qsS&yS8}#V#WCo5;_*>)tF@DO&!ANhT9e*uu2&or7ah`sChTWsJ>|-CdeP!^+ZRak?0zPlf6=|5>STs-@W`pLeTDM zq7n)Vp8llP6}F*NPS#N5#ijf|ShW08C=w}&9{|W9F^$cGKq>+`IQ?n78uD;8>0*D% z7Qmt&$D)5%tx2SVX}|GT7P<-qVY!Mg5BudTSD}lMLvHZCetRFIZ}ogC)H)~f zEL&5h&p=E{e092%JFI_0(n$8%HdV9lW+&_WU;Cb6z*FshB9FK>YVL(-|Q+#3KY(9DaE?9HpeM{a@;LdE5AB{iUx{54a#AJUVbU58L=D zb6wt=?s8X2)gaR;`_q;e6*H5Vnn?9YDOS6t16`Ea=cMQNTg$_Rf=*khOE z5$~OE+oaJwJF>h*NAA5`Za^CvAg&ZD1{>$tC)4i!5mZPYL!US47Qc1a6FNSV9={8G zYme)%@54vZ5ScD%qrXfkGm^qshU>+X+`%9t__+B*h6RIEG=4q;%$p>j#Y(ZnZ@Vuk z(^lvu$>~e}1ry+PIOdQIEQ{_SkzTZwpHSnSrgHN=f+;#er|YAjflPzNj7^vh0?AhJxxKKQ&LwNm+X13a{gOK&BsK^ixPC8I(k;+4<%)_`2A7j1wT{WxX+%L<1*>-eu+$ko*wYo0IvnPPfGtKh`eW0lGyUK8{OLQMOe52kH*C!f8@0(6I-VcqG1P<^6xmx zHRA17!L=ISBOW}+lpy5QuC5o=4P^Yd>{+!gZgVxz2q*x6LO@U*K$I~K%THR4v$uCK z@;LH5eJ&+k8~S(r`&K0$T?FJ$KX*5{qu&{78MzV-3nKG(xw4?ATB@=q=#|Ps?x{z% zl=fcYd%`AA9@AnUuj5}vbI$=-j3MN1L^}e5H!n%>9%P&4exU=+(O`2uDxg8TlqJ1$ z@rvGzY7Lqd&n=XVC|SSEchm6R(!srMFC)vvsp-qDHgL}t2@{HjfPY@Dzp#dQg-@(r z8FkJ1;uH;@`1<*sitUxw*2UKPV`st1Qs5{2gDgZ=nMeZqtkS`BIT$4}?gV#5Dv9D(>;Md7P-w zwc?1nsFYn)@`R57sHpc9LJ;InsJwMj-5%9dg8{eUYU}`VW z_XMw!Bz(EPipvGSdD7J+XozOW_8Lbk)SB{6<9=?fx4cuqNdR*$C=b9+x_c38tArPR%i!R=fLHC~MBT#k?><$ZP)oz)!w?b9Z{| z8Bd=czDeKNa0@Otx(-9U+lR-aETZ}$-W@WlO^|Qo9D*R%xj%bE;ZrCwc`BEqJ%fE? zh{~zhH+&naGo@Pn*nX5C<_$~V>+ylsdV)qR<~zY>pDON*e56Y*UG=HTbjjuCW;=({ z_zC*_fwh|LxU|q#6&oHT6{5dVln1;)AH{k>K3vYX$*pR9>EW;-4c->0|C6dZ$N&Rs z(BEZ9g>f)odfxx99gBXD%;YSe)QC~=B5v3wYT@KA6 zrS7H?wt9JaJsgH+h$wH??msHF*9B}>6&F)$Wgas7jK#@iLg<)i!Cm%+%3W5VoqYx* zOLf6H6F^zGc)p1>`Vsib51+PHu`()(AA8w7o>@Jqyo~V_K)v_g0@~=12Dwn~i2!g{ zV1z9vKhXGtWW4inTbK#wQFqFS2|GSHqg9{L#z~eX1ZaC!XhrgeS-sSPO%PT$T^+#r zI{l3B23V1=d1YjAdujwj`|eC9{2H2YBC2XYVXDUF`XzbB*I{nOJZF3l*e_A|JPW$Ii%N-;&-e!S-&2s&Y-Ej;$%iVL=4L!?Pu zNI2Yvv}P>`KGOT=d06X)z1$Ypi%N2dSgh>Z{dM22r(U5Z?gu56^Y}WuCT_eTOuA0w zlVkct$(R+*uyABaY~#$qXslX{beWW^F4HyZQc{;wzD+Ed*0_{Yi7x$n^0V*H>oo+- zsv_RD)4cH@AP|HAEkW@4l)%D3022n)euosX@u0`2^CcMawsYdNr}+GdkC(NDx!|SAL8zzL`U+QXtQhyN{Q*H^g zx=w!b5gT$%Ocx#)EbumjL_lDmkrzLxLNsx_D54?>hUwQIwb{&cKn~%YFQlESB^%f6aH>Areee{M5q!-Ak>U7JQ zUhKab_aAEZHrJ*V*2N|9)E8N zIXk)31-MkcBybC1d-_jukiO7gbtB>e{?0zb$jO^jJc?~SYJ4_ik`kxRab|> zK5hJdUp^z#zFkiA?}OcZl*>NDyve7&$J!bJ5Jti;dl#}hP8V}Lm<*&_gb6-+*f1^!%OoO)!56_uD zKf?9!DgG+Ah^_Hz9e#ld9hptO8muhY%9sqJ8^<6*hIdL41b`Fgz*mU85|kEL+FO9bE?NTuy2d_&3Adt#Hr3IE~g_7`uCpQG(^xJGMyDh^Hy(MsdwUCYdCD9nrzkm z;QkJ$;F8om6z}bdJy}%@4CLTHr<*>g8MvIj-}7jG3AEd8@fPEcs{bGS;t?PZABVAO}rMKA*Ksp0gc6UjR2i$iEN_kp7TH{X=_-t*c_TN9)O`Nvj&fg^d zB+N1eykXnYi3mPzgZ3R<+MMp(68n`M2@;nZ;CU3tXk2m>oC zWTJrUzp~ZOO?~l13!pT>WswhY0-I zd&=tQ0aZ-eE!b_hbRL1Osho_tD9?-rGo5n*<5N1K)5%(jbKQ+0Q&rGT3r!L{NOJYy zu?cM@)<7~eLZ-7boST3$vIYulj38w287PEN`7@)O40lL*GC9uf@At9g!16qJV`=gl z85tAmFI&E$Yp(xv$y{m@%HskQ(qTuAFxRH6N)_$Q>yr91Iyn~91aCCrie(TjV|yc^ zU=g^$X2df53(m*h!w7Ds`;kWm7Iv(uHuaq?6tUuQKFTMQ{Iv@lvGg@)bhV@{MMQ?I z=w=AerzpnXH+dPewAosG`c`6~6GVaMT+g^Lx_+P1Dxa|TFUcOr(P9paaBMyDKYN{R z>e5iwMS(}js(3#L)Z7|E2RK9V7l`x}sFszvmZa1d2R6jKtznf&)*HcB8t8-ATopvw z&>4jD#9Oru*6*FjZeI5$bW5(@ZM84#XfW}IiOs8PavP66Oelg?n244_ig4DPjLr=t zFv5&Vpy04?%0diaEzt}wl!O9zWhM-x91V~_;}kICOgh6z!per$A374Y9>57;hccX+ z$ccegmA^zNuN#4U_DjDXp-|Rd*Se9|)qc7Gy~( z1gS7HX681d^S8=Hs9cD!2pg-}K@bpyTvo&b5fTf)IK^+IKAyG$J|OMH87|~viC#t9 z)Xf&gdBiy9r4d?Gb)z1pjTGtfB_AbnOxJNrI&P4-W)`$7E^0h9*`olk`cY-G#NJ1m z$qOn<8g{DG0D}RDW@4&J7Pirdc{d1&M3Dc&K!LC1Fbt|e40t@(49&0j}Z=9u4K;RF=Xyo z-e9EP>TA6Re#4>(L0l#k@zXF<*$}`FOkw=-$Cy!vPT2K?uYf62WUBb6M4o4-->#)O zp|-U>6xpY=p;s3h!KJ?~*ZhJ(7()h6s6i**9iVOQSH6+G$SWiDNh})`(9J>Lq z%M-V`SgTM`MOvAZ>M1;nWU37Yh?~G5>-1Mj)z$Cf*d7Z70Rm}M$-kSv?%$rE0bdne zx%lxgaoCk|Xd?@1)a{z3(Qmtqo>D2*ty$15w|L%;}k+auC*9 zz7d`@SG1w7Aw2vh5f%_9=en~86)Y=Kb$u6|)Z6>N`J~j`n7)mE4a{9lJ#G$nKj$}N zq^ODrbf5qaz9=aI1Q3dtgA7PcHo^8!D&7_*d)lm~wN=MOBUSXyT%~a+f zp$c7&Wers#Ayrl7+@>L=NQqgRWz>{nAyy&aM1@e$CDaZ{jphxG1quWu2_#7@9HJsZ z!_266iA15y!!oQmhlzJ-NFs@qSu#+R7MVsxRg^?GLv%tQB}srONtptIj%*fC}XYRVdsEiSCl5~_(Y48~?~?(E8{u`>}fGczuzi89de%Btvc?+~!F49c$3?(s1Z zAN1%2&)|LX2A2$lyb|18r%9Ydtm$K$h>CRU=K7m(XwLBGnET7IG=U-=#&%X5bsA#S z!<9@8$!n7#J52?fiOq%Va{GO`;;|tunHD`$uEbk(?GmkTo2iA|9;Ta4a#Z0MLBnB` zC5aPwt-ON8;-4c)mrE9lHt4N26jN2YDYZ;pf8z+MzZ3xwg%Cv&hzcP=5D(>qQ3X%> z{tu<*O0a)LDyX0zkHwM$o&K+}_P+1sc)sWB{*Tt={rZXd4a5XPtc4XAzc1}RPy_%U z^}l=NL`DzS_PzX7P!D7Ac^`-U{s;0*{(oco)cGew{`7ue4pj$sezk9*=jr~xz4{#0 zL9gr9mG}NH=C$`M|8BOE2CGUXtShl&VS`Gq<^1^o`tsS|;!vYr5iRMY5A(ni5G_zG z>?T+*j*^3Gt@-t@ zzTVq?^t&kko$@E{LNCe1>?TDI!T(R}m-^fd2ejuU7x;;fm-Kuwd*4&~F~j*J@jjo^ zPp+BK;hW(7FWG0U{$Q#}f@!LweDp0Z{(H_zU9A7J_QijLU&oDKzlhI5w{P~yqc`KZ zbH@+c%6_(`vxRU@LHlKx!)N;E1e)1Y%Y?R{M#}N(J;%r% z=kq>gOy2MIy)Uc#@<>jR$&w-Yh@a2xzG0JSYZsyOTPN;z*jEJl!eNh>!ugCX+*d9Z zCG}Tqci&d4!8*ARx7j@-C$blIFN3pid4ysbs$%i#wHCdoz0@q3p&T5dEJ z6?s;)_;zvWQLk$9WDAz+!LdC8(WlIr3;Qwx<;;nYw_ZtzLK;X~U`o+$3@F+w$MU3# zpd8P$&ge%sXG*H5b#}X73n!;Yr@B0!d*7Z!qr31Z+w;E%^l3Tx!mHNXl8AI;b=LnT zHM$iByeig&Pk92bmk5|hy*~rrOaYlU;rriCRJrT^6ub{YW?bC>@^a%4A<^23JvjT6 zl+z;VE3f5t2$yAH)0~0;2nF@_b^UnNcEyXeR;~(@=g9$;J*2<|_Vz#p7Q0E2EMZ6h zl+p)#0O=q@iwbpzP>E5n$C`ORN#oA$cnS^ zGiMb|tN^mq$pmgdG(p!u3w#a#B=crwb2&VZlBI=w;+)w{D4fL7)?ZIQlfk3p~BDlb2&hxk;hN> z&$chW0A)gxH;~>%P{=*Xqg$1~Ta_-mAgHkhuB#Sx(L>%_Z}j2zJ;YDc{ELW1fh;et zniv@6VP-g|);G~fuEceug7NAj4J1wG$m+2Ttl(4P{gyU{@e5J3pGC8fNq6ACqZZ?em3*NRFg zh}w&_H8h*+yCFnmq1Sy^42?NrKCXT{RJyskcWyQNg>0J?3B@8GlZWVY*oli3OzcS5 zn!uH(lNY3qkDMq2#9 zd7fY{RT`Pp;$`ZpR%yMCBQXJSROJOZ#3)ybT7?<7p@IMb{!ZwR`9c^b2So2LI)d0?9Es* z32Wp$xu?OMF(5J^lsHj^6-E3SyH%!TrLkgr>69yOdu$2JM8aMys{z(hYErf z0Dzx<7oqXJ-JRcA&hlkpRazw^Q?`K+UH)JE30F|IZ$PUDXf6H$(-P`DRbM5)|Q~>W`1PnM( z4p%p6<#V)p&E-Dhy*rNczVCagxixb$=#|Z7zAr+U?GdpqCwOcVNbNB0nu8NdH6V`fa0DvHa(VR~$dNC6tt4AmgM*!b#BapWjhdvx> z8M@N|00)@|$z6+}Ytr!&Ul*Md;Eh4e?58=5+H`yOW~DiZloPgEXg~l0B&jPrB2rQT zX#+|CAcV!=X}FT;Wim%#m0&;>hXF}O1RM#rSQJpD^;r`kVtj0S`t+fd1!zj}V#LJ~ zXTrd-J9o(W2k!nq+Q;4`rd<8nNqyr?%!rEp<@eL?YvSIkBG=gaTgW~+{!37V`Na(v z-Jw6Bmb3j?9s}iBy-q#j#3L8<8e>O6{d@i3yHNO5;eI*ZgE^f~! z&iRF)#lV_wFAx_j9+}f7`#POXdsH)2|8XQZ>&*;g+?@RgJQq3m$U@ex5ydWGBsf z@4~DB=DjS>z}0!%*cMN`NQ@!&`}H97uU1_tYtH_?)%vV^JtF;IH7f)~TOLwwhVnQV zEL_y*7s*6Qd3^LUM_CbHV#;qav}{AYNN0F@r5(whqM7-$6yG$t0+|=s?~$8{AZp~Y zmZ8RSwA!va;HL##*V(>o0oOnP3eOQ_tAsu570HGnMcm!$CT3v`w%S`;CrZ`NDBHhY z3A=qM35_IQNg(blNm3jP)|#@NG2V`Y$}RXV`?D7CJY+uw-#DJnEvm{Y8ke6>D=qK0 zRXNv&?M1}W#0C=tOpaXTW_j=; zy5hR)HYE4zDa}+u;}NlTFsg%jfN+xeGI+polgx$phkJgJG>% zQ(nsO4+I-zGsqVsxyODagQh9sO2I~}lQPR5#WUD4zyaZJmNooH~ltQu$$x2@i7 zb*?nbT9!*WsMaSq)8w^iFNWI?C|#Q%i+Ntu^A#ML5CMq zmorp@IeA#u7~9V)!5ZFY%21BReP0<*PWjJcKNpJ;ltyN?nzXPX#vNSJN_Q8-<{%B;uPF8iShD&6tp<}hXhIti8?>x@853P*|X^RvPdt;nP!OqdcIHr*jQV~(j@|g^9Y&?@cv^q~p2>V(btbV@s>>&d;)VYYq4E4c{+GIE z*YUSC#c>tiIi6%wo4$~!X25YE2IejnG}mx6m{@boF*Q?DG}TI(s&{T$L`+3E$2DkF z=8j%@=Xt*j1y3r5qTe%xGRvnke46Vz&B1`Jr3K!InX6MmceSmwW|)%1tpkRyo?8;5 z-B>P~Y~y-;h*=e2wG<&($z(KjLR!6RiyDbw9Ol!DtgB&)Z#f$g=`!+T2HmInxQZdjRH-bBt^x?j3ZT zt60M42z7k^kf$Q$&|IpXPTCQbcbJ!QTifGTsJ1EM2ouT)Y(=7DTRSdrZUrw!;yXB~`5(A_)BAqF>#l>Yu3)^Il?X3V+~R{0@CGQY(?v@gewVYh(JU0}D0hOO>JLlEpR>?xqUaA`%= z+=^>!$5+^ziHPjVH>`*PxC82lFPLHpE?94QqqqdYW$KfS%%2==S)-f;&i0u|3PWuQ z@ar%_Go7v|*6_CXsnwr151MJMrVdYIo(5Ms6IW_zw|(fA*34A!3>?#Iq?wGkJchd% z$2=YwP}q^g|E?OH3%C8@DvY$YZl z)kT?y*rQp$se1~T%`zrN*LkWd_GZp^vL>VqVP8aAg(_==7FZm@2!`??JdO?oxxFo3 z_Ny-kd33bR$>yv&%nT_uA za4|BeiqQI?! zYNgrqFDZ>-NHD0nVd&2#H;m2Wt=lUihM!T~@wcR;+)%1)_>N(%Cc*f7BhNXV&6QjU zUn`FowdRT}Yapqy&)xkx=l}pd)*u}Kbs6!WJn&6RpD=)q?+{}U?LJy%o4F!G&Q-vK zOO$G!d_?ApDxg>&0vQmo2~|P^7ec5Y!er&Cz@P+KiMNM+2|gUKff~b$X_F@ba%Vej z#HPu1K>~XrtpCtI?^<3?p$$lwkCg+R1`khY zsh}IdUp6z54q-<--GLpL8fSPXbKvKi)=o%q3So{!J6xcVRc}NHtqw30N}1srL=l+D zwZAEb)<_*=YUEkK(qWPta<_POi+yoRZi7+F-6d}J<-RcIU^vd18todOT*j^%OCuY3 zji%YNG}x|PGv{1W;=Zv;)TWyoV9M!M@|@=#4J{H$_CdtVI`vQm8q{JiZJ@@BFcwRO z9tO0U$0)J5U6G=qYIl51t!}=@Tzsl(XuWf$zgMkoZnWS85I|hAtP#!S z6<4gQHIHw+wR=d9B<`f4V(j(jQ)>WPfePTRyY!>*_8r|Hx$bv+UE}dDqHqKFUK_&4 zFW;n-lB%ps_mNz1+L%W!?TCfv`ygFv>*ft*nn4AF&>BmM+ zxl$TBU@^Ag);WgIO){PblBc`Pl=TGr*$5D$JkE)o%BUNs+f@b68-MHt^{wD;Kz55-=Pf#x<;V1*be~B_`Ph!D@@Z zOhb$L!mcyXLUV@PD6};^s>>#^fS$H!^OeNol~ z$Z0Y%Q%KVhK@Y(T#yA|^p7>9U;}2P9Ai_Wyj`|V}^JTGSreI^7Swl-;s*0R8Ra&6Y z9~yuUHxsD}$nu41ZU`I<6}Guqryygn=0!J^Kz-WF$Us%DIc=j%>te3_#ZS@t=sE!E zNmg${siYhoj^;bjw|c4x&pI^~B+#+kqf+Bg9rvCN0d>YJ(kPDbx>cP9C9}rZz+A?z zc+IKE?T?R{t?mxgMlV|7^D+ud_#M>M8V5w*8(%B(B0A~RoM4a3O*DD>J+sXpRo@jVNva%L76=6!)-XA>pTe}g+31Q zWduDRkNCZxPvU>u?%sc1v#atr-_7-XkH^`2I?wLoelO?ypC2de{G$CsdgA+L5Bq(j zMeMvqI{yRUecugJhn}KGUMB$~kurK3`zf&afeE8Df2fmss5kfIp}-ihi67Zd@6`9^ zyx)AU()~MNBooyWfYUw5V2>WDn{`Vw3y>Z-a9gXXb8Nn2@|5GwY%gzV7Y1|ARb^ECx7J zPvA~q{!vazoDf3i=7pxRIVxj_H4g)Wrfi`w-TdxyuPvBUAxE-MQ3 zUg&C+Qd6#4_q|p;*^vpJ(|aw*Uv0snhyNKW-JZTy(Vj71Z`Xs203Zim=ex04$u8+x zn@J%cc1>Y9jV&{?m9hjOIqI&9z#Tq*&3fMNPjn{0z^mlwwmIt0nk+Cr|IXzEu?o#O z-#~)i9gjAd`N48=J{M5J3(yAHn^?AMXfX<}i`*#xwno!?t`V zpU5zG1W-pw5ev1hI*&L1udrpu{*GD-7L~ye^68T{W+G-Qmh!r+U8I{^+pIJy?jp%a zjPhSA>{DoIzh?$<49UoHCup8w*^6~8CaxUBTUpKJb)hYs)>M*#p=A(}0f%KMTn@wVw)O#`|;?iaIC2@4)L)~>oYNCt=kzyqBkmj!8~5GI2npu z(*g#p54jK^m<*kx@EWf{g@3)sb=UCB$PSoBcqt6F;p8K*ya*Dg!}1L^E|pYToHFJI zm;i3wQ(~}-u@cemauOxB^lDX~VmzrfPIjf5lthN0wKgOq7D%b`^HmDDbHbk*q}7(K z=pebG)xOC!TeV~ZDtN|@)uYOB(Kn9AzA+?WFEkn!0R%Hu{f{Iwg5av#uX5e1G+L5P zSraQ+%c7hxwIK261yhonGSn%m++uN;X(zqd$fIru9n=B0Desie_wHmHBs~+6X1lEp_X?LkV9rFg9!hy&>L`7g+chG-5 zYMl{FGRg1Z5#WtBzTppIgg}*~Ata~O+a%mIW~KI%uj1r7I=)jj+V;?=qd;q9fD!0dR$+o=%LP-R4 zd-s5ad0bpcvd$r>)TQ%AA;b_4pMHl# z(_maJDgZJCqXVZwD|U~6s}}sm@vkPWX{k>P(=3{%KU_bEofOK`UY+w_9ypYlH#4W{ z2Yfo1y$}wyBMXBPmjlMOiJnXw{MAZ6j(4AT^kTB`>9D^Ye6!QSF8~Uhbei|Kquy~} zT+{FMO*Ny69`*GtvrYZ6k6-4T0cq9zUD_OaZb&Ii;|M<~pV!&h2the7BM44{p?s;iGKDMo^F( zv(E@1wF(-LMho}|SaX}iL#cw~4GCE%AP=D6D4@AJrd_JimtZ1bwA`W9 z^2B@%X{@ZON~)f^N^P5>?DC2(lZsDy{m1Iamk(pPNc@;GPS!h_nK9ag%xXvC7!sU} z7neIkH6mYUM#^Ee(Cy!}Cjq@xah0#jvf(FDjWqOYw?gHnH;(mWwceO$JHqdMw@dp! zz~#;hb`1VJWN+_dNpT0yv!OZUGernO6WQb#n@|}A#sipRq!C1Ni)_&wo(>Kf4FN9t zn=2*sGCsaXO|*5ICj78U`F*@?`lqerlD#`bJ7?C`CSajNjW7%|Wgy)Ojz=ONHdIq; z9UH{D1=MVx0LElxF<6+dpmdT9n!t`YMo>eIjcL$n2KW@3yZZ9PB@Xz=+JQbRb=}dq zVk+l_EjWc0!kAYS)HYwrBVhMmJVCk{1U|)u^|ol;`WaBRCyUWyjc#uiPDo_+cR#4T zq~x9ucaT&v;cO*O9+I_wA(#=#b9ughn<5&tY`*DV+g`y92k(3Eno~^ z!m8-$!85*EMRP$d97^>D!c(rZSb2l|)slpa2ZI z35iG{1r|;%PN$_Kxe*Xdi@YyP=I&PzQk#hULmWs%--({D1EuVu_Y@15QlD+TbH3t$;3Rd{rL?FH~Z?lWEky-!L%88^N=)UmTSPD%2Gv`~fCXI7qr zOU~;h=~+2yEbu?!T`o$#`L2FL%CDj{B-B<#|)#csh1#y zH$&;(%Zco<)P#+|5>*#73l+Ilfn0H`i6)Mcq>xU;njI#M^z+n^MDUy0wLFB!RtTqx zZmSYMA3g=?;9V&*YVcLm?zfEdBulmh`aDhBnd~dmCW;Jo3-@hCQzP#1JXm{9tA&a> zU8J>NL(u6>{FS957siChQY9bSa%<1auFK`F`O2B}yIgt?ro2bI^~8xZrq8K3`b`|{ zjK~m^l>zZi8+&hGp5G6DFLdRLohJp=Lncz|UES{6$8(=?E=gzY1}1wZ#!OpULv3m8 zuOfG7YfNjF&5O|s+&WNtS?1gw_2QjyaATlvasSq4OIyw_dbJ@ho_etXv;01YZ;@%G z82I1T?oRrTuH0CSTP)Ryh)Dygb}XZ#pwu)|V0i=inwwQB`c2wbTUC$0&FO^;+1=#^ z1(Hmnk*j9oX|3+QrmK{)Ii4W5Ve(Tn}UZ7+cL^kG7uSl+%R3I?wR9}!v#_vU#1j}Iar zSAg{BT`9RGGD$$_WstW1NGFOm{ffGAh5eV6F7+#l7QBsdxWx>HoQhm^OT<1(eNr}> z;q<&3?SBpl?GC?*S$>z@nhv_0M@G*k%wwb0zsS|tV-DRo0{{7q!HvU*CrBh?GB;7< z=#F-;?8l_uw+eMgaof#D2rq9{MRC3_32}0>cMn=jjByNMVH31lO zLl%TuAvH)tWkyqt4&+8iLUl3(5l63B#yKfZ4>N0@I4Hmg2gWhg9!cR%jkBJSzsV%? zYDtVr@mU8)-?_Hta*M*27wqrr95plS$B~y4bZj#Rn-6m;C8&zJM<$J3c(dHJqE%Ic z33N$aO_e?&)+kgU9MHafbD=f`pA(+k=<>G{4`znF9&*BbI_0X}O2p0cIWZF%YxOo_ zM><#3N@D%Gt+l*dy82@^a^`Z+8#S2B!qHrPtC?#vQXeQ%LlV|Dbvn?aDGHn+bZ>E4 z(aR4~p(h^@v3T2frQ}JWG(k9f4V<}IX#yJHG{8*_&N;0DOKvnycQ8GmMQ?5{ZH0~% z-uYJ)#+2;Cz}i7~FgKeXCGN6su2|?c#&UKY^g7V6#?0)zayA8N!7&te`U}I~KGW%x z4Z4uqiX6+wMv<HnY&|g-OaCSX?!bG}f3f z?Y$9_)5@$E0|rWEy&Zay(<}Ysy|0ev+D-2`WyiH!MtavVmk&bnlJLjHCA}CDYRgI~ zcnqO<L&J5cR=?q3MzFZ%amMl@5U>dDjHX)=sKaJ*RL|J!t$J7Lz7#h+QgtPuBhf zXXxJ0tm6hHCIxeM=&SKXJx9DgPt!lzlz{{Y>Y>8)Jo*%u;&6`6a#$7yB&7b0jxGn& zRpk+Z$Q2pQy!8`RhtKDKYMw=U=PcxalIp%b3QjGGz;KxX5FwN~4K+>XDF(12QNhH+ zKoVsa=%hfCaRA7@pcSdKvRr_iaY0t&t7FB%yWVlXt)%z?%?pCc6}EoWwzV`c84v{O z7nmnJF#w)uMJZJQawahn8dUekD5H%L%W$XOHW2`B%Q%y<*q*%z6I*Pki5^W!nvDx3 zzZtQL)|B9c*GY#1Hamk5tWmzKJ4;%d`BoH-gicPF`Q(k&1O&Mxb4im0gitnYn{1)c zQY8pm6Z$Gqgi?wR4cIFSInWLmAynVvm-h!0P>U5z%>_{ZAF51_O!A>3f9KjVmPm zao6)FWAY4j1xpn<{HcvkETod9&1SPsFKXf;G0_peQU+;p2ex(5xcvf z=IgeqiEDaAZ0Ek_y6O?So*ECM4tSft5wiEjWEBmGNl|%lMig<@dVbNK&8v{{b1X-i zYL>k1j(6tM<-yv0ST7)BAuK5(hQ_|?_c}@tk@uNm%*sYEzCn%5Uq4|Ur)<$XT;8jT z#G`hN%*89I8zW^lZHh^Yw*dr>;!(HQ;lqoIG*SXw$+n6)8JhIp4jt3@)Jux|VIp3! zNAQR$8<`{*f;{{-`Pu8z!+1IjL}7JV2f|M}IW!A#cq37Wt8PB7tLdFek!o_waFmHG z9G7jbLLc;Ld$EpvOqvYhgB?-pde8VnP@4Lnt@JAK!)I%A+Pdu2*!v>1&{GBObtLWV zvZ%e*mzO3|hbG<>j27NMif+>%X9L65^|s$kjT|AuNJF_I-Q?}xNfS!MzH`|6rl!~6 z@U0?Jw_VmHUc$Q4f;)syoQ7W&L`FhlArL_#@v=HGO7|HGdkrsx)0E-G8@A;c>v%y% zGYT?AshgdglnNPoR~8t^Y&@oyGhbOQW!4wZsCRWEVi zAkB&wswE*C=+yb%X{R^`;q6i55xriUVG@NG^B+6rYuKtT=+{ULz%w7yAVR9?GXS^l@0lS$4t zF)+@#qwk`8zn8smhap7`nZTdU#{I;kHFEFM3U@bZ%L`%0O)EQne!Y3u)2_Q7s3aQD?!U&UN~YkA3z#Gk z%_C5WM<`@$jQ1?~=rSV%a`c%$Am5q_?Wn%OG!j4c^+{QmqKU8$ER2?3{wW zFihjLTH10Fp9X}LFvq`gl#N3DzXomUVrlNq0+{K!Nrq-;x*=E)*bJ;zhldCnZlZ!B z!O|cQN!sz55@_4I1EZkD;)U60c8hl}j&=K`N6b80rCTN8;m^B0e2n-6S^yg^vQ6D=a-HJL{HFnBjQAvT4g;5C#9g7U` z&G^|muij*B#UZv=YAtXI*iba@QIH@O4=YY)a8&;lG` z!Rfl=hoGvCb3lCSbOE6)8;o(l`JN!lAdGlNDu@{{j7Ks85U$j5E(!<{w@g@qZeUyx zxgBr>PHgCEJvwi!zD3=l3S6n0g?`u5%|%~4zR3;jmIkHL7&5Yw&tGr3It?VX`UhM4 z-%qs|^4or#8qIlQ@NQyveot?2da*-|eeD$q{7$H*h6~hG!ZAewNy!1EQ)_0}eNrM& zp%>CTWrZy3zokn?ZC$$xW3UoT!qfLh5xK@^?^%cu0Isf?BQ$Lco{}~6`-8OoPMR+!qbhZ2b>hIm zx;QXk1w{(kHFY7XlErW`>q%hg+)ZXkCZz&HCzW7X6Mjjrkm)B(FneC|B})D*pr{$j z4f`uTdeE^Lz4Frj`}Ugh;^tbk)Vi>#mxp~e7DcJ)@QiPJ5L5G1@m|t;firT$6+zdQ zwc=Q;#tX6;Eqp{JRp_9oGejEk9<;zI@}XI2Q^l?|Q?Wiq2Qq%}Z%q(qGyR#-BS^iD#I8kw$@ zXRfoVh!Q7Z<&nC=FOJT*&kfa88DXu9_TcB7ioua&(*Wv-Fm;=@mPq44Nf|I?sH2#v z$-_j}o5TysCisOaIxD6MpgE8y^|%_TdO-C8sQHKRdYAQ1sp(AD9<^a-h$qWlt!yc= zw+_WU!YzlT*OYv>j%mY7?)9oiAdgpRH)?%>NX`J*rJW5}_1O@avl#}*d^PAD%K`S$ z&uC-Y&Eq$BC5<+;L+QNzH7NTA@0_=bzY$P)B@RDaaDEMYiWMh?aHn^lnwNQm2mK1# z-Xls$>91YBL>v7a@zw5l;Y1r?VHaKMnQy{S*rt#?C)!n&BXg557#mc40L317ZS5ekcyHm9Z0)#%b>kQP$my`!}8Q9@5Pe zmZk5c+{u-{VuD6RUx0c<8}6BR;+@*4D#!thiBZf~QcmLB8IBF*Q#=Qea2cW50@0x> zBmPG%y^tc9yDr<~j~ zi6eNyU}y=1QZp^IjEBDiLN2`{ZY)*1TQXBOV_ThX99CbZN%1nm0Dm%f(Ob_Zm!|Gh zji81)q-m|6<^trycS-PUEW8a93DI|r9-fSAb#BsYvp`ZOLaqy~&xnX4EKZAZjE(*o zYssLHmQ0Zj<6h8IIyrIz%3>WL5OUsnJYXD{v7Qu$sz^N2#Zc}bS~IhYm1kwSqg0n` z)tomqtnBtxgxo*rlH={El1GSz+|HVU}0nbb|mj4z{KE z$-P`N!-?kYPKA|Xp`%2dlmmPO2^FM(+SwPzg$WoMd#sS4AbfkA8tfZ1>4vDsOz%1R zLM2zga9`JVVOw-e;87>9j{_M}tnE;Cs5BCsRF*;g60kvV9DAfAOjH=w zge%^;N4~*MwC9w(4)>XYb$5GQU+VCTx>ShU*q!#LRV&5t?q!r=0)iK&Nb9cb(D`Lk zb@8;61C&IQxjR+e!y-Wi$l$I8#IkRKQe`qi3CWt2vu-=BCQ1uK4kQd&wXEVC3X z1;tm9dXmbxSzWAzr5RBKv63*!v6c4q=a}G~d$;ZFPY|HnsH1dktomu{wUPg64!5?p zyR|}aU|kGoOa}8qJzFKYQ;e$u+cmf3L$?8~c+|_rXXj~KcMy}wUe?x_gf1%G%jxay zN?lI=qa}cTz5I%-3jG~Sz$Bvv=b=Lrx{KG&++pKX>aS^1j7h-J!J|rYTDo-=;8%C{ zc!or34J2^ z7q2~t-Fv(p>dkWO!%U-?{-reyuxYvJ=hyZHCuBn!@#{;BFjYNkY*n8`B)v{HkDF?Hd?Iq>*q5tdnw2$T>42iRYC|3v^sRVMe z6;-F_nQ|@Zy0KLzC9-F0%Bp0=Mul)Ib|PTXJflK+3oHa6sP<7w@QuHsw-rJ&}s%S~8a`u*%l3Hw)z$mtD?q z#PqGd0jSr`oaj_J;t892Xn1~|_zYU_udjs(NPRJJ*CZswXoR=xVqn83GLm!C$&x>@ zzN!yXp1&Y8AHF-?hF8;_CJ(P$g=G zeT)*-m*aPp)(FwSo_46FMMAy`j;7|QN3?ek#nA@%118j}%f`@_ z>NUP?4d8xTdhY&Xp>EAN)}M1W9QW3QqTn^tG@6do65J5X8QrcqHXBIjNs?5wX3J0* zqqrHcK3vV9bgHMtIZmL?8z2o2S`w+^I=D0GCEc|$wz39R9D`=Qbggu^%VNDR z(`7HkFpd>6p_8(z_d_A66=2LI%n7N%35#&%5yDj;!|A5H^YX=xm6+8fs)GjO2?@v~ z%?w3cYllvwW?lY^(bs++^(IWIeiTN(;*Gtm6MOd3qNHkC&f@Ef;G7f?a6wCE*IGgp za-Qs#x~VD;mfJiJ8ASJ+b|JGuT^z{)VjdRbOrerY$Ure;Dvg``k$nU0N)C$RuNUCP z^$)ENB{F&<7R``LM)mJZYI>Zl;WAXdON(GL<7`~1{$F8=h0q;6ppUz1SAn?#0RRP%LJhRjQ?%`6tD|~R^HI2s zx%1;*T%+Ao6f)E_9ymV17wR%JBgs%VW-8zlmcc4^&9nMks~d3e-leZzII&YmLd01J z_TG{jkDj+=rpGD0*p~b-sciXmNOSX9Nvda{I0`Z#r8|#37>(h9dxlM-(`@RE^K3@G z)_k}{1U6`>SE>;x#DQdz2`(@|SLc)R4!@c(_~k>Z$n|* zOBNO?4s*Do4Z&&Qgn4eBiPsN-(YiPECERi>QVlrpO-N}<5ereeNiu3tWECn-HAlR{ zm9;xIF_U^-`!8z6ZXO=~tU1W|YknH>eF2C%cg0jLKx#&02D2l$0A^frQbq)43b%`p zYW9&G+a8KFpzYsRcNG*)@KvaVI&7ON_5+9<##o}MBB(v9DSQ$R>6DiRWTn83*K^Bkd{JPZbhz!D6^XvtrT)9u*IodP_vLSe7&?O zd=BCEalo3uNF0M2w#aiV9Ra%`k%D5qj5M>@_V(j3664oiP*w^HdE$kD4atfV?hZHg zp)n54$n|3*$jK^aKw~;{uHLS%m1|c-E?s%1Cw1q+wHOTEd*7LmxaL`?g#U6ydOWEGagQG_8f>O0zB(!w*xu7mXt?Xb5DHO zECiGaQIax|^CPwyvAY+%gpjlXcwac=OljV_&K*UW=a8^U%jPc`FW+J`Ppec=6vXX*q78?HUY40=fJ~?^diP*sl4XK5dHf8YagnotSr@2yr z=FrlH27vC|3xXP%>K*OAN-a#e=|&^OCZ*?sLe=Ryxit}M%#H^Q4pkB8`CSegZo>JM zxma>2lsr8gLPphGaSq`(8 zf?DRX8_hGw_&HIl)7F-7PHtvNRY=_os^offm2e}xT_v&N?J+-Ol%=F zSmx>-%>63tu4hwDRo={}O1u-G7!gCzqvk-VCtdsuK@`B2*f$gnK-!4?b4aXPfjb1XdG@81g5U0*lL|AwIYNP2_mRA875%l(Y9OtSBH*tTvv>U z$t4UjWQIgqSkH-5$F7XU%-5Q^k&Ks)_oe4`XkZJ0fDQE?;EMLv%MOSMxnT;O506>V zjp+DM856E=MkyM_$P%|s35yrE7qpR=5QOBWh>cuw2rLplPV8berCd;Dl9M+|emAD3 zuD!Fan1mp9gDt3Wyzq^w4_;&38Vv3*RxU$8HsL|(+QHz3)CJUMJ3=8@Gi8B>DI|ud z$PocPk=ZgFq6O@v0|EemH5~Zg%K9VHTHq+R8>wqtQL~>`iv;NkN&%r2?QJM@eJNZ4 z7)LZ1OjTJ`ArH*cnjl02bO+E3jEU`v@Mx7NdJTZr$WJ7+=-x-S$ zUp1oZM^!9cXKEZPok=oVsYC@epp`LJn3Y-JUM$W8F2(5fpbEu4oFO{HM)<^29;< zJc*E4BLT?Bz#d`N$2dAIvO+=P2c&YM-a)AOF0ig;9VTLQ{YXN!e27u{Y!k_r`I_?R z`dJ$92$P3UvAAefHY${^MlL3ZB3Kb?h^myF&B{&7y&w|-mpus-ns7DoH1|xE1ZJ~S z+{Qz@X;L%QsYwEZ=ouJGl}O4!bGEg0*Xr=;nw6mN)P%yVeeCbK%moZ6uqE3;vnC;} z&OY`JHXxTE0w`s=+~`|UiIQrnjrLt$DqCdOuRWKR3Uhgu)bOIo@>ItYT6lR>loOTx zoR;l7z}i0ohY!a?CFY~GUpY*;dd_GNJQ&&mj=2epRPK%*FIQS-uq2wthS>t?41>7r z*^qxkQ3x`{!8!%#!c9G=2GEa83uNWNtgtI=MSsCS(+lq5KAy zRP|RjshJRxJLJ?Bz5MGK^uq5>!(|MYT|C%! z+aSUGXCf(Fm9#l}>sKH}V__H`gHDoNArpCUwyr1T9+}-Zv5zb!vPfX4T^!AZEC!aP z^<;TWFI2fG!W6j2UT``{Q_tbX?|zb7=kdPVAfz4=NR8WQmN@U*I9$C;9j?K zVR&cZv$kUgv6L|>QZ)h~?b}02vHP_l$HEq0LVcDr7kAHWW}71Ro15TD4Vu}&QUHMf zYIAz9B^E=}rVtn|&STASO`6UD%|nN4*RB?hsWMDS0Laed7tceqE5eyBhDr$@<{{Xy zw&W$(eG-JT(1z0{H>(lIrQLq&?KYcM2J8lpb$D-7!i;T#-Psgld)Y~3;-QKR)FK&X zmzbo@!=Qe6;S|qe5JVk!)!}M5C&FV^q$Yc&*hh^svI!5XwTYFJR}k!oT(7@INHyIW zfiVIqXpx-{z|Cw5GvM~Sno%ksZyU|jEygY)aE2!X(U%3NmNeLiw0aylnUj!jPrX#O z<^n9*$~p?i5{uMbVO;5XIy zr@QlAinT{x(wlK)ErQ^w&^+u$ox~(*B)Y0w27rm2ilxU*;<4?IW;S~7+ZTb6iNt9w zpIP0D=rFOgq4lxktY8SP-8EY(*5{Bwkppr>z4`OzK^qKXR%Y%e#g}nukxRws%rZCZ z)k2&@QL|ecX)U#XDrdA}<`(+d(k3E2#UdO8+yIBsrcBRlg~AbpQ| znBRi0%M3jY$)VJHl#fO?^qPVNtanVPw(c8@t%58ml33*)ouG}bAG0>g3O`*)!9M}a zi(`<|SJq|Rnb3?S5`|Ru@y&k~O#&yaY7MaHfF#l0?mej)Wu>DF+I#4fsFeGEK`s$W z9^5#lhDD=J2aJf-x%zFr9F~&K9UJv((wAoJ`@K-)gMta5LGzh)OC5#4JS}^*&0SGK z_eZIYIWD*|6g4woTSIa=;ofNnE*yz&_e_=eTAT5KueC11qVo_^CiF{s6hOgt2-~=Y z!Sk6#rDbSnKw|up0t!n~2ouRSz|~#Vx};QtSt+Wq7$cu(GL650>5ISW1^9dNx9qP< zk&2TaYukrcGr0SvR2M9R+AhhAqhqHg(l!JES>K6CGmW$3_#RxonCYa4CArfS%E!19 zFUdlzmNr{@+Km7syD&%T&KAayNO>lPh<^Ax!M~Clwyfk@>%k)7cR(xlwlT@i z#o#Zhh=zzQC%NvBWwJ?-L!A_mTE+QLprzqUpCm~#i7Jy3s(r_K;yvX0c+mr8TWZbR zR_t$znM$EUfPIn(k-FWc0(>&yWYNbUx3@*TP4$DZlwskahGt_Yt_Y1Lca6$ht&Y-J zwohQ}gaet5BK;RNBOhr^^n8?GInk4y87?_cQSaJX=TxXddPbC@8FUdnz zLXR!%OL|!^V&03ooU_RPr4;gLm@XL-{2??voj5m3I%a;oa?Ngy_PuIm!3`4Xo$r)f zc^)!FvLOeE?kW%@^eskc9?<9wfxKCf6k?QA99;W>~j=d-oc^%;i@c1zdlB0221HmfB;e*}6qPo0RAdHH`ZrtDzAYgu9 zc-{)kbv9vybMNa;d=w19TTt23TYNu5)W5pcxDDZVYi-R=bWpe^ou13kn}QpD_GF^&j#=yG^aDsDI{+18`C zI^u-1K#x)O(C!sI0!8w8wT7@vQ21Xi;Hd<4Qvy%rGV8vnbrQUbbXgUe>k2wTf_Fl* z-h_f*-$=FlEM>c!8s9ZZ`mMq|x66379d%#n!6_-TU^7*pn3y{QIvX;g$|;+&OE^bu zzd4E6_Z;6wDsgncbOY3=MYQ6?+v-(3G=^QyrVa;e#gcsn?^JyZ*VW*tLf(gNF3Vp! z6A=Y8>&lorkwUjHGFAjmH%z+R`I#yig7;Q21*CPSjAI&@sj=(Qg@g}C^E3AD`mI?Z zcsK++2w!F?Ugk=@<~v07>-H%jjIT&SeHMXSW1jCDVztAfPvCX4+7ecubUHV7Tt4uo` z3=xHpry*s5#%|qUrFa411mk+N=tB|6Ly$!b?Uv`lVZj7O4k9xU6LHj#{cU);cR!yV zXDj(ENtiO7UqMwz-D3#AF@pPzWo?(76b9VbQ&XZs7SLP+8z(9GBQ`O^M^ay7CCFfm zfG5b(LlK^4UC#9!r_|fEZ*wuoiL1EWqNBd=nwQM-8w?gfqWmo44d+e0r#c%SRTtDk z1S5eEVE1w-nZ+~(eW-(>PiKA$su&SkiPQ>&*Hw`)7J8D`Vz_9IhY>fB?u}CEM?t-X z2ceNdOrE{!kh3sx>16IP+}`5zrj>{rNO+){0Hhum-x$y#j7$$x80}RktluQ@A9ZRWiq!0fZO3(B^nacb?0KzaJdYs=SS?yP-c`oXOWUWN zS0nG?e)Y6olc@Nqga`UTkCw%H)m>%W#wDA&nAO}&RRq@wxlbQ~9VqMyQ}pg;QWY zQ#F^In#ffGTb$G=xTR0Ex4q5Og7C zRmO3`9(U3qH`I;+Ft{Ou7bA*bvI!{y5LL?tft0f}8=84gUe_$l%#}&vfTe0H=eOk) zJ-(6>MM3tPmE>ufwRj(-D8bZFAm6JsB3Y!i6Ao}l0Va%zz|9zRsVKM>>N?z397hs7 zJ*8pDwP1>~v^{V;gr@pDUvlBAhD2PMI^SD#-)yhCr0AAL#^=sgA}|nwHLI8a`ZNY* zl4e5I4Z||$!pXc=iDgiUIQv?RrcV;Qi5J<9C=yxj)%ZgabrF>tL=j@98FVS3Exuz= zqR<%_kp+})ByuZJKo}>w94wM%q$?=5BfArx%j;t>sO0dG0qn*uxxyz3tBDM)ORWkB zMzuhikhV*?ldV)yNudbcKyE5YX`q@W;ZsAXsvezTo2yu0n}hJkWs8x;1WvQg`P^6) z1R{;9krg*s+$k!7pp7L-fY1{JE0xqpu4{(?a#3AzH6r?4t`s6nkyCcSYDGvkTYx|$ z?&|`q$cSV9x6$=Czlw9bsRetu5Y&o|gcT$zhoJx{OoK8>R0=+~$KZPZq3eBJ)o};( zcmB-&;(<@E0HmNo%w9Lp>ShEMLR|%T6C#h=!u{W|i9!R^>|`i=PytmE5%eJuL*)70 z&&vFt+vj*Ve6rqkSvb~ycj@Nb=;e9$`kV^)vFp<}%JuODN$0O>nSYD)o(ZG@)W&uIa6c|w*WCA8U49bE+3Uwd>{D@FyNI*|J z^EDoZqZj6C?$-b|nf?U&U$IIhcp2e3=mzzH5FlU>KBkY=&i$>o(y2v%as?ci1p~l< zUzh6q`FgRT)9dM{;C^SniTeum`t5h59s*q@eQR?|H6x>dbK4H#rQLAGGvNV(U&s21 zD!WS=W`0V?VIosDVOaq;1wV}BIIgC=_tyDKex$Sj06%ReQVkZ1qRR(#aMO(8A!)jq z=*sgFbpQZN6~7E|r7k)MqyU6Pr$JTIfPSvC=T2QjLA@wQt(U3P)^F*ziKZ@9ZnfOG zSDH?q@?5J*H`9)8bT+)yLO)W~Oq9xNeX6(Cl`3__6$ z04w*pyGg`;oGPa;duCH08;C&);6Q@xFL%1@`|0@9?Yp|-nmPOHTuK$%_Aw9@0@on)){IszZkfg$3tt zE8IC6G)utfpq~~224NGQ2q~YLAxu~$3Z{%eAVPw78XG&C-SsLQs6ZkBD<5h~$_5Yz z=t&Vdn|Ubxm>rx6`Ot$@ED$Ofxn^zN>v5{l!Ir|nxEZrV!<7l7j1A11D3<2>K!w;M z6hDdt9k2As0G}E{gTu=2FEXeQ+(4C}1??eJ+;F5S^k5gggdchGI$v|n|txK^G&|^N;nS;ev2Of9#g*SAA(~tvXEWuaNh5j`Hun9 z9Ei_wdmhiF_&(3{?HD5j7pz|*@7?8`q^Kb}vAo8XKWLUt#jM7n|A!j@`}T~>_kKmI zZB+T}PGD)8_h`0!=Ad5}*Eg=>t;BYe)*#3}SinL9mKyr|#!%Qhu=y6?!{JA?R zD`&yU1OP5bkhKbw+A9*L&b?XVfpgPoF>y1Q5EQIQb6cPbB@p2xyJcLDt+GO~Z99&v zzFuX=M7Oa z@ptBc54Ks+^KNrMdSg@|j=1mh0LD!a(kI8Jno0|g4m6r=kis>-3!m9DK0G1@w{xwF zXgTLDo6y|Ukxu8m&{RDeuJ`4U(;yezWCAG1^ADp8Qawd?I^zXc$tt3(RUKn&0~^8u zAo$QJt%UmZe!O3GW&88`UXlIDzP3iH&<;(iXDn8=TNaB=^09xqFGqmvKc|EinyEI6 z9_JU9QJ6Q8zX?xdxeS;}gqkrp8cZ>i9yuFBTB@9YvS7E1tl(`h%GsFnDu&yhv{Lyv9SAE~L&DYKj^N}RTgx^C{7??J+vctr} zz@CbVYa&*o3bvOXEKUe)h#uuO9FfcGl@oEG*F?ydr(2~?ChuIvRiSu>moV&=4kfCVag7=45xa47cDPtsprTGBIIBe&eeH1s`C|`sOu-7;IDK< zI-%C3PWPS`fccQpE-w@%tm0`2xH%<*tPaSPt#kyMN>q1p-YSu&i29*iFAa6PHx^W8 zRJVEy-DY?@A=o~FAeZ`;(zaZ%u*2_FZ6^@#c^yXsNil2_5tL_^ z3N@#r&q~zgr6=D1fATjAt^)6Y>OHSR+I!D?;NGvl@Z|V0;aZn_#+D97E;%fGYB91* zn(CGEELW|i$t6XbG*&We#r-mN>+C^(R;=K{tbwrxc3eZQX@F5fgJ-31EVC+xjiB57 zB;Ej5TAML;Y@9SGU{zW`l86%IwY+U&_E_Rw+~^Ycql+rm;po#ifP{=vk~|X*+f?I% zp?cgwZ!5+sz&o5G7YNIeqphwevO3k|_otTAYjqGv8;wql5OvH6oNd5KcY?y`B$sWS zXtUODH=ZW7d&V%q#k*Uw<9E;pUpgMT0Px+!c+s`^-S-?m@|;9Xe4w?%&C~ff@};fw zk=~q-y$PM~zH-wE1yXp&F}$)RYECh!AYymL!3EAXL4!3HZR?RY#{#L7Eut>chk+t+ zzuEE~b=N`j;p?FI$!w2~dD^Z*y@8PHpz)7h*oPD6$&?+jp1osuHR0le6Y8*;^-1tn z=sIDILM2r+E%0(Im5hUNCAfRqt~Z#4mqb>|v=wkUt8y1dFAVd8PM!_<&fgFtSe>l99k{bF}$FNmDpiv>xM6@n3~@Ak$WT7vRctoy~6Fe zdXn2Q+?kU$4s`=-gU*q#ZY`d`IsldLM~Rtj4L>RtdnQ9QWE2B@1*%}UqM?+=php?E zSRk?j?;WPIx6btGcvJ()avn-9UFV!07IKPQ;}$u+Yne)EnkbT6fZb)&XSU|ZM_0M6 zl-RXyCciGAr6@m~p5+tDzu?LL+<6XhV#=gP8#lo#tl(vWP zz>{Oy#y$%KEAHhZ+!yyDx6a_ZoLD@L60X){-(KtuZ)~_7!=dT7s!*U4DYP~j%uE8x zFX=Xl&O8*MkRT@=MMZN!-a^m=Gqz?dfMROa&MUkQN1LjP7y4% zjTX>CGMF6Z*)oRF6`|#%KCfg{L5>P2IPsD|U?OC>cvstjL9W6>a)ZRyIn3m1zJsqo z=mVhu0B+Z|nefKEP8+!^@WgG3`KQiH!e}qBIL;Yx<+ZnI4-7EZnm!(tcQ4AtwUyfy zMmea=l3v6|Ps2T=OG4oFOt`$PW0)ZuoGaW?o`La!J=}k*zz1CbeaG;H*FYTrd=5WH zXQ#v3;*+SC<#}F5F&c%os@=c0Y!C(-SJ`wtW(p4Y-0tINxLgrwvgZ~T!>hlm%$npF z@l*t?U`zl9nFEsmrpQEueArcYo>yJB)KHB1e*$3tqzc#+2owSXs~~>l3blLgCQ;gS zwFGfz;`TNNS6;0CC%&X6&i2V)gcpSfOc0`Nm{f3K$&7k49sAHgId;G4sN!HSBYm7Oi z;j42k>r1|L+8(m6tqux-qZ$A$%u*XH!dwigV?-0_M>eW|TLHk!#~-q=ZsKwQpovj< zlieP+smC1b9yX{2r|}}gGe3l$k9wSO^yZggn#Z7LbItHcr#{p^{fdBEQ!s0P-L7`{ z3rd3y8a4>LxFczZyN$+OHVGa?$SY|aXDl|+8VDc=2S`AGrqi*K_ppHP2uKwn5aI2l z0GwVGTprqi;eQN{Le$@-U<#L?B^6Ab8M)^e8Fv2h@Sdfn22$09}3#>s326@+F-w?LsY3La6}x5ULUg1PyMZKWY7Y z91g!To1xe1`BU6}4T;H(k1AJ40fG=d?&d&`3IK>v3?0U#nccwSKq-3&hQbg|!WDPi zK0nFx{`Ob%{T@f9{NL;TKlp!EUpKz@U&n*R@z?qukK*-zlixr2KXak}Uv_WvWxCgI z=zZ*;>HiMLXXki6?uK{sehAm4@%f*Z+j`%p`M-7h7=3&`c8A=fnfU!*vLEAj@_K(7 z&(Tal{2*@|@GJEa>^GFh>w60(1@^4}b4daHkF567`+q6qgX#Y&YZv(7!enU1_1Y>; z{Ixzk_6EO~0ho!o1}E<7z4p^MB*SjbiP+^5`*2(rPV?jT5&UD}{q%KR`IFK}S7XZ_jL81W4HX}w0Ewge zzpIn^-m$Nz8#TTOyP5aMAJcW|8hWa$IfRxvCp0LZsTVeNn-P)Vl+m|qPmMT--jk5m zDgv?~6^9lkK^6IM&94MYTZWlzwp`4e)E$rKweQh;r!>fQ#1=S)iP|={I6-^Tl8B$h ztNmz;3+j(&@t_2X6|gu8hpw#!ya>tZ0scWhw$QFUD4KL9KRw_)Qg4qyW2N!AaBk+1 z%E>xEIZ}q0D^&v1VxPFQ{(I;}nT4rlozK>uWJOuf=3z-&qbwTZ6DAGq&p%lu9>J}P zZ}I2?*lb`a_fCDiV~=*bqPqYE=eWvH%5Af|49c)HBDxaAY7lUqho5$!QuofHYRQw^ z-Ev4^HXSRiw&#)PM-BZD6y4^Lg&LIx%Tnrv~myRL8PD5D8CAwV6~E(Bk~Hf+CC9As$mzp*tpkR33x zs@1V6odFl1OKL62kIzYBi>p~H2c9GIl}(C*nt?s_WlZp6`!ZAbz9 zOmTaqL=-Df0|F4Q4;Kt+b)dypVbO$496?72F5H;sxMUW>3!M zoepiQEF>}|yNi$bPsP(Q@GEuLh6S|l=pY|(o%i2c1cSKD3Ik1{wPm3d`aRx-<(*f^ z7`Jl1%{-cC*xYVshR!AJ7%56k-LxFp{_}KPSMW!XQ8@I-&G)h5b#f0$G48nl_io@p zgQFq)(dp-31$2!mjQy_Q`wP%L9i|~qqQ~Z@Xcnzy*G249-cIDqjYC41YV@DQY3VOY z)1ejD?9J0!|D$a4>4Z49deP`DuaFQWK6q_$BdtV#7#+>)x6#^XT+#7fq#c54%8-)N}i5Lr~{C*@-F{ zH@0NglW4IS6iHyhZ%k7t)d6o0&kBP(WbuVAr&a7|p!jaAkh{;ga?+)uHktT79sVRS zn(dA)kNRG8v=mg!zg%nD^*Uk?7FKN~QaLR<05>fhVV~E63d56EoPl-;sqqB@b|6Os zVRg2cZvIbx@j^*XOiSFnyLOdloQjrUsyxiQ$ILtDRl;i&}OrmZ`D;-#1t)+wh< zQ}hQwAD>G9U3u@k^W&N8W|4r^stLJV-s8ILrtR5{@NM=zE}jUC4(NgIt+TlvXCP#` z!+w-Whe~vo5On)qy5uX{2+pjy54nzvByres*-3W4JF{u0VoGh%o^0JnLeq*#7W_2! zVuSE#s@B^thmS0?b<3IG@gpdD?ApiVRt8(u60Jp+syXv#Ge@o(>kARj9B8{s#144V zHzhwC;L*YG@nQ=*V{HtRZ6$ zq#q4awVf*m=r6XXUBNXJU{_j;fP(|k-881Qs$4}JR@TJq4vi_vu2>;fkZ_1RI>NrT znnFs+$W0*}8RCDH!yY=k7a5a1X8C?~9(Jb0{)wR^W|M?6wOGiBl^9kA+9D9n(bbDD zh6)lZ`w{u6i^^vhG-+x@f{{(FI6x>-}%+dq&{#^k@e`^L)aZEca)67 z)00cLV+_~*ZPpt0&N&}`MYC#M;4LCb&IQg(9TmDTvRGBlhB*@y{bL5)w2di?aVMu_ zAlilJ6!K2o^FA9!)oy&Rrrm7>?61B?*4>4L4L~Y$Xs0}^;|ZLF@@anJYs@R^-Y?j1E;3t27drPQE)>r$%-6l!d!20?JK0c0&KGj3RF%M9 zg^Z(unNU4;n(l`t$)11lCbQJkm_i<4xVE+nWr$cE(wuYKdW}gm2G$+p@mlsr*J3qR zt2*L$oX4UU-A*uT>%6Qje4LSN#YiTW;<)S%oV=IcGl3G!7u4&aIBaGiPitOFj0SHL zrD_n;TLFZo0z*$}JI0VCcAin{Tt<>N5c%wgGhPfK3eKTudE-x{F(`8KTBQDKW!R{t zWl=7#IuUqh6{ZoFikdxD-=R^a*nHrOqR_2Exs*6L4bFw z6h!v$F)%9mm~?eDT^(LNACn5^z=y$w4j)!of4kO9jY_~No9}*%!jp;+r|QDLh6mV# zTMAJMhDLr5>tNNa%-O}-zh>>7owN&(FU^1h1~2VMKsEk_k5+rwEJ6Y8LWVQ|KRy-m zdVN1H^G_oWLsOscB`J9HsuMd{nD7A;`D~L_US^16QPoC#Xvw^s7MBuCcdd<%$&%4s zU(?ckS3ePGT@ly$?|BzEe(jMdWoUnpKScHVleq7ss#d}JIS9U4bLaE)JKdr{Hq3|U zXF`RlEnt~fXzw+{c{Db%Ji!Y;j*5#Tl)<6)3AOry+0sw{VUowyD;VxZj#V)07B;_ma zBpXT1T(I9@8*>+Uy3I=tqXvG;bw?n9bmBlEjen8B$j{eot0P<%5r034B_<(8e&=Ipd{cvY#O57~iC&J>HVtBGb>FacBx`yRZLkL$@@ydlB(emBCFtI+ej zUo+rVo__?r;PTNO&FhIJ$4QPsP2-!0Z*}>GZ9YF0)lbD-iZtr0CUp{MsKP(IhK;Uj zytH{bYrGuMa?Q_%kS{Bk1a*LZC`#CSQaRQ0dC`g!%BxPuOwIT(DZxu^xMA;c!kAze zulzXu?W>|6xL-O`!F-nX!BJy}DWSH^@m6|=jO)E#KNan^cvG8_yJSr@$bMoyA{F(h zCtFH%hnhQY$wfoa<;v!UTPx>8i4_4Xw!e&8d9#(^A3q{I$|&}iqBg$3btfCj-j_LH z1kyNNxuIJZK9zz53KI&b;XXu52pIbW`$6$5L|tuE0C@6hT3{>)a5Awv#8}a#*)I5( zLdP#%E(ul$%Z7x#)Y#@JvWxMtBB4xAQ#=?Z4A~<{a8_%f9WMOkvO~48SRzJE#B{$F z1>IZ&87!1bNrs$@XFo|+ia8=8oAHbR(n~dmJn!z`N>pZK{NvC+E7tbw%++Gp$D?PJ zTMKk3?>F+^!o8RTApZmd+{bquxK@A|1Uc{z24X-bKa z48$T}Wo&jDAK0U1lUrjdbLVWrK?)xyu$0jLtO()k@1_*Nga88v1E*JShhD{P9yb_V zOXn5I1F6W|=GEdUI|BveCCMqdPcSHMun8j>GGxv<07iKToebztfHM!{Bu|XQig+QA zA`$oeq|MBsQxx&6eXALPPh0F=hChN~jXbsJEWk!&&Y)((49-%i3sc_CD1qSP~l$ z2CK}GJOE?x@b+6xoVDCstNij&ohvpm?Ux`9SB)h$R`vTXP84qHU*3QY<^cX*|8KgS z9iEfnqj&IlKexAo>3x6VPy0VR(*JJr&+k_8J{det&t$%GUZRbqFZm>1!< zuA`op=MmP_ao-$XqmYczkx;Lj&neC|yQ8D-AS>2m(^S5r&%i~;(S|jgd43^FN}RLN zqy`VvQU4{hj32(``>x+qC?DkgS+KebAo}sY3_QYr@?UUit^c}=7Fz+UZheN z9$L4?0)zLnq?-9f1xhHnqlzsKV%X({{PuiWZ$MyDp}Hxk!SO~I?~9Q8e#~da*($Cm z!Ktr#O4!RQ-Xk%aGhWvFZ0@pt<(aFB+a~0G!Qd6U zMWdT-y_+>H3$jH1-~Ig`+I|?weFF-(R3qyM;4$_!V6ReVga*YUMwki=YsS(%N;L5fO3eva%`Tv6%7 ziciblglIx3CJrhqV>7dwfJ17x_gtFlf@rd#ckygthgr8xwyeg0pbcaI>G0a{JDh?q zGSxGxmGMyWOWV8|BjYNZG;NY}dAC*JjZH1!kC)+y4K4jBZ$kX#v?m#{#r_PuTfRF*v*6S%mKg$P-Irw`P=@%i-p81^ANh(nXq=3)0;_Ek3U{H-N2 z;P?Fuey63+ry&uu&vV{(jN>B?B*Hn`^BP;8>XYI*LDsyOgBnjNKWjbcbXS$%%m6~X zYf=;cFd2Ej3tKj-gs4ItWF>e~_ey#>d75eT`rOOfI^PL4L3`fxV1Rl80E>Q<70dt? z06`z6SC?0#bHD3Vq1c4)Uk6&(b>CiafbV!t;d_rkN15L;qwc<@4()-+YzL8I14dmU z1YYwn>%{Jn2vwz#-#Pv$o3xS0j=QJ(r3j6sDrD6AK2a#5^M~+*CQv zR_Y>^)sBETq~zaHNvt~CK>svEz5&axJa0TJWO>aj7}nF^;YypgCdjIdS%J2Yf1(qu zQWLZ64_kXXKoIVrOx?c2DB9ZE^>our6ae4|H&_@3PCMK@_#b@|{4t^WA8cW;BqN9& zsh>9CG;x%-^3CLLaWf$(w2=*kY=U|vNupcB!&i2ZwU`;PW$38|_gQyE`jmJ&XFp!D zIEKs9;A-e86XnS=ZDZA^o_Y9KF_Fma;)tq?4u=v!Lu+prBQMN=k{G;ih9U?0e!tKB zlLKSQ@w`p02WQ8wA(S|K+?KYpEDO}77EW8sdR|Kut?Cn z5n^##q|$hOVJI=NnM_s-(rEA$cGYEWP=^8YQbp>43_Wk1%%$$6miRMS^6zDuK672c zI#mU&l6$49lqPwa1GX_C4(W1n(;N#bd$JU(v)UZ4F}tn!bCw|*wWcBgQef%Fuo;c> z&Y`Q50Vxrg_B;lZW!5VlpWp_~B1gm^2cs~j-|*t$R;S0Tv|CTZtmEPz5&+!`Ii#i< zjOfY~$*$xLp9y8o*)}!Fs>c;L`kx<4kb7OUhy`uyYG{%VFSQe~LkA}=+Wn${)p~Aj zcu-5QUwQ07Q81gc42d^-EYVPSrFP2$m_k|+z)@1+NQRATgFU>7>{iYa%l5rDGyrONB|;GetJcfH{A zcKO}Ct~E&9fV+pO-0zoLsdHew{@))s8WaUONieEJaQ`>6#=s;XoL@J%=+*RS^K<lBG*~LIHFanjOOGzFc)f1dfuaa|%&_BXLDuzr67Y0L3A3NFj(ofFOW+*24tu2;oua p7kYmSc`{tCN{d+NyO+`bJJH~Sezr3w<##{-7ji{7P>?1(`fw9vtMdQ= literal 0 HcmV?d00001 diff --git a/data/EcoZUR_meta.rda b/data/EcoZUR_meta.rda new file mode 100644 index 0000000000000000000000000000000000000000..950deec7635e39130a319c055348d6c286588603 GIT binary patch literal 2033 zcmVF?8m zxdiAJbERMi-v|`F#DX@2000m$FGNWur{bsTo~F_~)Y_S-+NY`IKT2fDlVX^jqbI4R z)YJ7(Q`9u^Pg4oCHq{T+Jx%IO>O<6d8bP%K15-^hY*W)nJw|O%YIy*pnvYUx0LGew z(g4w*G{R&7qaXkSA?TV78fl?~4FJ)IF*Gs;jF|*e$xs78XwYe(4H`7aVgn;U84Umi zhME8}000000009(0V(>bs6nG9fB*mh0000041j1F41gLMXc`&-00000Xdn$VG|`Y6 zGz^(EX^;#NlR>6T5uu2|G|&u2fYGKQrkOCBXc{n07)B;c5hO`8o|0#(89bv;C70Y*LFr zgi5GMYg!hnc9i9o&J+~HRg6^1qFb`P<>@6=8V*S!kYL_Dx}%Y-v2hLrf^B<7YN$T`8&@76wTxv3KAa7k#$%j^p8 zKXQ#@_ZTIhW|HWc^jcW3k(VMm)@Vnjj1L4xM1<{j*U-lzXFAJ67K7?>`R(x|dv1KM z%u3m%I?Jmmmw7h}^=X?*dTtZPie2Z|jhK_iEOnfu>EtZzC<+g=g#Y2i;FYNb369-s zPVFC`lZmkfbu9@__e}4;zpbUQwPt?Y*{nCsF9tPz5*O>{NsOnxG2`@q2X-oHU7U5D z1n(7o00IR@00Cw3bhcl=>GAUWtlt~G+1XTgIVcs;@%G<+)dO(>7wg(^GqT}a49EhZ z{bn~b2qxEO?xC!R$EI^!A6?{Uy?Z@QZfk{j4_B)5Qwyz473?A zN=gt6P{_fC0(3w+B!Wl)lKX{hw=*+r`%dl7W@_2EmEA-;)S9cxsmHXhG9wv4$1|An z@ez*h)a_2`1rKyiW!;m9B|&>-Y($%MS)QsB!P%3pj8g1_0z0e@&l0MKtDeDsl z@tKI4Vc*$MbVxSnrWY<%DNt%!7LF>QBtTJ^$_q-h42C&S$P)IFf<;=5MRTHwiWiEa zfLLlLVAZ zr1w=b1d&9fj;FA`KrbQM}c0sq2WdvxFjLhn) zr6iIjnE+~JmV+rGQ3QsFhb9ZacOlu#L_rB4UJJP{QAESSLLv|&hPf4-IQaYbHx4^B zt`65y*v$lFMB7I2H}SND1Ek0z31=?B;l<2Blcf#`aNwO548So3L+Lt}4FnUYcu*fu zNg+fha(K}gXMWrGx5p7UlE*hjq4Jrb?VBqi{K-Vw;@$~oaErJlg91@e ze4Dt!fgp-Vf<}=bluFpWr-JXZGgCH3Te`cwpB`^6GPK2xkpz-R+f4D3l&5IhYjE{d zcE;|fHm1OdNF$_yD4?{8Nfy)in{J9!_^6rPH5)>?-IgRfyVcmBk|crxNhA>jk`9tY zt=y2&3Nf2z$AbfQ)wyOQ?A3Jm`3pg5BC4v1Kun32Ina=TiJmGpwtTt`#8l8k2HMY4 z3q)4{S4Q95G2v1sY@%9Hpkh=Bdr%_S6+7?n|8@9-i|wm%12l%Kh@jKGt23Px6#6y> zp4VQ0Cdb3YlR4G!)vdm#A0LJqtyRPm1VKnD9|r4Ws3r~d4r~sUP))ULFIv3oXk)B2 z5dqrc6cn+bXpoGOMle(u!yra6k^*F=B_=?mj3gqMl(7U!DIk=^Ato45nG+d=Q6$8K zW&k53s3|anVlYz~5J&`4OcsYc8M*0pn=-7?==8Z|S*FiK2D6Lk64> literal 0 HcmV?d00001 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 7eb24e7..fa8bb47 100644 --- a/man/est_nuis.Rd +++ b/man/est_nuis.Rd @@ -31,7 +31,7 @@ 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 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")`.} @@ -46,12 +46,20 @@ est_nuis( \item{verbose}{Do you want to receive updates as this function runs? Default is `TRUE`.} } \value{ -A list of estimated nuisance functions and some diagnostic information. +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 ed9e5ba..131d33b 100644 --- a/man/ndFit.Rd +++ b/man/ndFit.Rd @@ -40,7 +40,7 @@ ndFit( \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`.} @@ -63,13 +63,24 @@ ndFit( \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/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..db681b6 --- /dev/null +++ b/vignettes/intro_niceday.Rmd @@ -0,0 +1,120 @@ +--- +title: "Introduction to niceday" +output: rmarkdown::html_vignette +author: "Grant Hopkins, Sarah Teichman, and Amy Willis" +date: "`r Sys.Date()`" +vignette: > + %\VignetteIndexEntry{Introduction to 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} +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 covariate +- "preferential sampling": we assume that observed data are proportional to absolute data that we wish we could observe, with potential unknown sample-specific and category-specific scaling factors + +In niceday, there are two types of parameters that we estimate, and they depend on the assumptions that we are willing to make about our data. + +## 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"), # choosing single learner for the example to run quickly, + sl.lib.m = c("SL.mean"), # 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 count table, with observations on the rows and categories on the columns +- `data`: this is our metadata (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 + +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} +ndRes$coef$Psi2 %>% head() +ndRes$coef$Psi2g %>% head() +``` + +Finally, let's plot our estimates and 95% confidence intervals. + +```{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 preprint: + +Hopkins et al. "Nonparametric Identification and Estimation of Ratios of Multi-Category Means under Preferential Sampling." (2025+) From e432f2a895c58bb98895a57a7ad802e92f76d460 Mon Sep 17 00:00:00 2001 From: Sarah Teichman Date: Thu, 9 Oct 2025 11:26:41 -0700 Subject: [PATCH 4/5] add assumptions to vignette --- R/print.ndFit.R | 14 +++++++++++++- vignettes/intro_niceday.Rmd | 38 +++++++++++++++++++++++++++++++------ 2 files changed, 45 insertions(+), 7 deletions(-) diff --git a/R/print.ndFit.R b/R/print.ndFit.R index c946c95..6a4e2ee 100644 --- a/R/print.ndFit.R +++ b/R/print.ndFit.R @@ -16,7 +16,19 @@ print.ndFit <- function(x, cat("\nCall:\n", paste(deparse(x$call), sep = "\n", collapse = "\n"), "\n\n", sep = "") - cat("\nCoefficient estimates with the largest magnitudes:\n") + 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] diff --git a/vignettes/intro_niceday.Rmd b/vignettes/intro_niceday.Rmd index db681b6..33936b4 100644 --- a/vignettes/intro_niceday.Rmd +++ b/vignettes/intro_niceday.Rmd @@ -4,7 +4,7 @@ output: rmarkdown::html_vignette author: "Grant Hopkins, Sarah Teichman, and Amy Willis" date: "`r Sys.Date()`" vignette: > - %\VignetteIndexEntry{Introduction to niceday} + %\VignetteIndexEntry{intro_niceday} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- @@ -27,7 +27,7 @@ First we will install `niceday` if haven't already. Next we will load `niceday`, as well as `tidyverse` suite packages that we will use. -```{r setup} +```{r setup, message = FALSE} library(niceday) library(dplyr) library(ggplot2) @@ -38,10 +38,26 @@ library(ggplot2) 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 covariate +- "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 data that we wish we could observe, with potential unknown sample-specific and category-specific scaling factors -In niceday, there are two types of parameters that we estimate, and they depend on the assumptions that we are willing to make about our data. +In `niceday`, we identify our parameters with 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 @@ -90,12 +106,22 @@ There are more optional arguments for `ndFit()` that we won't cover here. Feel f 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() ``` -Finally, let's plot our estimates and 95% confidence intervals. +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 making a different independence assumption and targeting a different parameter. + +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 %>% @@ -115,6 +141,6 @@ Here we can see that our 95% confidence intervals range from approximately $-4.5 ## Citation -If you use `niceday` please cite the preprint: +If you use `niceday` please cite the paper: Hopkins et al. "Nonparametric Identification and Estimation of Ratios of Multi-Category Means under Preferential Sampling." (2025+) From 38c9b172762fb29782eab15bc3bf869fbde80945 Mon Sep 17 00:00:00 2001 From: Amy Willis Date: Fri, 31 Oct 2025 13:43:39 -0700 Subject: [PATCH 5/5] update readme and vignette --- README.Rmd | 24 ++++--- README.md | 123 +++++++++++---------------------- vignettes/intro_niceday.Rmd | 27 +++++--- vignettes/niceday_vignette.Rmd | 36 ---------- 4 files changed, 73 insertions(+), 137 deletions(-) delete mode 100644 vignettes/niceday_vignette.Rmd diff --git a/README.Rmd b/README.Rmd index a072c47..63d6194 100644 --- a/README.Rmd +++ b/README.Rmd @@ -19,28 +19,34 @@ knitr::opts_chunk$set( [![R-CMD-check](https://github.com/statdivlab/niceday/workflows/R-CMD-check/badge.svg)](https://github.com/statdivlab/niceday/actions) -`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 count matrix `W`, metadata `data`, the treatment `A`, and covariates to adjust for `X`. +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) - data(EcoZUR_meta) data(EcoZUR_count) -ndFit(W = EcoZUR_count[, 1:50], # consider only the first 50 taxa to run quickly +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, @@ -50,11 +56,13 @@ ndFit(W = EcoZUR_count[, 1:50], # consider only the first 50 taxa 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 6d5b07c..01e0625 100644 --- a/README.md +++ b/README.md @@ -8,54 +8,56 @@ [![R-CMD-check](https://github.com/statdivlab/niceday/workflows/R-CMD-check/badge.svg)](https://github.com/statdivlab/niceday/actions) -`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 count matrix `W`, metadata `data`, the treatment `A`, and -covariates to adjust for `X`. +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) -#> Loading required package: gam -#> Loading required package: splines -#> Loading required package: foreach -#> Loaded gam 1.22-6 -#> Loading required package: MASS -#> Loading required package: glmnet -#> Loading required package: Matrix -#> Loaded glmnet 4.1-10 -#> Loading required package: hal9001 -#> Loading required package: Rcpp -#> hal9001 v0.4.6: The Scalable Highly Adaptive Lasso -#> note: fit_hal defaults have changed. See ?fit_hal for details -#> Loading required package: ranger -#> Loading required package: xgboost -#> Loading required package: withr -#> Loading required package: SuperLearner -#> Loading required package: nnls -#> Super Learner -#> Version: 2.0-29 -#> Package created on 2024-02-06 - data(EcoZUR_meta) data(EcoZUR_count) -ndFit(W = EcoZUR_count[, 1:50], # consider only the first 50 taxa to run quickly +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, @@ -63,66 +65,19 @@ ndFit(W = EcoZUR_count[, 1:50], # consider only the first 50 taxa to run quickly 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 -#> -#> Call: -#> ndFit(W = EcoZUR_count[, 1:50], A = ~Diarrhea, X = ~sex + age_months, -#> data = EcoZUR_meta, num_crossval_folds = 2, num_crossfit_folds = 2, -#> sl.lib.pi = c("SL.mean"), sl.lib.m = c("SL.mean")) -#> -#> -#> Coefficient estimates with the largest magnitudes: -#> category est se lower_marg upper_marg lower_sim upper_sim -#> 13 13 -3.2741166 0.6118426 -4.47330614 -2.07492704 -5.2715008 -1.2767324 -#> 3 3 -2.6157685 0.7393277 -4.06482416 -1.16671283 -5.0293326 -0.2022044 -#> 22 22 2.1413861 0.9688809 0.24241444 4.04035780 -1.0215634 5.3043356 -#> 39 39 2.0849190 0.7705275 0.57471289 3.59512520 -0.4304980 4.6003361 -#> 30 30 -1.8174642 0.5879262 -2.96977845 -0.66515004 -3.7367723 0.1018438 -#> 14 14 1.7374620 0.4585679 0.83868546 2.63623844 0.2404493 3.2344746 -#> 21 21 1.6905260 0.5946739 0.52498671 2.85606538 -0.2508099 3.6318620 -#> 20 20 1.5157064 0.5558763 0.42620878 2.60520394 -0.2989736 3.3303863 -#> 33 33 -1.4345277 0.5776259 -2.56665373 -0.30240164 -3.3202101 0.4511547 -#> 32 32 -1.3580239 0.6794349 -2.68969173 -0.02635606 -3.5760655 0.8600177 -#> 15 15 -1.3004635 0.6041093 -2.48449595 -0.11643105 -3.2726018 0.6716748 -#> 31 31 1.2779538 0.6321374 0.03898736 2.51692033 -0.7856833 3.3415909 -#> 45 45 -1.2336750 0.4405577 -2.09715223 -0.37019774 -2.6718928 0.2045428 -#> 27 27 1.1777250 0.6579759 -0.11188403 2.46733410 -0.9702630 3.3257130 -#> 42 42 -1.0418668 0.5415186 -2.10322383 0.01949026 -2.8096755 0.7259420 -#> 12 12 1.0119062 0.8021781 -0.56033395 2.58414645 -1.6068355 3.6306480 -#> 25 25 1.0117772 0.3666362 0.29318338 1.73037094 -0.1851211 2.2086754 -#> 23 23 1.0073781 0.5866781 -0.14248979 2.15724602 -0.9078554 2.9226116 -#> 8 8 -0.8832539 0.4902861 -1.84419704 0.07768917 -2.4838120 0.7173042 -#> 4 4 -0.8691163 0.3683387 -1.59104691 -0.14718560 -2.0715724 0.3333399 -#> type -#> 13 TMLE -#> 3 TMLE -#> 22 TMLE -#> 39 TMLE -#> 30 TMLE -#> 14 TMLE -#> 21 TMLE -#> 20 TMLE -#> 33 TMLE -#> 32 TMLE -#> 15 TMLE -#> 31 TMLE -#> 45 TMLE -#> 27 TMLE -#> 42 TMLE -#> 12 TMLE -#> 25 TMLE -#> 23 TMLE -#> 8 TMLE -#> 4 TMLE -#> 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. ``` +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.” +Means under Preferential Sampling.” ## Bug reports and feature requests diff --git a/vignettes/intro_niceday.Rmd b/vignettes/intro_niceday.Rmd index 33936b4..d54afd8 100644 --- a/vignettes/intro_niceday.Rmd +++ b/vignettes/intro_niceday.Rmd @@ -39,9 +39,9 @@ This vignette provides an introduction to `niceday`, a statistical package that - "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 data that we wish we could observe, with potential unknown sample-specific and category-specific scaling factors +- "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 with assumptions about our observed data and parameters. We assume the following: +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. @@ -86,19 +86,19 @@ ndRes <- ndFit(W = EcoZUR_count, 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 + 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 count table, with observations on the rows and categories on the columns -- `data`: this is our metadata (which includes all covariates in `A` and `X` formulas) +- `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 +- `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! @@ -119,7 +119,16 @@ 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 making a different independence assumption and targeting a different parameter. +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`. @@ -143,4 +152,4 @@ Here we can see that our 95% confidence intervals range from approximately $-4.5 If you use `niceday` please cite the paper: -Hopkins et al. "Nonparametric Identification and Estimation of Ratios of Multi-Category Means under Preferential Sampling." (2025+) +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!]