From da26ad6546a6fff228017af4bd9aed368da5477d Mon Sep 17 00:00:00 2001 From: Brantly Callaway Date: Fri, 10 Jan 2025 11:12:17 -0500 Subject: [PATCH] update documentation --- R/classes.R | 224 ++++++++++++++------------- R/pte.R | 4 + R/zzz.R | 2 +- man/group_time_att.Rd | 10 +- man/pte2.Rd | 4 + man/pte_emp_boot.Rd | 9 +- tests/testthat/test-did-inference.R | 45 +++--- tests/testthat/testthat-problems.rds | Bin 0 -> 14769 bytes 8 files changed, 154 insertions(+), 144 deletions(-) create mode 100644 tests/testthat/testthat-problems.rds diff --git a/R/classes.R b/R/classes.R index 4b6704c..41a90b2 100644 --- a/R/classes.R +++ b/R/classes.R @@ -16,12 +16,11 @@ #' assumption #' @param Wpval p-value for Wald pre-test of ATT(g,t) version of parallel #' trends assumption +#' @param cband logical indicating whether or not to report a confidence band #' @param alp significance level #' @param ptep \code{pte_params} object -#' -#' @inheritParams pte_params -#' @inheritParams attgt_if -#' +#' @param extra_gt_returns list containing extra returns at the group-time level +#' #' @return object of class \code{group_time_att} #' #' @export @@ -39,21 +38,22 @@ group_time_att <- function(group, alp, ptep, extra_gt_returns) { - - out <- list(group=group, - t=time.period, - att=att, - V_analytical=V_analytical, - se=se, - crit_val=crit_val, - inf_func=inf_func, - n=n, - W=W, - Wpval=Wpval, - cband=cband, - alp=alp, - ptep=ptep, - extra_gt_returns=extra_gt_returns) + out <- list( + group = group, + t = time.period, + att = att, + V_analytical = V_analytical, + se = se, + crit_val = crit_val, + inf_func = inf_func, + n = n, + W = W, + Wpval = Wpval, + cband = cband, + alp = alp, + ptep = ptep, + extra_gt_returns = extra_gt_returns + ) # might eventually drop this, but add a few things to get access # to aggregations from did package @@ -64,7 +64,7 @@ group_time_att <- function(group, out$inffunc <- out$inf_func - class(out) <- c("group_time_att","MP") + class(out) <- c("group_time_att", "MP") out } @@ -83,26 +83,27 @@ summary.group_time_att <- function(object, ...) { # group time average treatment effects cat("Group-Time Average Treatment Effects:\n") - cband_text1a <- paste0(100*(1-group_time_att_obj$alp),"% ") + cband_text1a <- paste0(100 * (1 - group_time_att_obj$alp), "% ") cband_text1b <- ifelse(group_time_att_obj$DIDparams$bstrap, - ifelse(group_time_att_obj$DIDparams$cband, "Simult. ", "Pointwise "), - "Pointwise ") + ifelse(group_time_att_obj$DIDparams$cband, "Simult. ", "Pointwise "), + "Pointwise " + ) cband_text1 <- paste0("[", cband_text1a, cband_text1b) - cband_lower <- group_time_att_obj$att - group_time_att_obj$crit_val*group_time_att_obj$se - cband_upper <- group_time_att_obj$att + group_time_att_obj$crit_val*group_time_att_obj$se + cband_lower <- group_time_att_obj$att - group_time_att_obj$crit_val * group_time_att_obj$se + cband_upper <- group_time_att_obj$att + group_time_att_obj$crit_val * group_time_att_obj$se sig <- (cband_upper < 0) | (cband_lower > 0) sig[is.na(sig)] <- FALSE sig_text <- ifelse(sig, "*", "") out <- cbind.data.frame(group_time_att_obj$group, group_time_att_obj$t, group_time_att_obj$att, group_time_att_obj$se, cband_lower, cband_upper) - out <- round(out,4) + out <- round(out, 4) out <- cbind.data.frame(out, sig_text) - colnames(out) <- c("Group", "Time", "ATT(g,t)","Std. Error", cband_text1, "Conf. Band]", "") - print(out, row.names=FALSE) + colnames(out) <- c("Group", "Time", "ATT(g,t)", "Std. Error", cband_text1, "Conf. Band]", "") + print(out, row.names = FALSE) cat("---\n") cat("Signif. codes: `*' confidence band does not cover 0") cat("\n\n") @@ -113,7 +114,6 @@ summary.group_time_att <- function(object, ...) { cat(as.character(group_time_att_obj$Wpval)) cat("\n") } - } #' @title print.group_time_att @@ -124,8 +124,8 @@ summary.group_time_att <- function(object, ...) { #' @param ... extra arguments #' #' @export -print.group_time_att <- function(x,...) { - summary.group_time_att(x,...) +print.group_time_att <- function(x, ...) { + summary.group_time_att(x, ...) } @@ -143,11 +143,12 @@ pte_results <- function(att_gt, overall_att, event_study, ptep) { - - out <- list(att_gt=att_gt, - overall_att=overall_att, - event_study=event_study, - ptep=ptep) + out <- list( + att_gt = att_gt, + overall_att = overall_att, + event_study = event_study, + ptep = ptep + ) class(out) <- "pte_results" @@ -164,39 +165,39 @@ pte_results <- function(att_gt, #' #' @export summary.pte_results <- function(object, ...) { - overall_att <- object$overall_att$overall.att overall_se <- object$overall_att$overall.se event_study <- object$event_study event_study_att_egt <- object$event_study$att.egt event_study_se <- object$event_study$se.egt - + # overall estimates alp <- object$ptep$alp - pointwise_cval <- qnorm(1-alp/2) - overall_cband_upper <- overall_att + pointwise_cval*overall_se - overall_cband_lower <- overall_att - pointwise_cval*overall_se + pointwise_cval <- qnorm(1 - alp / 2) + overall_cband_upper <- overall_att + pointwise_cval * overall_se + overall_cband_lower <- overall_att - pointwise_cval * overall_se out1 <- cbind.data.frame(overall_att, overall_se, overall_cband_lower, overall_cband_upper) out1 <- round(out1, 4) overall_sig <- (overall_cband_upper < 0) | (overall_cband_lower > 0) overall_sig[is.na(overall_sig)] <- FALSE overall_sig_text <- ifelse(overall_sig, "*", "") out1 <- cbind.data.frame(out1, overall_sig_text) - + # Event study # header bstrap <- object$DIDparams$bstrap cband <- object$DIDparams$cband - cband_text1a <- paste0(100*(1-alp),"% ") + cband_text1a <- paste0(100 * (1 - alp), "% ") cband_text1b <- ifelse(bstrap, - ifelse(cband, "Simult. ", "Pointwise "), - "Pointwise ") + ifelse(cband, "Simult. ", "Pointwise "), + "Pointwise " + ) cband_text1 <- paste0("[", cband_text1a, cband_text1b) - cband_lower <- event_study$att.egt - event_study$crit.val.egt*event_study$se.egt - cband_upper <- event_study$att.egt + event_study$crit.val.egt*event_study$se.egt + cband_lower <- event_study$att.egt - event_study$crit.val.egt * event_study$se.egt + cband_upper <- event_study$att.egt + event_study$crit.val.egt * event_study$se.egt sig <- (cband_upper < 0) | (cband_lower > 0) sig[is.na(sig)] <- FALSE @@ -207,13 +208,15 @@ summary.pte_results <- function(object, ...) { out2 <- cbind.data.frame(out2, sig_text) c1name <- "Event Time" - colnames(out2) <- c(c1name, "Estimate","Std. Error", cband_text1, "Conf. Band]", "") - - out <- list(overall_att=out1, - event_study=out2, - alp=alp, - bstrap=bstrap, - cband=cband) + colnames(out2) <- c(c1name, "Estimate", "Std. Error", cband_text1, "Conf. Band]", "") + + out <- list( + overall_att = out1, + event_study = out2, + alp = alp, + bstrap = bstrap, + cband = cband + ) class(out) <- "summary.pte_results" @@ -228,8 +231,8 @@ summary.pte_results <- function(object, ...) { #' @param ... extra arguments #' #' @export -print.pte_results <- function(x,...) { - #summary.pte_results(x,...) +print.pte_results <- function(x, ...) { + # summary.pte_results(x,...) NextMethod(x) invisible(x) } @@ -244,25 +247,26 @@ print.pte_results <- function(x,...) { #' @param ... extra arguments #' #' @export -print.summary.pte_results <- function(x,...) { +print.summary.pte_results <- function(x, ...) { object <- x out1 <- object$overall_att out2 <- object$event_study alp <- object$alp - + # overall att cat("\n") cat("Overall ATT: \n") - colnames(out1) <- c("ATT"," Std. Error", paste0(" [ ",100*(1-alp),"% "), "Conf. Int.]","") - print(out1, row.names=FALSE) + colnames(out1) <- c("ATT", " Std. Error", paste0(" [ ", 100 * (1 - alp), "% "), "Conf. Int.]", "") + print(out1, row.names = FALSE) cat("\n\n") # event study - c1name <- "Event time"; cat("Dynamic Effects:") + c1name <- "Event time" + cat("Dynamic Effects:") cat("\n") - print(out2, row.names=FALSE, justify = "centre") - + print(out2, row.names = FALSE, justify = "centre") + cat("---\n") cat("Signif. codes: `*' confidence band does not cover 0") cat("\n\n") @@ -285,9 +289,8 @@ print.summary.pte_results <- function(x,...) { #' @return attgt_if object #' #' @export -attgt_if <- function(attgt, inf_func, extra_gt_returns=NULL) { - - out <- list(attgt=attgt, inf_func=inf_func, extra_gt_returns=extra_gt_returns) +attgt_if <- function(attgt, inf_func, extra_gt_returns = NULL) { + out <- list(attgt = attgt, inf_func = inf_func, extra_gt_returns = extra_gt_returns) class(out) <- "attgt_if" out } @@ -302,9 +305,8 @@ attgt_if <- function(attgt, inf_func, extra_gt_returns=NULL) { #' @return attgt_noif object #' #' @export -attgt_noif <- function(attgt, extra_gt_returns=NULL) { - - out <- list(attgt=attgt, extra_gt_returns=extra_gt_returns) +attgt_noif <- function(attgt, extra_gt_returns = NULL) { + out <- list(attgt = attgt, extra_gt_returns = extra_gt_returns) class(out) <- "attgt_noif" out } @@ -318,24 +320,23 @@ attgt_noif <- function(attgt, extra_gt_returns=NULL) { #' #' @param data data that will be checked to see if has right format for #' computing group-time average treatment effects -#' +#' #' @return \code{gt_data_frame} object #' #' @export gt_data_frame <- function(data) { - cnames <- colnames(data) - if( !all(c("G","id","period","name","Y","D") %in% cnames) ) { + if (!all(c("G", "id", "period", "name", "Y", "D") %in% cnames)) { warning("group-time subset of data does not contain correct column names") } - if( !("gt_data_frame" %in% class(data)) ) { + if (!("gt_data_frame" %in% class(data))) { class(data) <- c("gt_data_frame", class(data)) } data } - + #' @title pte_emp_boot #' #' @title Class for holding \code{pte} empirical bootstrap results @@ -349,10 +350,10 @@ gt_data_frame <- function(data) { #' @param dyn_weights list containing weights on underlying ATT(g,t) #' for each value of \code{e} corresponding to the dynamic treatment #' effect parameters. -#' @param group_weights vector containing weights on underlying ATT(g,t) -#' to group specific treatment effects -#' -#' @inherit attgt_if +#' @param group_weights list containing weights on underlying ATT(g,t) +#' corresponding to deliver averaged group-specific treatment effects +#' +#' @inheritParams attgt_if #' #' @return pte_emp_boot object #' @@ -361,19 +362,20 @@ pte_emp_boot <- function(attgt_results, overall_results, group_results, dyn_results, - overall_weights=NULL, - dyn_weights=NULL, - group_weights=NULL, - extra_gt_returns=NULL) { - - out <- list(attgt_results=attgt_results, - overall_results=overall_results, - group_results=group_results, - dyn_results=dyn_results, - overall_weights=overall_weights, - dyn_weights=dyn_weights, - group_weights=group_weights, - extra_gt_returns=extra_gt_returns) + overall_weights = NULL, + dyn_weights = NULL, + group_weights = NULL, + extra_gt_returns = NULL) { + out <- list( + attgt_results = attgt_results, + overall_results = overall_results, + group_results = group_results, + dyn_results = dyn_results, + overall_weights = overall_weights, + dyn_weights = dyn_weights, + group_weights = group_weights, + extra_gt_returns = extra_gt_returns + ) class(out) <- "pte_emp_boot" @@ -391,40 +393,40 @@ pte_emp_boot <- function(attgt_results, #' #' @export summary.pte_emp_boot <- function(object, ...) { - overall_att <- object$overall_results$att overall_se <- object$overall_results$se - #event_study <- object$event_study + # event_study <- object$event_study event_study_att <- object$dyn_results$att.e event_study_se <- object$dyn_results$se event_study_e <- object$dyn_results$e - + # overall estimates alp <- .05 # hard coded - pointwise_cval <- qnorm(1-alp/2) - overall_cband_upper <- overall_att + pointwise_cval*overall_se - overall_cband_lower <- overall_att - pointwise_cval*overall_se + pointwise_cval <- qnorm(1 - alp / 2) + overall_cband_upper <- overall_att + pointwise_cval * overall_se + overall_cband_lower <- overall_att - pointwise_cval * overall_se out1 <- cbind.data.frame(overall_att, overall_se, overall_cband_lower, overall_cband_upper) out1 <- round(out1, 4) overall_sig <- (overall_cband_upper < 0) | (overall_cband_lower > 0) overall_sig[is.na(overall_sig)] <- FALSE overall_sig_text <- ifelse(overall_sig, "*", "") out1 <- cbind.data.frame(out1, overall_sig_text) - + # Event study # header bstrap <- TRUE cband <- FALSE - cband_text1a <- paste0(100*(1-alp),"% ") + cband_text1a <- paste0(100 * (1 - alp), "% ") cband_text1b <- ifelse(bstrap, - ifelse(cband, "Simult. ", "Pointwise "), - "Pointwise ") + ifelse(cband, "Simult. ", "Pointwise "), + "Pointwise " + ) cband_text1 <- paste0("[", cband_text1a, cband_text1b) - cband_lower <- event_study_att - pointwise_cval*event_study_se - cband_upper <- event_study_att + pointwise_cval*event_study_se + cband_lower <- event_study_att - pointwise_cval * event_study_se + cband_upper <- event_study_att + pointwise_cval * event_study_se sig <- (cband_upper < 0) | (cband_lower > 0) sig[is.na(sig)] <- FALSE @@ -435,13 +437,15 @@ summary.pte_emp_boot <- function(object, ...) { out2 <- cbind.data.frame(out2, sig_text) c1name <- "Event Time" - colnames(out2) <- c(c1name, "Estimate","Std. Error", cband_text1, "Conf. Band]", "") - - out <- list(overall_att=out1, - event_study=out2, - alp=alp, - bstrap=bstrap, - cband=cband) + colnames(out2) <- c(c1name, "Estimate", "Std. Error", cband_text1, "Conf. Band]", "") + + out <- list( + overall_att = out1, + event_study = out2, + alp = alp, + bstrap = bstrap, + cband = cband + ) class(out) <- "summary.pte_results" diff --git a/R/pte.R b/R/pte.R index 8fc37d9..57ce0bd 100644 --- a/R/pte.R +++ b/R/pte.R @@ -514,6 +514,10 @@ pte <- function(yname, #' package will use the empirical bootstrap no matter what the value of this #' parameter. #' +#' @param process_dtt_gt An optional function to customize results when +#' the gt-specific function returns the distribution of treated and untreated +#' potential outcomes. +#' #' @param ... extra arguments that can be passed to create the correct subsets #' of the data (depending on \code{subset_fun}), to estimate group time #' average treatment effects (depending on \code{attgt_fun}), or to diff --git a/R/zzz.R b/R/zzz.R index 3f83bb9..46bb7ef 100644 --- a/R/zzz.R +++ b/R/zzz.R @@ -1 +1 @@ -utils::globalVariables(c("e","post","ciu","cil","G","Y","att","id","name","period","att.e","att.g","group","length.e","length.group","time.period")) +utils::globalVariables(c("e", "post", "ciu", "cil", "G", "Y", "att", "id", "name", "period", "att.e", "att.g", "group", "length.e", "length.group", "time.period", ".treat")) diff --git a/man/group_time_att.Rd b/man/group_time_att.Rd index 8ef8406..c71f42e 100644 --- a/man/group_time_att.Rd +++ b/man/group_time_att.Rd @@ -46,19 +46,13 @@ assumption} \item{Wpval}{p-value for Wald pre-test of ATT(g,t) version of parallel trends assumption} -\item{cband}{whether or not to report a uniform (instead of pointwise) -confidence band (default is TRUE)} +\item{cband}{logical indicating whether or not to report a confidence band} \item{alp}{significance level} \item{ptep}{\code{pte_params} object} -\item{extra_gt_returns}{A place to return anything extra from particular -group-time average treatment effect calculations. For DID, this might -be something like propensity score estimates, regressions of untreated -potential outcomes on covariates. For ife, this could be something -like the first step regression 2sls estimates. This argument is also -potentially useful for debugging.} +\item{extra_gt_returns}{list containing extra returns at the group-time level} } \value{ object of class \code{group_time_att} diff --git a/man/pte2.Rd b/man/pte2.Rd index 45ae172..8e609aa 100644 --- a/man/pte2.Rd +++ b/man/pte2.Rd @@ -124,6 +124,10 @@ their own choices for \code{gt_type}} this is a specific quantile at which to report results, e.g., \code{ret_quantile = 0.5} will return that the qte at the median.} +\item{process_dtt_gt}{An optional function to customize results when +the gt-specific function returns the distribution of treated and untreated +potential outcomes.} + \item{global_fun}{Logical indicating whether or not untreated potential outcomes can be estimated in one shot, i.e., for all groups and time periods. Main use case would be for one-shot imputation estimators. Not supported yet.} diff --git a/man/pte_emp_boot.Rd b/man/pte_emp_boot.Rd index 10bb15f..f080e63 100644 --- a/man/pte_emp_boot.Rd +++ b/man/pte_emp_boot.Rd @@ -31,8 +31,8 @@ for overall treatment effect parameter} for each value of \code{e} corresponding to the dynamic treatment effect parameters.} -\item{group_weights}{vector containing weights on underlying ATT(g,t) -to group specific treatment effects} +\item{group_weights}{list containing weights on underlying ATT(g,t) +corresponding to deliver averaged group-specific treatment effects} \item{extra_gt_returns}{A place to return anything extra from particular group-time average treatment effect calculations. For DID, this might @@ -45,6 +45,7 @@ potentially useful for debugging.} pte_emp_boot object } \description{ -Class for holding group-time average treatment effects -along with their influence function +pte_emp_boot + +Class for holding \code{pte} empirical bootstrap results } diff --git a/tests/testthat/test-did-inference.R b/tests/testthat/test-did-inference.R index 2c28854..f168936 100644 --- a/tests/testthat/test-did-inference.R +++ b/tests/testthat/test-did-inference.R @@ -6,33 +6,36 @@ library(did) library(pbapply) -skip_inf <- menu(c("run", "skip"), title="would you like to run or skip the inference tests?") -skip_inf <- skip_inf == 2 +# skip_inf <- menu(c("run", "skip"), title="would you like to run or skip the inference tests?") +# skip_inf <- skip_inf == 2 test_that("tests for inference", { - skip_if(skip_inf) - cl <- as.numeric(readline("how many clusters would like to use for the inference tests? ")) + # skip_if(skip_inf) + # cl <- as.numeric(readline("how many clusters would like to use for the inference tests? ")) + cl <- 1 mc_sims <- 100 rejs <- pblapply(1:mc_sims, function(mc) { sp <- did::reset.sim() data <- did::build_sim_dataset(sp) - - res <- pte(yname="Y", - gname="G", - tname="period", - idname="id", - data=data, - setup_pte_fun=setup_pte, - subset_fun=two_by_two_subset, - attgt_fun=did_attgt, - xformla=~X) + + res <- pte( + yname = "Y", + gname = "G", + tname = "period", + idname = "id", + data = data, + setup_pte_fun = setup_pte, + subset_fun = two_by_two_subset, + attgt_fun = did_attgt, + xformla = ~X + ) # truth is that att = 1 - tstat <- (res$overall_att$overall.att-1) / res$overall_att$overall.se - rej <- 1*( abs(tstat) > qnorm(.975)) + tstat <- (res$overall_att$overall.att - 1) / res$overall_att$overall.se + rej <- 1 * (abs(tstat) > qnorm(.975)) rej - }, cl=cl) - + }, cl = cl) + rej_frac <- mean(unlist(rejs)) - - expect_equal(rej_frac, 0.06, tolerance=.05) # make test fail if reject 0 -}) \ No newline at end of file + + expect_equal(rej_frac, 0.06, tolerance = .05) # make test fail if reject 0 +}) diff --git a/tests/testthat/testthat-problems.rds b/tests/testthat/testthat-problems.rds new file mode 100644 index 0000000000000000000000000000000000000000..1c243746061e2a10fe47e96e948be0cc5d44940f GIT binary patch literal 14769 zcmY*=c|26_`@VgTcVrKheJMruWlBZJk|j%+5JEz-Wf&?XTiKT^Wy?C*qA(--E@Y=M z7_!b7+YGb(&h+{I^P4|(UT2;;=XuV#@9Vkl>$($*XQBD`q2CE-hn$-ZJl?*=o$|gx zuY>!-BU@I9sP_y4uiY3}6kMPFHxrp?^0d~zx~s#>=8vP`TW=R>f6*v*=Lr( zm9)`Lldms2|!RUioIcGiaA^9a|3u@-PT_2vhr>gNo73o!~92y{@{6X{q z+NNZ{Cy8x;2X{R~$tcTru;XSKB(csxRrm9_O{Q}v`Iv_{f8EP&X!^WE#m|+_-2tV! zmZZl$+>W`Et_vxZB27vg3Y#4OLQh zE=}pRl2R^x=TULbLR%Rsfw~-Z`MzQB2d%u0KwsGRsTGSfE9P>&`Ol^rd{Wu^b8#Zx1o@*_|Rl|60aTv#Gc{oy(;Rlg%vo@a35lxzN zuAr64H*+8Naos!8YPyfBFJ{3<9YQnXU>QMr@x$DdO6BDd=r z+bZu|rIb~jf5dq3mS7()UtSO4x%F=OE@376=cm-y4Bs3ssm*?HB@n-KSu?6WAsscp z(3>RMzH8MCH19c^!at`>dZM~N8%I2=oUL5`J;k}?!tj!PrPcG$J74r3oKc)k`L6RJ zvBV0@V7-_Xr`x>X+mY3>(x(EAFXm_ zg0mL$7BkClilZA+VjJ=qBt8kF#_cVBfg>5{oc*c%@PY3_WJi+s~S>Dy= z9qd2HeUATVceZK&k+!0tJ;<3Gm$$O^YvFF&&(yM)a@ms=cWSP1sO3#>?U+b6k5>s9 z=HKPul1LhS95VM@Fxf32PTo4WfuF|B_f_AypgDoJ+@_|c5uSt|3k_h&+5isIo>`WZ2Lc5vnXRVbg!( zp_z2x{^qouiOm;nY{ItBfGKeDUfCbxKNXHlrb{7<$T@tvf>GWM5^`x}z##`AD*Kc9 zxcLsqnf*7;Uy!=I-|8JUP3QR2CscWZo*Z;jcG@%PQ*;~A^45Q+v7eE!+GHZ<5vT*g zu$@{k{pD*sIz1uaC&k-(rV_YEOC-!aA(GN!E$|Eb8B_?2bvi<15bTS7@heb{%l#tJ z0RocHRZ2uOy-;{9lm>6q6B;u?_BBUc7~G;MAtAq93HAh*esgJywgmFb0_A>2@ucgF z(1^P9-Hqcfek^!nf!sofQ*lSN0+>UFw8e+zI+P?m{29>}BgQ}~Nsw0$ui$EY&hpeP zRK73`c@^uEuv3Nne3h1HbVCqi5Q(8lBSUXH9zD-XyXtbPQK9z^m{Be1yd^3q5xxpy zgh;ybZlyw9R4<-~=;xl}1oLhPd;Em4vR$pWTg>9{r? zd}zXb8Kh;q4Xxq^&#_U48od$W4Ej0cDWffT#Rlp@f4rjtM>Dv{g zwjpd#_cL;V%YpDf>?bQByT~enBPZ{PVdwFoz-f;tn8#CEif&m4j{=wz_~9Swh78bqwVjN|e$K$9AEJE9inX z!tmqop3N9v#3)f|jdhQ_ZV*Fy@n$tWA%wE;8a*i#Zb*@C7at%VSq2I|4m~4&L@6l8 zasclINW6)HN-TmpLU3tIfd4fK{N#ZDwXXkY9fY;8m}0z?sq_svNeiZLoqwU2@^g`;@+zqo{|U62iuWuZl6 z>{*mSo^Z0I)j@PFWGbS^ewg7ljZaG2!9C5q+nJsRO|Rg~bmKLrG?Jv}OHJ3D*pK8R z!R)Owus1{-Z@iodba z)+mCI(*ZBP?I6zUh8!qROsJUw&vN0I5}Ct8zo-CnOCw~5G>WqIiMr|({;fQ``zTu< zOnk`Fy6Q@sRv^SqU)S6G-}_U8!hQn)Ij4%bj-g$Qns$E}1eWUJ=oG`U{V2TV#}-j& z1#hM_qNVW8ufN78gHgWDPoqy9f9@A&@{=C?@}*D1c}^4~nNCFBBZ6*>fvd}*zcCED zbSE9bw@z~^QPW;?Alz0nR7KhBEh+xUM$Ic^4B2;L6?E~&K*5hDPKTe`=xX|qu`<$R zwE`DC?nL_QtRF@!^sq_E3}>OgSFpKjv@}Q?O5*v;3{Sa#0zy;dcx&84c?*9dIjM=} z)!9aoKIxDjX{7u`^Lp&~);!8L?@`0~7O*kf>So@bf?g+&k=-qd6py%Qz@03oxlzJa z>Q9=Fr5-_*Z=B{9QQ`AH{L3E3M-(Yu&RWRD!Wo2J^hE#%Kn%o7J(A|SiWaB`8H$W) zCoSE;qe4~~itU2m%CC+xoICY5A!HSTD$2!HQynNiLp&X1!wc0N>io9&IyFoIZCr3ONy5*biE0O30Y0G${XHB>iDaj6Z*wS-^^XR;)&+Q<~5b{-xr+ z`Vn(Ge7R~JEWt(h4hntLTrkS2cmiuuT;w9grZ)IEcHj35C0Js3$dEPLOP5Dip>m_E zO_-p259d0iBd3%D=NA+v@=~h1Jdc{GG z60L9n&yU9r8JKw}Y3WFZ{`81{Iq>gtv_0lYUIy2$g7xku!i(|<2@lr_=J;zNrv7z4 z{RbHC5jj~eyD*`_?X}}Md4jT}579CBC}WDSUswnr?ZO>T-{tcCUKK=xIWRl(5kAPI zh-=rr`uZ7AUI-OxVSK8)^Yb#je%yGL?{C(YwNv%V`bcC(O7;lMj3}M&huZAitE7;U z1X&yNH5Gj**XySR1LG$Sx--hl9yFEkwMhmSPjY*{yAK%{XK4&_yf>;yJSXJkIBV}9 z&+kg$S@6GNoLB#)H{i*3S!Ivbp=Vy^AKfDN()cl%DPKT^$)p7{7kThb=C6>P#dg}1&2{jgFg_*{2m zZ!>1Ym=)ZuftyVR-{dEp`#Yv(ZS^-=VpaZq>_2zr}x zzpaN8n%HXUngp>&mp$-;K1NTnPe|KDXxfS;F|oA?S=S`jwN{`Rr_2x?1lz$%hld9n z83iA39R?0}IM3;B;3h7wB)7EYL|Hz+SO77Zea9edN|CjVJ50==m`nI@2)Pw@^OSSL zIrcZ3LY7~6E#Ljz?V13VqIX>HuYUKLYu0o^KKX>-^254(Jlj>(wEcmPri|%J1slWM zizgpPMs2fn>vJnElr2Io%rh=(wbA)12T3*Qu``WX8IFD!J%q$--2>B{pb@{xO23*y z9*G_(dVyT7MA-bze4|F+EiN1$^VjArIvU2)<|eSYonPmYTe>If>d&b8pz84GqhTkB zBo+}AE4IgbX&^H;ofsajt$f9p6~FqY(mKTD^XJDr;z0@j!D}2#?`6JxM#ommE=?06 zc44shUgNif270;SrPUYXO?s%;eH}mgBK}yiITXf5Pbv3b1}UwIX+E?SnGI@F{dng? z_jj)I0!3!V1%cN~B)p{Z#{JW2BUZ;~5#FN>9BdEP1$WJajX2g-7D7FICQtPu;!xTS zFfH%O33cXWg{qg}5jkfE_t^S|&2Y(*E8H}er&$Lw zGsgFA-m$BX==nMfC+htSSy{c#UC*c4UGiv=GX%~8Dkntd{!#Jo6`JgxLnWq9cgZzQ zpZvL>ZG+!!{wdMxozw98bi()erO0N7T}6wEUo>aQ!&dcUi0W=t|INP&y*ja7eR2!R ze+qEF;|`zydSTB7!(6w#n`t^z%2PAN`tUqT$n-^=n*X%Y%-)1?@DESS$F=^?#ihof zJH=9h4o}QrzpB#9q((5)YMfeeTH-}B=2(%8@p8n&A)fUeX8i1 zt-MCK-sR{-(YWaT^b?p0=WGY35K zgS(Bv<2(I)eb2vH(VBG)ZZ1IZ6i}AgRP~7#`o^gR`na)r5xj)%@EnNK; z){TztY~xrNSl0NXF+X_qr0vD*mnTl>lC|Y);F#LZ!OS^>aKWC*A#=08*)P;C1jG8A zVgGGFG`SVBlS-I2zNDG9&;{h*Ejjw`9HcjWZ|;*!6t31*H==kVU*y%of-JX+V<<|k zz0}i?XRXv&UoN4;xCwsIr5p7*z_&9t)+M@{d1ihB?DF_@&S!o$>qU)`MRo}4L7-7> zqs4PIc`rJg(!)53nirPQaGpV_@#hG~pfIJ(^G0!*=!~O~DnI3iDV9w!f9%?{b2a|R zZrH-67E8!#8X5#+-xkZ%yCd+2Df}N^=c6tuTvpG0S7N}tZzsL;?x~!{YRVhjt&~JI z`^>owCWr1fe6Fu{d1#lx49$`Yxwj74+>KHRM;T)5I0^FJW643|V7(jSw z>u}FoJGhDBYrh}yPIdB0PQv;U1ENQkv&mcIbuXBb@LX+8ug&j&F<8l1<&{OZ6U3T&-r&Y7u)W5EL6WRxU!IYcq00uX{h9W z^X09}nmaXPydN-qFU>?0u1uKRDQhE$dQ3alL`@)HYE3t3eHYcu7QM$=Q265_d))kF zv1-M$Wlml*?DCY~WmOgF-rT<-e0&h0FV^RRUQ(ZZ!SPK8th`%%4paGVC zmjkW`56Rj`%AKk4D7&q1U@cZ3Z5E; zh{5ZjpXT`*pH@=-+MSbBPkY<@S5olh%^rJP?)BxdSF-v2?w-8qh7XJFV@)1UXzUet z_9wN1i#&(j9U413*LS7}a@_%raZzb=XaSBfHsc$JM*gZcVT_aVa)g7eU z;M1L3HuO8h-62MCU0}X#>f*PBkl$~f=={2uBm?T{`JIzgY-RAf_{WbS@mgoeev>bH zp#0)Y5q)kcqCFQ`L9)O+4U0_=rwG5OhBE)Y+0M3W>ow2k31rb z1S`6i8D*3Fr&3sU8p8sogF;K6{8dfXYs=;1^6Fv|if(BqC6IAeRd1D(ln7xgI;LLD z_+-byKt!QnP^RpiIoyU0JS>C0_05JB%-LFO=AbqxOj8o^fF8*rvs1={f4Uy|1t-SRem*qw=$O;rU`@|TFi?A_(SZBkmXXKfG;O}N&iR?{X5iY+<5`MW z(ev~WNor6fRXNZa!$0OWzDg3xT{wU@*eR=g&lA-FGX366nY z@~`k^sL{NGD&ofYDbJkY??{z6yc&acK16}s4afp%mSX*&ml&Tn)OsJ>0q6q4&Tt{R z{EC+x_?I{K6d7ATK3QNPk8P<6c7}hh9^zUBQ#ESJL+uWhF`8p6+;6I%v> z>$T((?2lkm{RL>Ge55bKIf@f2k*+h`ONVeQ$GfX?P&Sncg&hG)*maA@M|~7hJiJA2 z;|tLiv?qp%U=tIh@CIe%)}MX_>T>(cgkS_!E@4{R;wbaIdn!%&!$5=pkA!wfk+%4t zd>8q2M%rWAAI{KC0X;xnj@kccXTAml(Vl)$BqJEN!Jkn`ujXE8t)jNSIwgeAv`17e z2;JyZqO?wYHG@lqCf6l)O)@d8WKmZh^on!iT7?Z?lgK5aj zSF_zCoC`&~5XmZY1Y(1zt>|}pQVM9u&u;+#Fu-H?*svH&IqlD=ZiK-NC{JMn#Z0(+ z$AjlDXkcfbdyUnNWQ?XzI8Il6QU>Sf%A2@P*R)=ICjz9_5!pPA`3`5B1u`ZPLLIkJzHi(}1vd5`Y5Cml;x&tz30b>Nas~GykX zYqBc!^mGHpCkY-Xr+!$#IHc&eIjuSA&$-W@EuD61V~vQmtVhEP|LXvEvwo?a(zKqc zRNMf>%b>^5n*`sy3dcd*TNJsHr<6^lebhDWkT*|Gh!)NkVhHw+^0xiM(zo(DP_7Bf zE`YW&w>F(BU+7SzT|H?uMe{|-io?vRc_t?z`qJJlrVUA&JA;~`_GOA>wX29eKQi5g zlW7Jb_l-T@Ygl8!F@$8020vMI=~6P(7|1kGDnaT9*1v*$zLW`MJD^GpG}95i8|ULB zeHqVD<{&^a2CebeL*p-o*DbOVEpF_=*h5ll1&Ld#+mOCO8rj!7&FI9<~ zqDlh<0dbzCm7yFAr~#J^q`~~A4E*OPFR9m(&jJsSFGu4pQohmo0%}EgyCVJSX<4LD z<=@BN#-hM2n}nzOXEIs^-}fTg(;tZ)c2gq0BJOp98A&j|rbI+CtupOn_dGaDwR>0TFvl3-(rr>fGuv)3iR-&zM#oAXQc{Z&~@E6rIyZ zXEwC}Xl8FMhN*A;ev+a{m$tWm0c)aV|rBFVy zHF+O&M~LEIb_#ZOk>BWxO^tSq^)J-a(o)AEJ8~2JpuTx0;p47iyO?TAge0YsrX(C8 zfS*jrj0zW`$G8xfZg@f!_r@|0Fyp^RnU=sKhsYuAi!~T0H(Tl8hd4HU2V|Rl^gqL( zt=;C@^2zS8b;y}P@q}i03m2fTV~Hc0yhIhz)2~othox%jhH-Riza%V(f^jlBm35F5 zz;rT)$1lYIbv0tiqzLN|b%vbWT6#$1200O5%2Xa8b0wPjTA+p%Fbru&dy!OZbziEy z0S@VAB!5IZMvjBX=kU+;pCX`45MHbc5bD8o5c&(VWwoT6OKM8bm`7nZgDH!W!jI74sFbXZ6p6kYydW zV)psf7)v7lz&eB#vJR|jp+2FOeS9aJvt#rI1ZE3E&{@U!i;~e9|D!a zhVw>r;l9W0o2CyTKs4~2v+GDPIx5M5HMfs)c=(!`Z?5;YtI+`qX2qi0<4FacxMZTT z1AJ$wNx0_NgAQSZvdb8CenUO9I6+fg5(R8}0+xJlNu!f&+K@MxY(53Df+Lu8JD%tf zbh)Tx2fqBuEhYeSu@7zRSlezO=tlaqdNcZFzYY6^Oz73WMWY+3&j3Rx95r|y9ShNz z9+g{TMUk=#10#IF$ZEGUz80pdZ6N5IB!uohwL}r7fCngd{NN zGt{e!{fcWSeIefC1n@!zv%tgT>LYv;9hDTpc=0c(kM3?@bIL!y(?aDC7=@{mnqIvZwAUi zWjL|66pImLNj^@D5l(B43ueIU>(3%e$puy@fx%qbZb}h#IoH;v>!EMdbrb2qfZw#p zM1<(~7vvbcw&Q&+=j<80i4i{Y$WCZDv=s!7E9s@eiRW%Hrw&nR?hbr=xaqAR6149=fWQT`zBgn zY!?y7{$~(Z0&u2EfHOH$l<;0l;y@&aUHpfWWL<0}Z3(IMLG-j|NG(P47ZsrPUz(8G zNgDnSNWPEr|J&ty&GiM_OHIHNXgwR37vQUk(VK0}59Sdcl8()+LgWXvXu2s)z_TPZ zy&UHcM~STx55Jn+goC)qlKaoWvaNsRM}ANwLjME&yD}6U#OHK?9G=FRPDz;&bc;Kf z?VM<3(5JM4S|@x!R8uYhDWp{h0Tp+2Djy~%8Qtd6vb|4?up+60*@En7{RuG}PLnUrSh6px%uR1h0Te|DC zvez1aesEA=kGxrbVB-W-PLzOlo2$#D$#rWy9rahleDjk*D$=kI`J!Gl37DQ15q1iI ztV|>90CdhJZ7+9jI#mu5hYpM!hr|CCDU!qE2-Vxj2gy(h6GXBQe~baPk}EV6@~vnX z)mFhDTT=}i@1f&bXt*c^e&*jG5x~-|C6Xq_}1I08pv`f{2wib)3Y#oPv@M+BJ&Bz${9qO|-iVJlf+78;JO&WI z&M%1r#f%1aUjRi)elNH6%K$I|qUHdmT z8w~I;hCL#ffm~-LUuK{>zz-Q94DT^8z0G{}h6MD!0GgaY1mElB5IFiMOt(sqZ8*1G z1jqZnBE5|PuL)8e&q1C2=BbB=(^at<|s`Q#nt`#6onQy~_J5V9Y|-HEgSWdfK+I?Oq6 zro=bDugcD$A&@WqrYE}eJiubusF@Mo9sTWV>lyRyV>;wWGv~H1lP+EB43h3JK(_Dz zQQZz>mdKxr5gOIAsV6GwG$Jn17if;apE@dS0f0XNMDGJhk`1jGZ8!?^ z`Y1r7@BjmB+05o^1=E#S{S5HW=1ISkX_8#)%}(PJUmB9`4Nnlo0C*RO`xBs(NND6K zl4k$Sw9v@OmRUBklPKCbD<(m{xpu6Cm@c9Ot-)OL|2y)~!!}$}$zm(`4MeaQGd%z7QSYsS zUw=Ev>i&!rJt6YhE0l1mscQQdV;zzPZ2ezBd?(B$zElJT8XT@tzq!Us{(!cQ9FL&D z#=)O~iBbB95vIWi)EtskhZDe<1+fS|yqx|Fl68#yherWMpSHKPPfcT%rP>y&r41ly&JO#E!urrzI2m>^7P-y?bmGwWkdP3O&_!!X00pF$K z-fM#JEJFCqG%y=bQ~?!s_<_g<=K7h`RAvVUK>yU23Rj+MTy@Yl)HPhh@&X%qa%90^ z-=W7LLaY*NnK2F2Kp<-U6fr z(@PxQMh;Vlou|>I>m%Dixo-eu9k~7jR}0&1Rr)$8>T9y@slg9E=v$Kb8&a(wBenSg zEQxLWH!6^?#+{%CaGC)iFDl*1ceceYZYM_r$`GKkYyxGo#fg{C`3He~PJOulUDum% zl~6SGZVIp1xAodje8wPr4@x_e`~)AJzX#rCM|3GuLPAs5n;zPkHc!h`!8(SAxRlIT zj?J9qic~XzYKVBQRB3k_c-=k#V4=ve_uLxZq0dCR6dYlusBcw{T2uh1eZHikq7Ik40{K6fq9I<=pAln6vOrX*g$?Rr28>4H)T@dM$p8Uri_zSufAjt- zWY*b4Z5r9AtW|=<1FSMw9y?j6*t*iarIlZYT$yMLY@P$=hT1heqbeld{_kWMwivvhVH6 zsgK;MEu6)C4hj&{?e*&$EBKmbY&_8Z0YDc6HeS#l0%HBY$D60%OMmV(m5V9wj6@;o@({P9TkQVfl%wuR3VHEgWfWFOa^RF9o$|IH(+m&;U~DjF z0&l{@Gi6W`obA{%#E_!<3jHf33j+9Uo7f}}nLwM>$NRgA zUx+}>NNpiG3ol?ph~d5T&x?uC!Ow>n7W^p)%L0Pd0wl%&taz?BHd2upP~ocjK{zoq zPQ*E)%P$cuY92?+Lb|ZVLfa%4)bIe`dYPQWj}ejup(X0kjCobc{HV3MMj0 zGkP#9weZ5OD1ey&Hh2tFM!6R$4U}Vn!G{1L6FUSdH!*rR_dzt&0T62QNOzq$XHq6wE!1r2M0a$c7D! zq=iGFY(PP?0lrQqRy{lJ8h$_2SY22&@JMsU;IQPBX7LATepzX{gna3LUcV$t$SbaLOIv74)q@ z`}KeY_89 z*#?jm9?}!N+9Sg$_v@UZ2d}rzAUJ?-vu-leQn%SGxbQ~{9exmCKF{djv1gCgjb2*M zQ~?6+a3O*f;7)udr%@!e<*K7>H<6oco8+>`%=>jU{7whKxRZj_n}E0>G)jsCx=9~_ z#{*KndsQ`3oW46(N9kqG;wlVVb|-zw6Wa5~+kX~8U!nG(t>eXpu=#w)*ES*oGjQm7 ze$>yM@Rp!{vAt4dgy(57U4J>bJ^ly2Gkvu4m!F6JlqvYYtnA5x?zy_Fi5*$2KwHn5 z{6MCW&zthLZpbK~6clo?#QuvNcm$n4A6vppvX) zKH;M1TO2PBNdHN@ynhhXAQ_@x)YUDM`G^uY(Lcgo^Now3vi-2sYeu=aB!o{w**OC8 z=?p=6cqBoj-a~IiQ!%=zyj+&u&7OXM-DPM?eJkeLJk6EW1l+49+Usr|+PyLSUTp!J zT@9)pm=cc=-~6^r^5r_^fZazPNgpPv1{;nd#1Emn#}^gfet-Hy&5z$tF1x4Vmew>$ zf4wQbc1uU=+vrm*Wslo=8v0kdbdJA+{C)_u1?JH2W+x7Onwnf!x&a;ei)p(;KHA}s>xvHshlxJiP6*=&Pcl|ADbO;Audj^ePs+AeqHFKUm+gwC;4M*wEN>uU=kD}(o z<>!sZ-RQjG%Fzjs>dddFIj2gltjYy_=zz7mQr@DXrkrAM^NsY9lv}g5M-mi{hP~*E z=t6G(n#E4mj8xcbf_QCfLHhUmZawGUp;Pu8tS$wHTw>f^ihuOXp+>d(=vaDYLY(t4 z&-dGK<;?T@IZty3j2)Qmj-;IqGKm+y&XnG?GQ2(sFLZMagtVs9MiHE8_EI zjqIdE-WjK{w|`b2jz|{0OcWsth355q{2Ov9cZlP(B3&vr@FnWWIGA~>bL_?^n<69Y z_AY#(+~=Pk;u>rkTiZ`IY*wXWtR(WD`;bH2{r;|Ah!Pnby&I}fAT~KCV`K8Y^B^_1RrAyBDyfuZc-i`FmU3|iA+>8*q-!vAyYv6yRRgAbRFZfk=6`Sk~b4`r@x z4qJY0GQYOS--<6oxn#>7jku`M#~!$yYwqry4XCIzKmMKG*gBiw(bMBjoc(Ua+x^Jj@NOJsgzE$Uqc`Q|hg~Gf&}M~% zs;0d_LVU%5BdObNwBWp=TvgrQC#uTVMDyg>Cj^__7ijDqM&N{^I?HgfnxhDM%z4wd zAAcaf3OVDZxV<;+2ql;BWyWifTtckNoC-&irFwc&XCG^4D8(veL#`AgCXUwp4peWk zP^?LwP1gjUKd`(kVbS{#r+qwIbn@07~7$d(}aT*HH8yXvOt zrY)!45semZ%j}_R%IVq}SIwkcWQLs!EUFsy$b|5NqO{em!xA-$v|; zP*+b+tjLsN=6Yy?LSm^^Q~~O$9j}O?djMesqI!0)02OH`RPtGS0`>ga_5AF+>t+=U zq?WHiGcPhXF@Kr`O%pvzb!gc+ujh_8PKHY!3Mp&2tN{H?&0fd$Ky@;>>`HHl9M$_M(hxt6?W7~dV)4fLKi zQTcUFNqN$JfH1t{nK2dHYL_FtQ#J#*`dP{iCzp#ykC}$=3M|IY2WkSwf&=F z_9)h1>7aN*{z1oWbSmy%miZ_V?iwVNm#6PhHa6>d5d1<&N9_3J;O7mF#lR$!Vwb5! z4{g`sY@-h+?rSsKQhA&nY8&siR+?URR`Vq$Nt^FaArjd}UCJv;a~ZXFizL^Kt<#>mbS2${CM{}OHjUs$Su*DA zRgC0s8{8Tnn|Z+z{C4F!E>L`MYa36?3&pIyuqghT=9=>-D=PP?i{07UPN~m#zU?@9 z)(<0s zm?lw=Dy4uUE$SZ|v$?36FyKh{*hG`O`*>C*Vfs1bB>YU*d^hpI2kY6LjuZZ`8kJ;u zTT$$rw~~rNeJCaEVyUy!jStyHBTAL_I}Us=`El6MEeWH(ALh?W`pT?gXMs;dbMFv`NJ;O7>g=b_U--Z$U#cp7-tCWd9&+j? zrID{&(u38$LP*(3vlBUzBi`>6x9#-w{f`@*LS)b!^9w)GMEvZAThPmF*$LykHNnko zZDp;|khSdL1!sT_f1S;6j0|0CEe)9qmAco0)wN9A`V<*4&N{LCV>HbYwf#v;JU9xS zh^xUrEU5CSdlhbAbF;a{VHSL^f?c5^~m;b f=OI%OW!Ftx&F+