diff --git a/NAMESPACE b/NAMESPACE index 2f2fffb..faa1358 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -3,16 +3,21 @@ S3method(print,group_time_att) S3method(print,pte_results) S3method(print,summary.pte_results) +S3method(summary,aggte_obj) S3method(summary,group_time_att) S3method(summary,pte_emp_boot) S3method(summary,pte_results) +export(aggte_obj) export(attgt_if) export(attgt_noif) export(attgt_pte_aggregations) export(compute.pte) export(compute.pte2) +export(crit_val_checks) export(did_attgt) +export(dose_obj) export(ggpte) +export(ggpte_cont) export(group_time_att) export(gt_data_frame) export(keep_all_pretreatment_subset) diff --git a/R/ggpte.R b/R/ggpte.R index 74da757..efca0ac 100644 --- a/R/ggpte.R +++ b/R/ggpte.R @@ -8,12 +8,57 @@ ggpte <- function(pte_results) { plot_df <- summary(pte_results)$event_study colnames(plot_df) <- c("e", "att", "se", "cil", "ciu") - plot_df$post <- as.factor(1*(plot_df$e >= 0)) - ggplot(plot_df, aes(x=e, y=att)) + - geom_line(aes(color=post)) + - geom_point(aes(color=post)) + - geom_line(aes(y=ciu), linetype="dashed", alpha=0.5) + - geom_line(aes(y=cil), linetype="dashed", alpha=0.5) + + plot_df$post <- as.factor(1 * (plot_df$e >= 0)) + ggplot(plot_df, aes(x = e, y = att)) + + geom_line(aes(color = post)) + + geom_point(aes(color = post)) + + geom_line(aes(y = ciu), linetype = "dashed", alpha = 0.5) + + geom_line(aes(y = cil), linetype = "dashed", alpha = 0.5) + theme_bw() + - theme(legend.position="bottom") + theme(legend.position = "bottom") +} + + +#' @title ggpte_cont +#' +#' @description a function for plotting results in applications with a continuous treatment +#' +#' @param dose_obj a `dose_obj` that holds results with a continuous treatment +#' @param type whether to plot ATT(d) or ACRT(d), defaults to `att` for +#' plotting ATT(d). For ACRT(d), use "acrt" +#' +#' @export +ggpte_cont <- function(dose_obj, type = "att") { + dose <- dose_obj$dose + if (type == "acrt") { + acrt.d <- dose_obj$acrt.d + acrt.d_se <- dose_obj$acrt.d_se + acrt.d_crit.val <- dose_obj$acrt.d_crit.val + plot_df <- cbind.data.frame(dose, acrt.d, acrt.d_se, acrt.d_crit.val) + ggplot(plot_df, aes(x = dose, y = acrt.d)) + + geom_line(size = 2) + + geom_ribbon( + aes( + ymin = acrt.d - acrt.d_crit.val * acrt.d_se, + ymax = acrt.d + acrt.d_crit.val * acrt.d_se + ), + fill = "lightgray", alpha = 0.5 + ) + + theme_bw() + } else { # att(d) plot + att.d <- dose_obj$att.d + att.d_se <- dose_obj$att.d_se + att.d_crit.val <- dose_obj$att.d_crit.val + plot_df <- cbind.data.frame(dose, att.d, att.d_se, att.d_crit.val) + ggplot(plot_df, aes(x = dose, y = att.d)) + + geom_line(size = 2) + + geom_ribbon( + aes( + ymin = att.d - att.d_crit.val * att.d_se, + ymax = att.d + att.d_crit.val * att.d_se + ), + fill = "lightgray", alpha = 0.5 + ) + + theme_bw() + } } diff --git a/R/process_dose_gt.R b/R/process_dose_gt.R index 55f6b23..b2fba98 100644 --- a/R/process_dose_gt.R +++ b/R/process_dose_gt.R @@ -1,9 +1,165 @@ process_dose_gt <- function(gt_results, ptep, ...) { - browser() - # make the call to att, to get same format of results att_gt <- process_att_gt(gt_results, ptep) o_weights <- overall_weights(att_gt, ...) - 1 + 1 + # main dose-specific results are in extra_gt_returns + all_extra_gt_returns <- att_gt$extra_gt_returns + groups <- unlist(BMisc::getListElement(all_extra_gt_returns, "group")) + time.periods <- unlist(BMisc::getListElement(all_extra_gt_returns, "time.period")) + + # check that order of groups and time periods matches + if (!all(cbind(groups, time.periods) == o_weights[, c("group", "time.period")])) { + stop("in processing dose results, mismatch between order of groups and time periods") + } + + inner_extra_gt_returns <- BMisc::getListElement(all_extra_gt_returns, "extra_gt_returns") + att.d_gt <- BMisc::getListElement(inner_extra_gt_returns, "att.d") + acrt.d_gt <- BMisc::getListElement(inner_extra_gt_returns, "acrt.d") + att.overall_gt <- unlist(BMisc::getListElement(inner_extra_gt_returns, "att.overall")) + acrt.overall_gt <- unlist(BMisc::getListElement(inner_extra_gt_returns, "acrt.overall")) + bet_gt <- BMisc::getListElement(inner_extra_gt_returns, "bet") + bread_gt <- BMisc::getListElement(inner_extra_gt_returns, "bread") + Xe_gt <- BMisc::getListElement(inner_extra_gt_returns, "Xe") + + # point estimates of ATT(d) and ACRT(d) + att.d <- weighted_combine_list(att.d_gt, o_weights$overall_weight) + acrt.d <- weighted_combine_list(acrt.d_gt, o_weights$overall_weight) + + # values of the dose + dvals <- ptep$dvals + degree <- ptep$degree + knots <- ptep$knots + bs_grid <- splines2::bSpline(dvals, degree = degree, knots = knots) + bs_grid <- cbind(1, bs_grid) # add intercept + bs_deriv <- splines2::dbs(dvals, degree = degree, knots = knots) + bs_deriv <- cbind(0, bs_deriv) # intercept doesn't matter here, just placeholder to get dimensions right + + # since we are picking dvals over a grid, the only randomness comes from + # estimating the \beta's + n1_vec <- sapply(Xe_gt, nrow) + acrt_gt_inffunc_mat <- gt_results$inffunc + n <- nrow(acrt_gt_inffunc_mat) + keep_mat <- acrt_gt_inffunc_mat != 0 + if (!all(colSums(keep_mat) == n1_vec)) { + stop("something off with overall influence function") + } + + att.d_gt_inffunc <- lapply( + 1:length(Xe_gt), + function(i) { + out_inffunc <- matrix(data = 0, nrow = n, ncol = length(dvals)) + this_inffunc <- (Xe_gt[[i]] %*% bread_gt[[i]] %*% t(bs_grid)) + out_inffunc[keep_mat[, i], ] <- (n / n1_vec[i]) * this_inffunc + out_inffunc + } + ) + + att.d_inffunc <- weighted_combine_list(att.d_gt_inffunc, o_weights$overall_weight) + biters <- ptep$biters + alp <- ptep$alp + cband <- ptep$cband + boot_res <- mboot2(att.d_inffunc, biters = biters, alp = alp) + att.d_se <- boot_res$boot_se + if (cband) { + att.d_crit.val <- boot_res$crit_val + att.d_crit.val <- crit_val_checks(att.d_crit.val, alp) + } else { + att.d_crit.val <- qnorm(1 - alp / 2) + } + + # influence function for acrt + + # acrt influence function - same as for att.d except use derivative of basis functions + acrt.d_gt_inffunc <- lapply( + 1:length(Xe_gt), + function(i) { + out_inffunc <- matrix(data = 0, nrow = n, ncol = length(dvals)) + this_inffunc <- (Xe_gt[[i]] %*% bread_gt[[i]] %*% t(bs_deriv)) + out_inffunc[keep_mat[, i], ] <- (n / n1_vec[i]) * this_inffunc + out_inffunc + } + ) + acrt.d_inffunc <- weighted_combine_list(acrt.d_gt_inffunc, o_weights$overall_weight) + acrt_boot_res <- mboot2(acrt.d_inffunc, biters = biters, alp = alp) + acrt.d_se <- acrt_boot_res$boot_se + if (cband) { + acrt.d_crit.val <- acrt_boot_res$crit_val + acrt.d_crit.val <- crit_val_checks(acrt.d_crit.val, alp) + } else { + acrt.d_crit.val <- qnorm(1 - alp / 2) + } + + # placeholder for tracking `call` + call <- NULL + + dose_obj( + dose = dvals, + att.d = att.d, + att.d_se = att.d_se, + att.d_crit.val = att.d_crit.val, + att.d_inffunc = att.d_inffunc, + acrt.d = acrt.d, + acrt.d_se = acrt.d_se, + acrt.d_crit.val = acrt.d_crit.val, + acrt.d_inffunc = acrt.d_inffunc, + pte_params = ptep, + call = call + ) +} + +#' @title dose_obj +#' +#' @description Holds results from computing dose-specific treatment effects +#' with a continuous treatment +#' +#' @param dose vector containing the values of the dose used in estimation +#' @param att.d estimates of ATT(d) for each value of `dose` +#' @param att.d_se standard error of ATT(d) for each value of `dose` +#' @param att.d_crt.val critical value to produce pointwise or uniform confidence +#' interval for ATT(d) +#' @param att.d_inffunc matrix containing the influence function from estimating +#' ATT(d) +#' @param acrt.d estimates of ACRT(d) for each value of `dose` +#' @param acrt.d_se standard error of ACRT(d) for each value of `dose` +#' @param acrt.d_crt.val critical value to produce pointwise or uniform confidence +#' interval for ACRT(d) +#' @param acrt.d_inffunc matrix containing the influence function from estimating +#' ACRT(d) +#' @param pte_params a pte_params object containing other parameters passed to the function +#' @param call the original call to the function for computing causal effect parameters +#' with a continuous treatment +#' +#' @return dose_obj +#' +#' @export +dose_obj <- function( + dose, + att.d = NULL, + att.d_se = NULL, + att.d_crit.val = NULL, + att.d_inffunc = NULL, + acrt.d = NULL, + acrt.d_se = NULL, + acrt.d_crit.val = NULL, + acrt.d_inffunc = NULL, + pte_params = NULL, + call = NULL) { + out <- list( + dose = dose, + att.d = att.d, + att.d_se = att.d_se, + att.d_crit.val = att.d_crit.val, + att.d_inffunc = att.d_inffunc, + acrt.d = acrt.d, + acrt.d_se = acrt.d_se, + acrt.d_crit.val = acrt.d_crit.val, + acrt.d_inffunc = acrt.d_inffunc, + pte_params = pte_params, + call = call + ) + + class(out) <- "dose_obj" + + out } diff --git a/R/pte.R b/R/pte.R index 28b6fdf..7c097a7 100644 --- a/R/pte.R +++ b/R/pte.R @@ -421,7 +421,15 @@ pte <- function(yname, max_e <- ifelse(is.null(dots$max_e), Inf, dots$max_e) balance_e <- dots$balance_e - event_study <- pte_aggte(att_gt, type = "dynamic", bstrap = TRUE, cband = cband, alp = ptep$alp, min_e = min_e, max_e = max_e, balance_e = balance_e) + event_study <- pte_aggte(att_gt, + type = "dynamic", + bstrap = TRUE, + cband = cband, + alp = ptep$alp, + min_e = min_e, + max_e = max_e, + balance_e = balance_e + ) # output out <- pte_results( diff --git a/R/pte_aggte.R b/R/pte_aggte.R index 4112aba..95e4425 100644 --- a/R/pte_aggte.R +++ b/R/pte_aggte.R @@ -110,21 +110,7 @@ pte_aggte <- function(attgt, } selective.crit.val <- mboot2(selective.inf.func.g, biters = biters, alp = alp)$crit_val - if (is.na(selective.crit.val) | is.infinite(selective.crit.val)) { - warning("Simultaneous critival value is NA. This probably happened because we cannot compute t-statistic (std errors are NA). We then report pointwise conf. intervals.") - selective.crit.val <- stats::qnorm(1 - alp / 2) - dp$cband <- FALSE - } - - if (selective.crit.val < stats::qnorm(1 - alp / 2)) { - warning("Simultaneous conf. band is somehow smaller than pointwise one using normal approximation. Since this is unusual, we are reporting pointwise confidence intervals") - selective.crit.val <- stats::qnorm(1 - alp / 2) - dp$cband <- FALSE - } - - if (selective.crit.val >= 7) { - warning("Simultaneous critical value is arguably `too large' to be realible. This usually happens when number of observations per group is small and/or there is no much variation in outcomes.") - } + selective.crit.val <- crit_val_checks(selective.crit.val, alp) } # get overall att under selective treatment timing @@ -157,7 +143,7 @@ pte_aggte <- function(attgt, if ((selective.se <= sqrt(.Machine$double.eps) * 10)) selective.se <- NA } - return(AGGTEobj( + return(aggte_obj( overall.att = selective.att, overall.se = selective.se, type = type, @@ -237,21 +223,7 @@ pte_aggte <- function(attgt, } dynamic.crit.val <- mboot2(dynamic.inf.func.e, biters = biters, alp = alp)$crit_val - if (is.na(dynamic.crit.val) | is.infinite(dynamic.crit.val)) { - warning("Simultaneous critival value is NA. This probably happened because we cannot compute t-statistic (std errors are NA). We then report pointwise conf. intervals.") - dynamic.crit.val <- stats::qnorm(1 - alp / 2) - dp$cband <- FALSE - } - - if (dynamic.crit.val < stats::qnorm(1 - alp / 2)) { - warning("Simultaneous conf. band is somehow smaller than pointwise one using normal approximation. Since this is unusual, we are reporting pointwise confidence intervals") - dynamic.crit.val <- stats::qnorm(1 - alp / 2) - dp$cband <- FALSE - } - - if (dynamic.crit.val >= 7) { - warning("Simultaneous critical value is arguably `too large' to be realible. This usually happens when number of observations per group is small and/or there is no much variation in outcomes.") - } + dynamic.crit.val <- crit_val_checks(dynamic.crit.val, alp) } # get overall average treatment effect @@ -271,7 +243,7 @@ pte_aggte <- function(attgt, if (dynamic.se <= sqrt(.Machine$double.eps) * 10) dynamic.se <- NA } - return(AGGTEobj( + return(aggte_obj( overall.att = dynamic.att, overall.se = dynamic.se, type = type, @@ -401,7 +373,6 @@ overall_weights <- function(attgt, min_e = -Inf, max_e = Inf, ...) { - browser() group <- attgt$group time.period <- attgt$t att <- attgt$att @@ -468,5 +439,223 @@ overall_weights <- function(attgt, # quick sanity check if (sum(out_weight) != 1) stop("something's going wrong calculating overall weights") - data.frame(group = group, time.period = time.period, overall_weight = out_weight) + # following convention of package, return groups and time periods in their original scale + data.frame( + group = t2orig(group, originaltlist), + time.period = t2orig(time.period, originaltlist), + overall_weight = out_weight + ) +} + +#' @title crit_val_checks +#' +#' @description A function to perform sanity checks and possibly adjust a +#' a critical value to form a uniform confidence band +#' +#' @param crit_val the critical value +#' @param alp the significance level +#' +#' @return a (possibly adjusted) critical value +#' +#' @export +crit_val_checks <- function(crit_val, alp = 0.05) { + if (is.na(crit_val) | is.infinite(crit_val)) { + warning("Simultaneous critival value is NA. This probably happened because we cannot compute t-statistic (std errors are NA). We then report pointwise conf. intervals.") + crit_val <- stats::qnorm(1 - alp / 2) + dp$cband <- FALSE + } + + if (crit_val < stats::qnorm(1 - alp / 2)) { + warning("Simultaneous conf. band is somehow smaller than pointwise one using normal approximation. Since this is unusual, we are reporting pointwise confidence intervals") + crit_val <- stats::qnorm(1 - alp / 2) + dp$cband <- FALSE + } + + if (crit_val >= 7) { + warning("Simultaneous critical value is arguably `too large' to be realible. This usually happens when number of observations per group is small and/or there is no much variation in outcomes.") + } + + crit_val +} + + + +#' @title aggte_obj +#' +#' @description Objects of this class hold results on aggregated +#' group-time average treatment effects. This is derived from the AGGTEobj +#' class in the `did` package. +#' +#' @title Aggregate Treatment Effect Parameters Object +#' +#' @description An object for holding aggregated treatment effect parameters. +#' +#' @inheritParams aggte +#' @inheritParams compute.aggte +#' @param overall.att The estimated overall ATT +#' @param overall.se Standard error for overall ATT +#' @param egt Holds the length of exposure (for dynamic effects), the +#' group (for selective treatment timing), or the time period (for calendar +#' time effects) +#' @param att.egt The ATT specific to egt +#' @param se.egt The standard error specific to egt +#' @param crit.val.egt A critical value for computing uniform confidence +#' bands for dynamic effects, selective treatment timing, or time period +#' effects. +#' @param inf.function The influence function of the chosen aggregated parameters +#' @param DIDparams A DIDparams object +#' +#' @return an aggte_obj +#' @export +aggte_obj <- function(overall.att = NULL, + overall.se = NULL, + type = "simple", + egt = NULL, + att.egt = NULL, + se.egt = NULL, + crit.val.egt = NULL, + inf.function = NULL, + min_e = NULL, + max_e = NULL, + balance_e = NULL, + call = NULL, + DIDparams = NULL) { + out <- list( + overall.att = overall.att, + overall.se = overall.se, + type = type, + egt = egt, + att.egt = att.egt, + se.egt = se.egt, + crit.val.egt = crit.val.egt, + inf.function = inf.function, + min_e = min_e, + max_e = max_e, + balance_e = balance_e, + call = call, + DIDparams = DIDparams + ) + + class(out) <- "aggte_obj" + out +} + +#' @title Summary Aggregate Treatment Effect Parameter Objects +#' +#' @description A function to summarize aggregated treatment effect parameters. +#' +#' @param object an \code{aggte_obj} object +#' @param ... other arguments +#' +#' @export +summary.aggte_obj <- function(object, ...) { + # call + cat("\n") + cat("Call:\n") + print(object$call) + cat("\n") + + # citation + citation() + cat("\n") + + # overall estimates + alp <- object$DIDparams$alp + pointwise_cval <- qnorm(1 - alp / 2) + overall_cband_upper <- object$overall.att + pointwise_cval * object$overall.se + overall_cband_lower <- object$overall.att - pointwise_cval * object$overall.se + out1 <- cbind.data.frame(object$overall.att, object$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) + cat("\n") + # cat("Overall ATT: \n") + if (object$type == "dynamic") cat("Overall summary of ATT\'s based on event-study/dynamic aggregation: \n") + if (object$type == "group") cat("Overall summary of ATT\'s based on group/cohort aggregation: \n") + if (object$type == "calendar") cat("Overall summary of ATT\'s based on calendar time aggregation: \n") + colnames(out1) <- c("ATT", " Std. Error", paste0(" [ ", 100 * (1 - object$DIDparams$alp), "% "), "Conf. Int.]", "") + print(out1, row.names = FALSE) + cat("\n\n") + + # handle cases depending on type + if (object$type %in% c("group", "dynamic", "calendar")) { + # header + if (object$type == "dynamic") { + c1name <- "Event time" + cat("Dynamic Effects:") + } + if (object$type == "group") { + c1name <- "Group" + cat("Group Effects:") + } + if (object$type == "calendar") { + c1name <- "Time" + cat("Time Effects:") + } + + cat("\n") + cband_text1a <- paste0(100 * (1 - object$DIDparams$alp), "% ") + cband_text1b <- ifelse(object$DIDparams$bstrap, + ifelse(object$DIDparams$cband, "Simult. ", "Pointwise "), + "Pointwise " + ) + cband_text1 <- paste0("[", cband_text1a, cband_text1b) + + cband_lower <- object$att.egt - object$crit.val.egt * object$se.egt + cband_upper <- object$att.egt + object$crit.val.egt * object$se.egt + + sig <- (cband_upper < 0) | (cband_lower > 0) + sig[is.na(sig)] <- FALSE + sig_text <- ifelse(sig, "*", "") + + out2 <- cbind.data.frame(object$egt, object$att.egt, object$se.egt, cband_lower, cband_upper) + out2 <- round(out2, 4) + + out2 <- cbind.data.frame(out2, sig_text) + + colnames(out2) <- c(c1name, "Estimate", "Std. Error", cband_text1, "Conf. Band]", "") + print(out2, row.names = FALSE, justify = "centre") + } + cat("---\n") + cat("Signif. codes: `*' confidence band does not cover 0") + cat("\n\n") + + # set control group text + control_group <- object$DIDparams$control_group + control_group_text <- NULL + if (control_group == "nevertreated") { + control_group_text <- "Never Treated" + } else if (control_group == "notyettreated") { + control_group_text <- "Not Yet Treated" + } + + if (!is.null(control_group)) { + cat("Control Group: ") + cat(control_group_text) + cat(", ") + } + + # anticipation periods + cat("Anticipation Periods: ") + cat(object$DIDparams$anticipation) + cat("\n") + + # estimation method text + est_method <- object$DIDparams$est_method + if (is(est_method, "character")) { + est_method_text <- est_method + if (est_method == "dr") { + est_method_text <- "Doubly Robust" + } else if (est_method == "ipw") { + est_method_text <- "Inverse Probability Weighting" + } else if (est_method == "reg") { + est_method_text <- "Outcome Regression" + } + + cat("Estimation Method: ") + cat(est_method_text) + cat("\n") + } } diff --git a/R/pte_params.R b/R/pte_params.R index 6362f66..12d2a22 100644 --- a/R/pte_params.R +++ b/R/pte_params.R @@ -92,7 +92,7 @@ setup_pte_basic <- function(yname, #' and \code{attgt_fun} in the next steps. #' #' @inheritParams pte_params -#' @param required_pre_periods The number of required pre-treatment periods to implement +#' @param required_pre_periods The number of required pre-treatment periods to implement #' the estimation strategy. Default is 1. #' @param ... additional arguments #' diff --git a/README.md b/README.md index 3dcf314..0f3a3e6 100644 --- a/README.md +++ b/README.md @@ -82,23 +82,26 @@ did_res <- pte(yname="lemp", subset_fun=two_by_two_subset, attgt_fun=did_attgt, xformla=~lpop) +#> Warning in pte(yname = "lemp", gname = "first.treat", tname = "year", idname = "countyreal", : 'pte' is deprecated. +#> Use 'pte2' instead. +#> See help("Deprecated") summary(did_res) #> #> Overall ATT: #> ATT Std. Error [ 95% Conf. Int.] -#> -0.0323 0.0126 -0.0569 -0.0076 * +#> -0.0323 0.0126 -0.057 -0.0075 * #> #> #> Dynamic Effects: #> Event Time Estimate Std. Error [95% Conf. Band] -#> -3 0.0269 0.0168 -0.0159 0.0697 -#> -2 -0.0050 0.0125 -0.0368 0.0269 -#> -1 -0.0229 0.0126 -0.0550 0.0093 -#> 0 -0.0201 0.0131 -0.0534 0.0132 -#> 1 -0.0547 0.0179 -0.1003 -0.0091 * -#> 2 -0.1382 0.0409 -0.2422 -0.0341 * -#> 3 -0.1069 0.0341 -0.1936 -0.0202 * +#> -3 0.0269 0.0154 -0.0131 0.0670 +#> -2 -0.0050 0.0129 -0.0385 0.0286 +#> -1 -0.0229 0.0144 -0.0603 0.0145 +#> 0 -0.0201 0.0117 -0.0506 0.0103 +#> 1 -0.0547 0.0165 -0.0976 -0.0118 * +#> 2 -0.1382 0.0396 -0.2410 -0.0354 * +#> 3 -0.1069 0.0341 -0.1954 -0.0184 * #> --- #> Signif. codes: `*' confidence band does not cover 0 ggpte(did_res) @@ -145,6 +148,9 @@ covid_res <- pte(yname="positive", xformla=xformla, max_e=21, min_e=-10) +#> Warning in pte(yname = "positive", gname = "group", tname = "time.period", : 'pte' is deprecated. +#> Use 'pte2' instead. +#> See help("Deprecated") #> Warning in compute.aggte(MP = MP, type = type, balance_e = balance_e, min_e = #> min_e, : Simultaneous conf. band is somehow smaller than pointwise one using #> normal approximation. Since this is unusual, we are reporting pointwise @@ -154,43 +160,43 @@ summary(covid_res) #> #> Overall ATT: #> ATT Std. Error [ 95% Conf. Int.] -#> 14.8882 58.1647 -99.1125 128.8888 +#> 14.8882 65.6876 -113.8573 143.6336 #> #> #> Dynamic Effects: #> Event Time Estimate Std. Error [95% Conf. Band] -#> -10 -3.7266 3.3841 -13.4426 5.9893 -#> -9 2.6607 1.5017 -1.6509 6.9722 -#> -8 0.8290 2.1465 -5.3339 6.9919 -#> -7 5.2843 1.9978 -0.4516 11.0203 -#> -6 2.8555 1.8211 -2.3730 8.0841 -#> -5 1.3589 4.0879 -10.3778 13.0956 -#> -4 0.3294 3.9431 -10.9916 11.6504 -#> -3 -4.2227 6.9336 -24.1295 15.6841 -#> -2 -3.8447 2.4340 -10.8329 3.1434 -#> -1 -0.2234 3.3015 -9.7021 9.2553 -#> 0 -10.8156 8.8849 -36.3249 14.6936 -#> 1 -13.7998 12.2977 -49.1073 21.5077 -#> 2 -7.8432 11.9865 -42.2575 26.5711 -#> 3 -4.5541 11.7315 -38.2362 29.1279 -#> 4 -3.5368 12.7993 -40.2847 33.2110 -#> 5 8.5221 11.1569 -23.5102 40.5544 -#> 6 1.1140 17.8634 -50.1732 52.4012 -#> 7 6.6384 24.7105 -64.3072 77.5841 -#> 8 7.1288 23.8552 -61.3613 75.6189 -#> 9 10.8758 30.7177 -77.3170 99.0687 -#> 10 17.5057 34.2477 -80.8221 115.8335 -#> 11 40.8318 46.6352 -93.0614 174.7250 -#> 12 48.6134 57.4288 -116.2692 213.4960 -#> 13 52.4228 64.1772 -131.8350 236.6806 -#> 14 50.2000 47.2508 -85.4608 185.8608 -#> 15 68.2960 66.0155 -121.2395 257.8315 -#> 16 44.7305 85.7855 -201.5661 291.0272 -#> 17 61.4670 76.4826 -158.1205 281.0544 -#> 18 50.4635 114.1011 -277.1294 378.0564 -#> 19 47.3392 129.9716 -325.8193 420.4976 -#> 20 28.6326 121.7077 -320.7994 378.0646 -#> 21 4.3445 133.7291 -379.6020 388.2910 +#> -10 -3.7266 3.2311 -12.3090 4.8557 +#> -9 2.6607 1.5006 -1.3252 6.6465 +#> -8 0.8290 2.2907 -5.2556 6.9136 +#> -7 5.2843 2.2970 -0.8168 11.3855 +#> -6 2.8555 1.6468 -1.5186 7.2296 +#> -5 1.3589 4.6708 -11.0476 13.7654 +#> -4 0.3294 3.4672 -8.8801 9.5388 +#> -3 -4.2227 4.1922 -15.3579 6.9125 +#> -2 -3.8447 2.9210 -11.6033 3.9138 +#> -1 -0.2234 3.4245 -9.3195 8.8727 +#> 0 -10.8156 9.3885 -35.7531 14.1218 +#> 1 -13.7998 15.7824 -55.7205 28.1209 +#> 2 -7.8432 11.4501 -38.2566 22.5703 +#> 3 -4.5541 12.2315 -37.0430 27.9347 +#> 4 -3.5368 12.3862 -36.4367 29.3631 +#> 5 8.5221 11.9115 -23.1167 40.1609 +#> 6 1.1140 16.5244 -42.7776 45.0056 +#> 7 6.6384 19.8568 -46.1044 59.3813 +#> 8 7.1288 28.8276 -69.4422 83.6998 +#> 9 10.8758 25.2019 -56.0647 77.8163 +#> 10 17.5057 32.3052 -68.3023 103.3137 +#> 11 40.8318 41.9009 -70.4639 152.1275 +#> 12 48.6134 41.2350 -60.9137 158.1404 +#> 13 52.4228 57.2359 -99.6052 204.4508 +#> 14 50.2000 71.5414 -139.8257 240.2257 +#> 15 68.2960 58.4186 -86.8735 223.4655 +#> 16 44.7305 77.2877 -160.5585 250.0195 +#> 17 61.4670 97.5688 -197.6919 320.6258 +#> 18 50.4635 96.0637 -204.6976 305.6246 +#> 19 47.3392 117.3176 -264.2759 358.9542 +#> 20 28.6326 159.2878 -394.4622 451.7275 +#> 21 4.3445 120.6413 -316.0987 324.7878 #> --- #> Signif. codes: `*' confidence band does not cover 0 ggpte(covid_res) + ylim(c(-1000,1000)) @@ -249,23 +255,26 @@ did_res_noif <- pte(yname="lemp", subset_fun=two_by_two_subset, attgt_fun=did_attgt_noif, #this is only diff. xformla=~lpop) +#> Warning in pte(yname = "lemp", gname = "first.treat", tname = "year", idname = "countyreal", : 'pte' is deprecated. +#> Use 'pte2' instead. +#> See help("Deprecated") summary(did_res_noif) #> #> Overall ATT: #> ATT Std. Error [ 95% Conf. Int.] -#> -0.0323 0.0107 -0.0533 -0.0113 * +#> -0.0323 0.012 -0.0558 -0.0087 * #> #> #> Dynamic Effects: #> Event Time Estimate Std. Error [95% Pointwise Conf. Band] -#> -3 0.0269 0.0151 -0.0026 0.0565 -#> -2 -0.0050 0.0115 -0.0275 0.0176 -#> -1 -0.0229 0.0153 -0.0529 0.0072 -#> 0 -0.0201 0.0102 -0.0400 -0.0002 * -#> 1 -0.0547 0.0172 -0.0885 -0.0210 * -#> 2 -0.1382 0.0391 -0.2148 -0.0616 * -#> 3 -0.1069 0.0379 -0.1813 -0.0325 * +#> -3 0.0269 0.0140 -0.0005 0.0544 +#> -2 -0.0050 0.0134 -0.0312 0.0213 +#> -1 -0.0229 0.0125 -0.0473 0.0016 +#> 0 -0.0201 0.0121 -0.0439 0.0036 +#> 1 -0.0547 0.0170 -0.0881 -0.0213 * +#> 2 -0.1382 0.0332 -0.2032 -0.0732 * +#> 3 -0.1069 0.0324 -0.1704 -0.0434 * #> --- #> Signif. codes: `*' confidence band does not cover 0 ggpte(did_res_noif) diff --git a/TODO b/TODO new file mode 100644 index 0000000..c669201 --- /dev/null +++ b/TODO @@ -0,0 +1,7 @@ +TODO: + +* make sure that the knots are the same across g and t +* (I think) in the function overall_weights, work with original groups and time periods instead of the adjusted ones, but just follow what `process_att_gt` does. +* I am just going to ignore the untreated group for conducting inference + * it is difficult to operationalize in code + * It doesn't matter for ACRT(d) at all; for ATT(d), it doesn't matter if we interpret estimates as being nonparametric. \ No newline at end of file diff --git a/man/aggte_obj.Rd b/man/aggte_obj.Rd new file mode 100644 index 0000000..d919633 --- /dev/null +++ b/man/aggte_obj.Rd @@ -0,0 +1,53 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/pte_aggte.R +\name{aggte_obj} +\alias{aggte_obj} +\title{aggte_obj} +\usage{ +aggte_obj( + overall.att = NULL, + overall.se = NULL, + type = "simple", + egt = NULL, + att.egt = NULL, + se.egt = NULL, + crit.val.egt = NULL, + inf.function = NULL, + min_e = NULL, + max_e = NULL, + balance_e = NULL, + call = NULL, + DIDparams = NULL +) +} +\arguments{ +\item{overall.att}{The estimated overall ATT} + +\item{overall.se}{Standard error for overall ATT} + +\item{egt}{Holds the length of exposure (for dynamic effects), the +group (for selective treatment timing), or the time period (for calendar +time effects)} + +\item{att.egt}{The ATT specific to egt} + +\item{se.egt}{The standard error specific to egt} + +\item{crit.val.egt}{A critical value for computing uniform confidence +bands for dynamic effects, selective treatment timing, or time period +effects.} + +\item{inf.function}{The influence function of the chosen aggregated parameters} + +\item{DIDparams}{A DIDparams object} +} +\value{ +an aggte_obj +} +\description{ +Objects of this class hold results on aggregated +group-time average treatment effects. This is derived from the AGGTEobj +class in the \code{did} package. + +An object for holding aggregated treatment effect parameters. +} diff --git a/man/crit_val_checks.Rd b/man/crit_val_checks.Rd new file mode 100644 index 0000000..d86b0f6 --- /dev/null +++ b/man/crit_val_checks.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/pte_aggte.R +\name{crit_val_checks} +\alias{crit_val_checks} +\title{crit_val_checks} +\usage{ +crit_val_checks(crit_val, alp = 0.05) +} +\arguments{ +\item{crit_val}{the critical value} + +\item{alp}{the significance level} +} +\value{ +a (possibly adjusted) critical value +} +\description{ +A function to perform sanity checks and possibly adjust a +a critical value to form a uniform confidence band +} diff --git a/man/dose_obj.Rd b/man/dose_obj.Rd new file mode 100644 index 0000000..89db99a --- /dev/null +++ b/man/dose_obj.Rd @@ -0,0 +1,55 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/process_dose_gt.R +\name{dose_obj} +\alias{dose_obj} +\title{dose_obj} +\usage{ +dose_obj( + dose, + att.d = NULL, + att.d_se = NULL, + att.d_crit.val = NULL, + att.d_inffunc = NULL, + acrt.d = NULL, + acrt.d_se = NULL, + acrt.d_crit.val = NULL, + acrt.d_inffunc = NULL, + pte_params = NULL, + call = NULL +) +} +\arguments{ +\item{dose}{vector containing the values of the dose used in estimation} + +\item{att.d}{estimates of ATT(d) for each value of \code{dose}} + +\item{att.d_se}{standard error of ATT(d) for each value of \code{dose}} + +\item{att.d_inffunc}{matrix containing the influence function from estimating +ATT(d)} + +\item{acrt.d}{estimates of ACRT(d) for each value of \code{dose}} + +\item{acrt.d_se}{standard error of ACRT(d) for each value of \code{dose}} + +\item{acrt.d_inffunc}{matrix containing the influence function from estimating +ACRT(d)} + +\item{pte_params}{a pte_params object containing other parameters passed to the function} + +\item{call}{the original call to the function for computing causal effect parameters +with a continuous treatment} + +\item{att.d_crt.val}{critical value to produce pointwise or uniform confidence +interval for ATT(d)} + +\item{acrt.d_crt.val}{critical value to produce pointwise or uniform confidence +interval for ACRT(d)} +} +\value{ +dose_obj +} +\description{ +Holds results from computing dose-specific treatment effects +with a continuous treatment +} diff --git a/man/figures/README-unnamed-chunk-3-1.png b/man/figures/README-unnamed-chunk-3-1.png index 44456b1..6dc0fb8 100644 Binary files a/man/figures/README-unnamed-chunk-3-1.png and b/man/figures/README-unnamed-chunk-3-1.png differ diff --git a/man/figures/README-unnamed-chunk-6-1.png b/man/figures/README-unnamed-chunk-6-1.png index 57924c4..e1180b5 100644 Binary files a/man/figures/README-unnamed-chunk-6-1.png and b/man/figures/README-unnamed-chunk-6-1.png differ diff --git a/man/figures/README-unnamed-chunk-8-1.png b/man/figures/README-unnamed-chunk-8-1.png index 3fc4579..9f0a1ab 100644 Binary files a/man/figures/README-unnamed-chunk-8-1.png and b/man/figures/README-unnamed-chunk-8-1.png differ diff --git a/man/ggpte_cont.Rd b/man/ggpte_cont.Rd new file mode 100644 index 0000000..c143787 --- /dev/null +++ b/man/ggpte_cont.Rd @@ -0,0 +1,17 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ggpte.R +\name{ggpte_cont} +\alias{ggpte_cont} +\title{ggpte_cont} +\usage{ +ggpte_cont(dose_obj, type = "att") +} +\arguments{ +\item{dose_obj}{a \code{dose_obj} that holds results with a continuous treatment} + +\item{type}{whether to plot ATT(d) or ACRT(d), defaults to \code{att} for +plotting ATT(d). For ACRT(d), use "acrt"} +} +\description{ +a function for plotting results in applications with a continuous treatment +} diff --git a/man/summary.aggte_obj.Rd b/man/summary.aggte_obj.Rd new file mode 100644 index 0000000..6bfe48b --- /dev/null +++ b/man/summary.aggte_obj.Rd @@ -0,0 +1,16 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/pte_aggte.R +\name{summary.aggte_obj} +\alias{summary.aggte_obj} +\title{Summary Aggregate Treatment Effect Parameter Objects} +\usage{ +\method{summary}{aggte_obj}(object, ...) +} +\arguments{ +\item{object}{an \code{aggte_obj} object} + +\item{...}{other arguments} +} +\description{ +A function to summarize aggregated treatment effect parameters. +} diff --git a/tests/testthat.R b/tests/testthat.R index 96a4188..66f8706 100644 --- a/tests/testthat.R +++ b/tests/testthat.R @@ -7,6 +7,6 @@ # * https://testthat.r-lib.org/articles/special-files.html library(testthat) -library(pte) +library(ptetools) -test_check("pte") +test_check("ptetool")