From 1bfd7bf44bb454ea68a640bc797f924ae4b02142 Mon Sep 17 00:00:00 2001 From: jonathanroth_MSFt Date: Thu, 22 Apr 2021 10:03:36 -0700 Subject: [PATCH 01/13] Added initial FRT code. need to optimize and modify the CS and SA commands --- R/compute_efficient_estimator_and_se.R | 62 +++++++++++++++++++++++--- 1 file changed, 55 insertions(+), 7 deletions(-) diff --git a/R/compute_efficient_estimator_and_se.R b/R/compute_efficient_estimator_and_se.R index e9b3ca7..8e0b268 100644 --- a/R/compute_efficient_estimator_and_se.R +++ b/R/compute_efficient_estimator_and_se.R @@ -80,12 +80,13 @@ balance_df <- function(df){ ##Check that (i,t) is a unique identifier - it_counts <- - df %>% - dplyr::group_by(i,t) %>% - dplyr::summarise(n = n()) + # it_counts <- + # df %>% + # dplyr::group_by(i,t) %>% + # dplyr::summarise(n = n()) - if(max(it_counts$n) > 1 ){ + #if(max(it_counts$n) > 1 ){ + if(anyDuplicated(df[c("i","t")]) > 0 ){ stop("There are multiple observations with the same (i,t) values. The panel should have a unique outcome for each (i,t) value.") } @@ -832,7 +833,8 @@ staggered <- function(df, FALSE), return_full_vcv = F, return_matrix_list = F, - use_last_treated_only = F){ + use_last_treated_only = F, + compute_fisher = F){ # If estimand is provided, force to be lower-case (allowing for non-case sensitive inputs) @@ -858,7 +860,8 @@ staggered <- function(df, beta = beta, use_DiD_A0 = use_DiD_A0, return_matrix_list = T, - use_last_treated_only = use_last_treated_only)) + use_last_treated_only = use_last_treated_only, + compute_fisher = compute_fisher)) resultsDF <- purrr::reduce(.x = purrr::map(.x = eventPlotResultsList, .f = ~ .x$resultsDF), .f = bind_rows) @@ -1010,6 +1013,51 @@ staggered <- function(df, se = se, se_neyman = se_neyman) + ## Do FRT, if specified + permuteTreatment <- function(df,seed){ + #This function takes a data.frame with columns i and g, and permutes the values of g assigned to i + + #Subset to first period so that we have one observation per unit + #XX this assumes balanced panel. Make sure you balance before doing this + first_period_df <- df %>% dplyr::filter(t == min(t)) %>% dplyr::select(i,g) + + #Draw a random permutation of the elements of first_period_df + set.seed(seed) + n = NROW(first_period_df) + randIndex <- + sample.int(n = n, + size = n, + replace = F) + + #Replace first_period_df$g with a permuted version based on randIndex + first_period_df$g <- first_period_df$g[randIndex] + + #Merge the new treatment assignments back with the original + df <- dplyr::left_join(df %>% dplyr::select(-g), + first_period_df, + by = c("i")) + + return(df) + } + + if(compute_fisher){ + FRTResults <- + purrr::map_dfr(.x = 1:50, + .f = ~ staggered::staggered(df = permuteTreatment(df, seed = .x), + estimand = "estimand", + A_theta_list = A_theta_list, + A_0_list = A_0_list, + eventTime = eventTime, + return_full_vcv = F, + return_matrix_list = F, + compute_fisher = F) %>% + mutate(seed = .x) + ) + + resultsDF$fisher_pval <- mean( abs(resultsDF$estimate/resultsDF$se) < abs(FRTResults$estimate/FRTResults$se) ) + resultsDF$fisher_pval_se_neyman <- mean( abs(resultsDF$estimate/resultsDF$se_neyman) < abs(FRTResults$estimate/FRTResults$se_neyman) ) + + } if(!return_matrix_list){ #If return_matrix_list is not specified, then we return results DF unless return_full_vcv =T if(!return_full_vcv){ From 553cf5d02e37bccdf0ec714ea2095113af19d2fb Mon Sep 17 00:00:00 2001 From: jonathanroth_MSFt Date: Thu, 22 Apr 2021 12:55:49 -0700 Subject: [PATCH 02/13] Have made some improvements for speed --- DESCRIPTION | 3 +- R/compute_efficient_estimator_and_se.R | 152 +++++++++++++++---------- 2 files changed, 93 insertions(+), 62 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 0419198..c84d130 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -18,7 +18,8 @@ Imports: Rcpp, magrittr, MASS, - stats + stats, + Rfast LinkingTo: Rcpp, RcppEigen diff --git a/R/compute_efficient_estimator_and_se.R b/R/compute_efficient_estimator_and_se.R index dda3202..a53bad3 100644 --- a/R/compute_efficient_estimator_and_se.R +++ b/R/compute_efficient_estimator_and_se.R @@ -2,27 +2,32 @@ #' @title Calculate group level summary statistics #' @description This function computes the mean-vector and covariance matrix of the outcomes for each cohort, where a cohort g is a group of units first treated in period g #' @param df A data frame containing panel data with the variables y (an outcome), i (an individual identifier), t (the period in which the outcome is observe), g (the period in which i is first treated, with Inf denoting never treated) +#' @param is_balanced If true, the df has previously been balanced so this does not need to be done internally. #' @return Y_bar_list A list of the means of the outcomes for each cohort g #' @return S_g_list A list of covariance matrices for the outcomes for each cohort g #' @return N_g_list A list of the number of observations for each cohort g #' @return g_list A list of when the cohorts were first treated #' @return t_list A list of the the time periods for the outcome. The vector of outcomes corresponds with this order. -compute_g_level_summaries <- function(df){ +compute_g_level_summaries <- function(df, is_balanced = TRUE){ #Balance the panel (and throw a warning if original panel is unbalanced) - df <- balance_df(df) - + if(!is_balanced){ + df <- balance_df(df) + } g_list <- sort(unique(df$g)) t_list <- sort(unique(df$t)) + + #Reshape so that Y_{it} is a column for all t + df <- df %>% + reshape2::dcast(i + g ~ t, value.var = "y") %>% + dplyr::select(-i) + compute_Ybar_Sbar_g <- function(g){ #Filter to observations in cohort g - dfg <- df %>% dplyr::filter(g == !!g) - - #Reshape so that Y_{it} is a column for all t - dfg <- dfg %>% - reshape2::dcast(i ~ t, value.var = "y") %>% - dplyr::select(-i) + dfg <- df %>% + dplyr::filter(g == !!g) %>% + dplyr::select(-g) #Order the columns ascending in time tVec <- as.numeric(colnames(dfg)) @@ -36,7 +41,7 @@ compute_g_level_summaries <- function(df){ #Compute means Ybar_g and covariance S_g Ybar_g <- base::colMeans(dfg) - S_g <- stats::var(dfg) + S_g <- Rfast::cova(dfg) return(list(Ybar_g = Ybar_g, @@ -44,6 +49,7 @@ compute_g_level_summaries <- function(df){ N_g = N_g)) } + resultsList <- purrr::map(.x = g_list, .f = compute_Ybar_Sbar_g) @@ -65,6 +71,7 @@ compute_g_level_summaries <- function(df){ } + balance_df <- function(df){ ## This function creates a balanced panel as needed for our analysis @@ -819,6 +826,61 @@ calculate_full_vcv <- function(eventPlotResultsList, resultsDF){ +processDF <- function(df, i, g, t, y){ + + #This function processes the df inputted to staggered (or staggered_cs/sa) + # It checks that the columns in the user-inputted values of i,g,t,y are actually in the data + # It also renames these columns to "i", "g", "t", "y" + + # Let's make sure we have columns with name i, t, y and g + colnames_df <- colnames(df) + if(!i %in% colnames_df){ + stop(paste0("There is no column ", i, " in the data. Thus, we are not able to find the unit identifier variable.")) + } + if(!t %in% colnames_df){ + stop(paste0("There is no column ", t, " in the data. Thus, we are not able to find the time identifier variable.")) + } + if(!g %in% colnames_df){ + stop(paste0("There is no column ", g, " in the data. Thus, we are not able to find the group identifier variable.")) + } + if(!y %in% colnames_df){ + stop(paste0("There is no column ", y, " in the data. Thus, we are not able to find the outcome variable.")) + } + + # Sanity checks + if(i %in% c("g", "t", "y" )){ + stop(paste0("Unit identifier cannot be labeled g, t, or y")) + } + + if(t %in% c("i","y", "g" )){ + stop(paste0("Time identifier cannot be labeled i, g, or y")) + } + + if(g %in% c("i", "t" ,"y" )){ + stop(paste0("Group identifier cannot be labeled i, t, or y")) + } + + + # Re-label i, t, g, y + if(i != "i"){ + df[,"i"] <- df[,i] + } + + if(t != "t"){ + df[, "t"] <- df[,t] + } + + if(g != "g"){ + df[, "g"] <- df[,g] + } + + if(y != "y"){ + df[, "y"] <- df[,y] + } + return(df) +} + + #' @useDynLib staggered #' @importFrom magrittr "%>%" #' @title Calculate the efficient adjusted estimator in staggered rollout designs @@ -837,6 +899,8 @@ calculate_full_vcv <- function(eventPlotResultsList, resultsDF){ #' @param return_full_vcv If this is true and estimand = "eventstudy", then the function returns a list containing the full variance-covariance matrix for the event-plot estimates in addition to the usual dataframe with the estimates #' @param return_matrix_list If true, the function returns a list of the A_0_list and A_theta_list matrices along with betastar. This is used for internal recursive calls to calculate the variance-covariance matrix, and will generally not be needed by the end-user. Default is False. #' @param use_last_treated_only If true, then A_0_list and A_theta_list are created to only make comparisons with the last treated cohorts (as suggested by Sun and Abraham), rather than using not-yet-treated units as comparisons. If set to TRUE (and use_DiD_A0 = TRUE), then beta=1 corresponds with the Sun and Abraham estimator. +#' @param compute_fisher If true, computes a Fisher Randomization Test using the studentized estimator. +#' @param skip_data_check If true, skips checks that the data is balanced and contains the colums i,t,g,y. Used in internal recursive calls to increase speed, but not recommended for end-user. #' @return resultsDF A data.frame containing: estimate (the point estimate), se (the standard error), and se_neyman (the Neyman standard error). If a vector-valued eventTime is provided, the data.frame contains multiple rows for each eventTime and an eventTime column. If return_full_vcv = TRUE and estimand = "eventstudy", the function returns a list containing resultsDF and the full variance covariance for the event-study estimates (vcv) as well as the Neyman version of the covariance matrix (vcv_neyman). (If return_matrix_list = TRUE, it likewise returns a list containing lists of matrices used in the vcv calculation.) #' @references #' \cite{Roth, Jonatahan, and Sant'Anna, Pedro H. C. (2021), @@ -898,57 +962,21 @@ staggered <- function(df, return_full_vcv = FALSE, return_matrix_list = FALSE, use_last_treated_only = FALSE, - compute_fisher = FALSE){ - - # Let's make sure we have columns with name i, t, y and g - colnames_df <- colnames(df) - if(!i %in% colnames_df){ - stop(paste0("There is no column ", i, " in the data. Thus, we are not able to find the unit identifier variable.")) - } - if(!t %in% colnames_df){ - stop(paste0("There is no column ", t, " in the data. Thus, we are not able to find the time identifier variable.")) - } - if(!g %in% colnames_df){ - stop(paste0("There is no column ", g, " in the data. Thus, we are not able to find the group identifier variable.")) - } - if(!y %in% colnames_df){ - stop(paste0("There is no column ", y, " in the data. Thus, we are not able to find the outcome variable.")) - } - - # Sanity checks - if(i %in% c("g", "t", "y" )){ - stop(paste0("Unit identifier cannot be labeled g, t, or y")) - } - - if(t %in% c("i","y", "g" )){ - stop(paste0("Time identifier cannot be labeled i, g, or y")) - } - - if(g %in% c("i", "t" ,"y" )){ - stop(paste0("Group identifier cannot be labeled i, t, or y")) - } - - - # Re-label i, t, g, y - if(i != "i"){ - df[,"i"] <- df[,i] - } + compute_fisher = FALSE, + skip_data_check = FALSE){ - if(t != "t"){ - df[, "t"] <- df[,t] + #Process the inputted df by checking the inputted columns and renaming to i,t,g,y, and balancing on (i,t) + #We skip this if skip_data_check = TRUE + if(!skip_data_check){ + df <- processDF(df, + i=i, + g=g, + t=t, + y=y) + #Balance the panel (and throw a warning if original panel is unbalanced) + df <- balance_df(df = df) } - if(g != "g"){ - df[, "g"] <- df[,g] - } - - if(y != "y"){ - df[, "y"] <- df[,y] - } - - - #Balance the panel (and throw a warning if original panel is unbalanced) - df <- balance_df(df) # If estimand is provided, force to be lower-case (allowing for non-case sensitive inputs) if(!is.null(estimand)){ @@ -974,7 +1002,8 @@ staggered <- function(df, use_DiD_A0 = use_DiD_A0, return_matrix_list = TRUE, use_last_treated_only = use_last_treated_only, - compute_fisher = compute_fisher)) + compute_fisher = compute_fisher, + skip_data_check = T)) resultsDF <- purrr::reduce(.x = purrr::map(.x = eventPlotResultsList, .f = ~ .x$resultsDF), .f = dplyr::bind_rows) @@ -997,7 +1026,7 @@ staggered <- function(df, return(resultsDF) } } - g_level_summaries <- compute_g_level_summaries(df) + g_level_summaries <- compute_g_level_summaries(df, is_balanced = TRUE) Ybar_g_list <- g_level_summaries$Ybar_g_List S_g_list <- g_level_summaries$S_g_List N_g_list <- g_level_summaries$N_g_List @@ -1163,7 +1192,8 @@ staggered <- function(df, eventTime = eventTime, return_full_vcv = F, return_matrix_list = F, - compute_fisher = F) %>% + compute_fisher = F, + skip_data_check = T) %>% mutate(seed = .x) ) From de65e1c99ad5318cf4ddcd0e50e284c4a783ae04 Mon Sep 17 00:00:00 2001 From: jonathanroth_MSFt Date: Thu, 22 Apr 2021 15:46:49 -0700 Subject: [PATCH 03/13] Sped up FRT calcs. Added FRT options for SA and CS. Standardized data processing for CS,SA,and staggered --- DESCRIPTION | 3 +- R/compute_efficient_estimator_and_se.R | 165 ++++++++--------------- README.Rmd | 2 +- README.md | 4 +- man/compute_g_level_summaries.Rd | 4 +- man/figures/README-unnamed-chunk-7-1.png | Bin 32830 -> 6619 bytes man/staggered.Rd | 11 +- man/staggered_cs.Rd | 11 +- man/staggered_sa.Rd | 11 +- 9 files changed, 96 insertions(+), 115 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index c84d130..170317f 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -19,7 +19,8 @@ Imports: magrittr, MASS, stats, - Rfast + Rfast, + tidyr LinkingTo: Rcpp, RcppEigen diff --git a/R/compute_efficient_estimator_and_se.R b/R/compute_efficient_estimator_and_se.R index a53bad3..2554a58 100644 --- a/R/compute_efficient_estimator_and_se.R +++ b/R/compute_efficient_estimator_and_se.R @@ -20,7 +20,8 @@ compute_g_level_summaries <- function(df, is_balanced = TRUE){ #Reshape so that Y_{it} is a column for all t df <- df %>% - reshape2::dcast(i + g ~ t, value.var = "y") %>% + tidyr::pivot_wider(id_cols = c("i","g"), names_from = "t", values_from = "y") %>% + ungroup() %>% dplyr::select(-i) compute_Ybar_Sbar_g <- function(g){ @@ -900,6 +901,7 @@ processDF <- function(df, i, g, t, y){ #' @param return_matrix_list If true, the function returns a list of the A_0_list and A_theta_list matrices along with betastar. This is used for internal recursive calls to calculate the variance-covariance matrix, and will generally not be needed by the end-user. Default is False. #' @param use_last_treated_only If true, then A_0_list and A_theta_list are created to only make comparisons with the last treated cohorts (as suggested by Sun and Abraham), rather than using not-yet-treated units as comparisons. If set to TRUE (and use_DiD_A0 = TRUE), then beta=1 corresponds with the Sun and Abraham estimator. #' @param compute_fisher If true, computes a Fisher Randomization Test using the studentized estimator. +#' @param num_fisher_permutations The number of permutations to use in the Fisher Randomization Test (if compute_fisher = TRUE). Default is 500. #' @param skip_data_check If true, skips checks that the data is balanced and contains the colums i,t,g,y. Used in internal recursive calls to increase speed, but not recommended for end-user. #' @return resultsDF A data.frame containing: estimate (the point estimate), se (the standard error), and se_neyman (the Neyman standard error). If a vector-valued eventTime is provided, the data.frame contains multiple rows for each eventTime and an eventTime column. If return_full_vcv = TRUE and estimand = "eventstudy", the function returns a list containing resultsDF and the full variance covariance for the event-study estimates (vcv) as well as the Neyman version of the covariance matrix (vcv_neyman). (If return_matrix_list = TRUE, it likewise returns a list containing lists of matrices used in the vcv calculation.) #' @references @@ -963,6 +965,7 @@ staggered <- function(df, return_matrix_list = FALSE, use_last_treated_only = FALSE, compute_fisher = FALSE, + num_fisher_permutations = 500, skip_data_check = FALSE){ #Process the inputted df by checking the inputted columns and renaming to i,t,g,y, and balancing on (i,t) @@ -1156,36 +1159,40 @@ staggered <- function(df, se_neyman = se_neyman) ## Do FRT, if specified - permuteTreatment <- function(df,seed){ + permuteTreatment <- function(df,i_g_table, seed){ #This function takes a data.frame with columns i and g, and permutes the values of g assigned to i - - #Subset to first period so that we have one observation per unit - #XX this assumes balanced panel. Make sure you balance before doing this - first_period_df <- df %>% dplyr::filter(t == min(t)) %>% dplyr::select(i,g) + # The input i_g_table has the unique combinations of (i,g) in df, and is calculated outside for speed improvements #Draw a random permutation of the elements of first_period_df set.seed(seed) - n = NROW(first_period_df) + n = NROW(i_g_table) randIndex <- sample.int(n = n, size = n, replace = F) #Replace first_period_df$g with a permuted version based on randIndex - first_period_df$g <- first_period_df$g[randIndex] + i_g_table$g <- i_g_table$g[randIndex] #Merge the new treatment assignments back with the original df <- dplyr::left_join(df %>% dplyr::select(-g), - first_period_df, + i_g_table, by = c("i")) return(df) } if(compute_fisher){ + + #Find unique pairs of (i,g). This will be used for computing the permutations + i_g_table <- df %>% + dplyr::filter(t == min(t)) %>% + dplyr::select(i,g) + + FRTResults <- - purrr::map_dfr(.x = 1:50, - .f = ~ staggered::staggered(df = permuteTreatment(df, seed = .x), + purrr::map_dfr(.x = 1:num_fisher_permutations, + .f = ~ staggered::staggered(df = permuteTreatment(df, i_g_table, seed = .x), estimand = "estimand", A_theta_list = A_theta_list, A_0_list = A_0_list, @@ -1247,6 +1254,10 @@ staggered <- function(df, #' @param eventTime If using estimand = "eventstudy", specify what eventTime you want the event-study parameter for. The default is 0, the period in which treatment occurs. If a vector is provided, estimates are returned for all the event-times in the vector. #' @param return_full_vcv If this is true and estimand = "eventstudy", then the function returns a list containing the full variance-covariance matrix for the event-plot estimates in addition to the usual dataframe with the estimates #' @param return_matrix_list If true, the function returns a list of the A_0_list and A_theta_list matrices along with betastar. This is used for internal recursive calls to calculate the variance-covariance matrix, and will generally not be needed by the end-user. Default is False. +#' @param compute_fisher If true, computes a Fisher Randomization Test using the studentized estimator. +#' @param num_fisher_permutations The number of permutations to use in the Fisher Randomization Test (if compute_fisher = TRUE). Default is 500. +#' @param skip_data_check If true, skips checks that the data is balanced and contains the colums i,t,g,y. Used in internal recursive calls to increase speed, but not recommended for end-user. + #' @return resultsDF A data.frame containing: estimate (the point estimate), se (the standard error), and se_neyman (the Neyman standard error). If a vector-valued eventTime is provided, the data.frame contains multiple rows for each eventTime and an eventTime column. If return_full_vcv = TRUE and estimand = "eventstudy", the function returns a list containing resultsDF and the full variance covariance for the event-study estimates (vcv) as well as the Neyman version of the covariance matrix (vcv_neyman). (If return_matrix_list = TRUE, it likewise returns a list containing lists of matrices used in the vcv calculation.) #' @references #' \cite{Callaway, Brantly, and Sant'Anna, Pedro H. C. (2020), @@ -1288,56 +1299,21 @@ staggered_cs <- function(df, A_0_list = NULL, eventTime = 0, return_full_vcv = FALSE, - return_matrix_list = FALSE){ - - # Let's make sure we have columns with name i, t, y and g - colnames_df <- colnames(df) - if(!i %in% colnames_df){ - stop(paste0("There is no column ", i, " in the data. Thus, we are not able to find the unit identifier variable.")) - } - if(!t %in% colnames_df){ - stop(paste0("There is no column ", t, " in the data. Thus, we are not able to find the time identifier variable.")) - } - if(!g %in% colnames_df){ - stop(paste0("There is no column ", g, " in the data. Thus, we are not able to find the group identifier variable.")) - } - if(!y %in% colnames_df){ - stop(paste0("There is no column ", y, " in the data. Thus, we are not able to find the outcome variable.")) - } - - # Sanity checks - if(i %in% c("g", "t", "y" )){ - stop(paste0("Unit identifier cannot be labeled g, t, or y")) - } - - if(t %in% c("i","y", "g" )){ - stop(paste0("Time identifier cannot be labeled i, g, or y")) - } - - if(g %in% c("i", "t" ,"y" )){ - stop(paste0("Group identifier cannot be labeled i, t, or y")) - } + return_matrix_list = FALSE, + compute_fisher = FALSE, + num_fisher_permutations = 500, + skip_data_check = FALSE){ - - # Re-label i, t, g, y - if(i != "i"){ - df[,"i"] <- df[,i] - } - - if(t != "t"){ - df[, "t"] <- df[,t] - } - - if(g != "g"){ - df[, "g"] <- df[,g] - } - - if(y != "y"){ - df[, "y"] <- df[,y] + if(!skip_data_check){ + df <- processDF(df, + i=i, + g=g, + t=t, + y=y) + #Balance the panel (and throw a warning if original panel is unbalanced) + df <- balance_df(df = df) } - #Balance the panel (and throw a warning if original panel is unbalanced) - df <- balance_df(df) #Drop units who has g= < min(t), since ATT(t,g) is not identified for these units TreatedBeforeMinT <- df$g <= min(df$t) @@ -1354,7 +1330,10 @@ staggered_cs <- function(df, use_DiD_A0 = TRUE, use_last_treated_only = FALSE, return_full_vcv = return_full_vcv, - return_matrix_list = return_matrix_list) + return_matrix_list = return_matrix_list, + compute_fisher = compute_fisher, + num_fisher_permutations = num_fisher_permutations, + skip_data_check = skip_data_check) return(results) @@ -1376,6 +1355,10 @@ staggered_cs <- function(df, #' @param eventTime If using estimand = "eventstudy", specify what eventTime you want the event-study parameter for. The default is 0, the period in which treatment occurs. If a vector is provided, estimates are returned for all the event-times in the vector. #' @param return_full_vcv If this is true and estimand = "eventstudy", then the function returns a list containing the full variance-covariance matrix for the event-plot estimates in addition to the usual dataframe with the estimates #' @param return_matrix_list If true, the function returns a list of the A_0_list and A_theta_list matrices along with betastar. This is used for internal recursive calls to calculate the variance-covariance matrix, and will generally not be needed by the end-user. Default is False. +#' @param compute_fisher If true, computes a Fisher Randomization Test using the studentized estimator. +#' @param num_fisher_permutations The number of permutations to use in the Fisher Randomization Test (if compute_fisher = TRUE). Default is 500. +#' @param skip_data_check If true, skips checks that the data is balanced and contains the colums i,t,g,y. Used in internal recursive calls to increase speed, but not recommended for end-user. + #' @return resultsDF A data.frame containing: estimate (the point estimate), se (the standard error), and se_neyman (the Neyman standard error). If a vector-valued eventTime is provided, the data.frame contains multiple rows for each eventTime and an eventTime column. If return_full_vcv = TRUE and estimand = "eventstudy", the function returns a list containing resultsDF and the full variance covariance for the event-study estimates (vcv) as well as the Neyman version of the covariance matrix (vcv_neyman). (If return_matrix_list = TRUE, it likewise returns a list containing lists of matrices used in the vcv calculation.) #' @references #' \cite{Sun, Liyang, and Abraham, Sarah (2020), @@ -1417,57 +1400,22 @@ staggered_sa <- function(df, A_0_list = NULL, eventTime = 0, return_full_vcv = FALSE, - return_matrix_list = FALSE){ - - # Let's make sure we have columns with name i, t, y and g - colnames_df <- colnames(df) - if(!i %in% colnames_df){ - stop(paste0("There is no column ", i, " in the data. Thus, we are not able to find the unit identifier variable.")) - } - if(!t %in% colnames_df){ - stop(paste0("There is no column ", t, " in the data. Thus, we are not able to find the time identifier variable.")) - } - if(!g %in% colnames_df){ - stop(paste0("There is no column ", g, " in the data. Thus, we are not able to find the group identifier variable.")) - } - if(!y %in% colnames_df){ - stop(paste0("There is no column ", y, " in the data. Thus, we are not able to find the outcome variable.")) - } - - # Sanity checks - if(i %in% c("g", "t", "y" )){ - stop(paste0("Unit identifier cannot be labeled g, t, or y")) - } - - if(t %in% c("i","y", "g" )){ - stop(paste0("Time identifier cannot be labeled i, g, or y")) - } - - if(g %in% c("i", "t" ,"y" )){ - stop(paste0("Group identifier cannot be labeled i, t, or y")) - } + return_matrix_list = FALSE, + compute_fisher = FALSE, + num_fisher_permutations = 500, + skip_data_check = FALSE){ - # Re-label i, t, g, y - if(i != "i"){ - df[,"i"] <- df[,i] - } - - if(t != "t"){ - df[, "t"] <- df[,t] - } - - if(g != "g"){ - df[, "g"] <- df[,g] - } - - if(y != "y"){ - df[, "y"] <- df[,y] + if(!skip_data_check){ + df <- processDF(df, + i=i, + g=g, + t=t, + y=y) + #Balance the panel (and throw a warning if original panel is unbalanced) + df <- balance_df(df = df) } - #Balance the panel (and throw a warning if original panel is unbalanced) - df <- balance_df(df) - #Drop units who has g= < min(t), since ATT(t,g) is not identified for these units TreatedBeforeMinT <- df$g <= min(df$t) if(sum(TreatedBeforeMinT) > 0 ){ @@ -1483,7 +1431,10 @@ staggered_sa <- function(df, use_DiD_A0 = TRUE, use_last_treated_only = TRUE, return_full_vcv = return_full_vcv, - return_matrix_list = return_matrix_list) + return_matrix_list = return_matrix_list, + compute_fisher = compute_fisher, + num_fisher_permutations = num_fisher_permutations, + skip_data_check = skip_data_check) return(results) diff --git a/README.Rmd b/README.Rmd index 874e2d6..1c1f12e 100644 --- a/README.Rmd +++ b/README.Rmd @@ -164,4 +164,4 @@ The syntax for staggered_stata is very similar to that for the staggered R packa ![Stata examples.](man/figures/Stata_screenshot.png) -The staggered, staggered_cs, and staggered_as commands all return a vector of coefficients and covariance matrix in ereturn list, and thus can be used with any post-estimation command in Stata (e.g. coefplot). +The staggered, staggered_cs, and staggered_as commands all return a vector of coefficients and covariance matrix in ereturn list, and thus can be used with any post-estimation command in Stata (e.g. coefplot). diff --git a/README.md b/README.md index 6cbbc7e..cc146b4 100644 --- a/README.md +++ b/README.md @@ -4,6 +4,7 @@ # staggered + The staggered R package computes the efficient estimator for settings @@ -45,7 +46,6 @@ for modifying and plotting the results. ``` r library(staggered) #load the staggered package library(dplyr) #load dplyr for data manipulation -#> Warning: package 'dplyr' was built under R version 4.0.2 #> #> Attaching package: 'dplyr' #> The following objects are masked from 'package:stats': @@ -55,7 +55,7 @@ library(dplyr) #load dplyr for data manipulation #> #> intersect, setdiff, setequal, union library(ggplot2) #load ggplot2 for plotting the results -#> Warning: package 'ggplot2' was built under R version 4.0.2 +#> Warning: package 'ggplot2' was built under R version 4.0.5 library(purrr) df <- staggered::pj_officer_level_balanced #load the officer data diff --git a/man/compute_g_level_summaries.Rd b/man/compute_g_level_summaries.Rd index dcf7cba..b68e2c9 100644 --- a/man/compute_g_level_summaries.Rd +++ b/man/compute_g_level_summaries.Rd @@ -4,10 +4,12 @@ \alias{compute_g_level_summaries} \title{Calculate group level summary statistics} \usage{ -compute_g_level_summaries(df) +compute_g_level_summaries(df, is_balanced = TRUE) } \arguments{ \item{df}{A data frame containing panel data with the variables y (an outcome), i (an individual identifier), t (the period in which the outcome is observe), g (the period in which i is first treated, with Inf denoting never treated)} + +\item{is_balanced}{If true, the df has previously been balanced so this does not need to be done internally.} } \value{ Y_bar_list A list of the means of the outcomes for each cohort g diff --git a/man/figures/README-unnamed-chunk-7-1.png b/man/figures/README-unnamed-chunk-7-1.png index c2b024ceacaf6a532ce1428e4b75c6f0e36e52f4..8915d7711b7f8d2db346d643287aeb4270f876dd 100644 GIT binary patch literal 6619 zcmcIJ3p~`@)_Z0!g;6@FIOr*lu+@jgBESoK#~NkCew@JpNq=MI_}N{wK<0 z48<57&r{w+qfpG1REiOj_pAB-Q|J4>d(XZ1yZ8IP-<{dB_J6Ov_S$Q&z5aWxwdXR) z#)2TJC3 zu$oz{w?ZxzaH&LYB3h_il#0ve^3fvX3i({2usOlHSt#Tu68VW#J{K*Ce3XZX1pNq6 zDirXA!nZ7+w?g6Q=qSn#9fJ-Kjeg4#qHKgB1}KS8DD-z8>O~c^`rKj13joN}iGG+S zU*;JAHh5YdG_m{X=2+I%2!3O5jK_|I$6WgH&J=xLx{^NE>quM53Cp|y$`8)5m1XhT z$4kNog~I(3uK;$_xmU_R?X3dY6Z!fkA6H&7@P<{cXX-8C#2>9IihyYVmbA?Y#2%D_ zuF){+>$HLSEpg;92gqGc`38ETcA5KVhfvZ(u!@AA1V^Ha-Nuj7yN2b}V1jQ(5@ zHe*tx@Vt_rY-`Joy=G3CxHDGiS`gmp;9$k}5iU7gZ`n~xem5-{KHG)p4J|+2d4zIS zYr=Z&yv*>HkQ3YPrC{5F_BOkvHhf6W55&;ySU#AQdA(rY&l}zry>gehO|eOylDe)} zmdvg&z`m8>o#ZRXx=0aDcf;fIu3-lEAG8$#N>2W%0opBoUCaB;8LG7j&Mpjb2f3P{ znt8_B<%u71KlUBEBF4MU%vo^g?mf8j=zw{c=|>;p!gPSYs@nA%*M=LIb$eTgcQ|R+ zmwZjfJJjYPP9C~#zOAva>p%=#FdD9;s==>3_o8yHQWGkD@oIMTyWzGxdqJpq(m@YP z^2SF(;fRC$lh>cdSN1TZZnIlr>Iuh+t^5-y4{Pn7Yp_jU2JeeEefCKwpFeTu)MVVq z+o?I*()9&b@~)A|$HuMCTZ`Q_hq{L^Oq3~d9z_}=ant6?NY2r^@iX_G4ifK}TTEVg z^(@~UI?Z<^{=mfP6%@9m1WmZQD|Eyf|ngOexgxrW}UMUGg(1AxjMpO~=jo@y^RQ zT|0AEB*@1<+xx$9$T_Qj(&K_cQWrg4J=FJI@sDS-?BwmU7CxD!)Jj{QI?+`c{S0Q$ zl}DVkFRlv|$^2&?>VuG0!Cg(h}SV$TWwh9Oy@8 z8vSG^O?K~Uut+8!1J=dpWaWpQcKTJ`9A`_Y9zrmbMu$91B^h0rPQ_uh5El*F+;jBP6a;{*=oU52m1c7v29pA z8RC(f83imxX_cDWHc4J#T09^$aa@J(h0CXd2*^X%oPvd7Bcwq(6YF#{vJYhJi&$K? z(1ilX$BM#3tm~OZknYa-nEBm^yjke|u6aSu+n5)6cfP}~SZeEXFWkTXd19dM1KQ`2 zP7>p!nc!*dYh3~7LF{vxNro1=Wk0&R`;HBI4gqp+*V&1%j;F?VeXixXIIJ(acQ`h-*@}?5X;Q`3@PyDtS^L*;RqtQ(o*yV8SJ}L1 zYuKe$yZSvXLb}c9h{6MT%RjkA{i#)&2+%8n=zo z{f3oD`R4ocOP0sOC6Z;V_#@Z)f*kv)1vtGEz0?AA{XS}ILdwaeR!o7hbZqOB1~n5u zc8T{Lx1?lAGL?sMDs8tD(|fQhJPEb%TQ*U_M;0xg+Se?Tp#$j!jqeRoT45^xsCKy8 zDYgQJuq)?I&N>UFo=%gWwm!W&tdsL@B1ZF$7!5U>Ee@W0T^FZyCMlmihzrZ_y{^gm zo`_&1O?Z_k3-#>1ts*8}x&`n3e+boYv5G3#O~b^uEroFte|@}15CoxuD|4~Tc+KQFI3UmqBkIr0nfiZyt}4J=I?rQX8%v_!!Z&C>{nEt9v1nSbo4mI;sqof68uOLNs?3JiqS+yd9A|@qN5GCL z50N*(l4V$?B#!ipJs_U~nq88JrZg1ggt3dlFiBc)!FNF3Q=B9(ftdaq20fPMWv(yH zntI*YUHWV}7zB8sQ&Gi+`#ViO4q6TlyUot)1*t$3XH=Yxn2;%f-Ar0rZu!;lRl@cT z9Fs1RZ!1)sor5>fxf*A&Ijy%e<`>U6DCM=gwq!x=x9j(7TI!Bvi$N}jG3yU z$#kLkE)O#a_y3hfbXXWa;3R`1x!MD7H*u2hzpGM~Kr%8lMq6>TZ?xc$Qta0%7OUr4 zQP{><5|$*g6&5_g$Uzp;%-r7@mVBME?JEOqEqKp;EVKOo-|6`%fC)ra>gtH$ zT$YEb&x-TZ?7>WmJd}9?Phbui1>TT{C{CEJV0#eYAx?Vli0Qg54`sOr;D*ISR$kE| z!1642Qii-r6zwD+3`gx-XtWM=1=!D-zjEkis{p#;mk2IL^Ro*v}u@;I<#%A%z%rq^qB zUJeSdk30tg!0)y^u5L6r3#)KCJlDz~GN3v*K=Gu}{`Q2Vn97L+XH3`Z9n5euraB~C z<=Up%+RF0>K&i+S-#16D)G* z^4j?6pn=K(sz4vnE+^;u;WE4!hr*~Hs?(k4HxD42fh021umS6tCqclF09!`_Vu>R? z4JRZ~1s(#PJ!^qDIvwCR0Z+=6+ag7XA=E#T;xov(E}#!#gTC@YMT|P!+FT$DW0_a} zqA^T$42aMagD`fG{feqHmIUMz<)Km3o}N3>Kfl@mV`+2f*6-1R1GAGi*X$D2U;3GF zp_gv_deU{qgfZR1S1KOZ?L?@#`J;u@22LV*D=MJG*N0?D=AcZ%o5JWn;JtCA9a_%O z^JLaV(1mPHn3cfzcYAIE`m{@>?R0d*IpUCVxG@1im0?=qnOC`(;;5}p`yE%2t{LY5 zE5@L~#|(c<|8lE^!2#H9W6fl7Xj;qXtakz_xoKotw02)R&}HV#5NFOBbfIreCEl*~ z4!TZp-ub=g>WVh&^Y2hVUoRfj$}kPnohL>L8_5;pD@scIq0Uh{y)1{^vx3Vrx-0*X+9FkPqZm)wAV9vrxP%9|t&Cn)g3 z&y?Z(v51S(7D@R9xag=e@IT1}v?Vy(T8-hp<|U36QQ`zDDgt^AK*v%RR{8Hod*+e`Zu0`UUx~4gD!;zozj6agKqr&goW9hY&uT;y)$O zdk_+sp)9Z31Bm5(DwcE+fMlEpk3&yCkLeUO`6LzG1_C1E+6tXsaZ6{ zlo>JG_iMA$+xL4P4cll$BYiHu)-#8UgR_@X)4GeMQ!v`ToxVi8IIqUNOsjfE#b2IR znY%Ki8Sbh*Ftx@9lU3TjO@~Dp@zzOdb)NRD)fMbi`Vh$-NZfqu);BOwSaNuYfZ^F< z6vv(!*ao4G;pG8u%-6vg>)!|}i;0$=tTa(ep|nK&zKMO-TTaz2(_`60(pg4Gc?>^}QJi~G1*&W-cn0{MmF;pX2A0kbS07O0 zd)Q}1zp!_Pl&o^jrlqbyU%Y#p&8j?zg!}A0aHHzhhS9AGmRZzpk4xU-Je9v&D;pyH zXUeKPByIzmCFnEMLeIuvGh(wO(hn+@4;jHeC$FGKlMIiL_XLZCV3-CF{fuMPuBv`F zuTm$4(~cEqs|vZqR?qb%_KW2gSuRm3H&o=H_)^rTWf;_}iTIr|lYWHpN8Pj}675K^ z!tjy+C%DHL^j)pR0OZj77(DY{?Jc*J=f&?pg{3_xll1*KhkGx2o`idrt;K*o{gTOE zcu$9>B(J+wdL8LVyfVqw%n;5%W1;Au>Y_LNj2l*S6!+G}F#URDF@3hqtS4(CGFy(9 z_+*%9;P9N5*XyC@?g5;jA5jmu9W%W6MHs^UuC*9_$C(W`@uX6eNzfWJn2;ExmpTb| zi=2dq$ct*PMm@g%D%aHRkEKTc--zViN?IkHBJppYdQ}{%vC?!^9`)H8alEY zr`J1sW%2P-2I+=W5&kI&EJ8(10oxz zX`$vQ%9c=2g0A=$+p7zQV!w3#6M$x3!X+^!xvfJMa*y*Olva%=hPdCxi<8jA>#u<= zEYOrj-jscr=3f9ZvOGl$ZQ!fGU!?w@K?>|7(#a`pG{&++^S^%uaISulvgaBuWXg?u z!Nsrw^*RUEApK5Fgg@p|`?(29yiAAuf~J!y)nU_S^S*&9&==(Iw)fh06?OK*DQ4dm zQswz+UG1^YCz>Ua3(G~$iE1H>^p+kDn81~mOZ7A=`mR|F+`}lTK|N6d`;EBX@vrHW zW+TxDnn;*UWzA2kV|5u4NZq|zIo_FC-EqBMymCodF|2rZSG=S5y+AY#qo0~VBlbT+ z@Aaf_Mm;ycZtcoAX4$<()Pnv9iP)~)bjnu}nls^{08MSY^~G3d^26<3TWEvk9rkAT zi1)GoxRZmF37|CL#{Hq`-7UYb)&BsE^hY9pO`HWQOVja%>MR;Mz9dP96iN4 z_u+EWrrKDa6yfY4$Ne&Q za869%^T@_ey)QG@FBwmay?Te-Ct|uH7?ogJu*SB?TOt_G zE=M)ZY8o4H0}_{zt5r|5D*Ev=FJ97~>pR^boKlthrTzE&>i^Wom&>a+4lC*2bA9KI*}Nen5vF@tp4zL3q{{#&D- z09d=k2`vfzU9L>fm3%%fuOTu=K6g_00<^Ke4NG!p#e}Rp8Cd>geeXd&tS6zTkK9D* zKJkz4xXvCv_+o1`@AmQ#HMJYQ6TQn*#~ui*Mz&~Ry=1)=8{EgMT|VXCizgFkp3uI( hy@&eS%l{c+k-$sh+kHLgzYYLcn%f*KGX3#S{{b_V`$7N! literal 32830 zcmcG$bySpF^fyjQiJ+8%N(%!bN=lb>cY}0;bazTiHzP=Q=YUEmFmy>P-6-9>XT0}% z?{}?tt>3%W?|uJZtzl-K^E_vtefIwB&;Fcws~{(ii%p7+goK1EDe+1P2?+&?goOMZ z^A>o=BXWTa2?_0vg@}lPm58{Aowc2V@*4vq6LAw;69)?;C2?USq?b`qs=DS6RSEgi zn`#(o1_H#rUDdjNXhIiXu_g^qu{wFIGbk4Cb<`oL-_OGPSv{EgjwMWA z^eKIDwNN`9aX_?jPehGi|Rj;G_fUjTg#@KcccHqrCf1g~N#@GQKvU*{HJldhY;U(Yjte>(nsd*L(W>>8;-#q#>= z*81>s;rtuwpF_o;_smtPTW|Zgo0m8bVVS<69(y6lEEyMg`j9Gq`S}CNT#7il9%bW5 z<6&;8JM3SQs2q(Cnuf>Rk*SrkIU`yl)?J0Vug#28zSjj8PdMwjn@aD7r+GM5v8mHFnYYZP$->n(72KFsBxEi*1 zvz>OQ_yM1sF@GMDs62&jl9-JnoBO_dFBF z>{~A{>3B}G>bdZ~m!b|`nD=q`nJvNgkV^UU-Hw&q=sFD-P+pD)c6?HoCmk#y-#XHE zv!$I!`sM}GyScdX@IF8tTR>`r-k{#ly5V#gP4s~P@zO*?(o|L!i2=OEL_!X=Ktcnr zkij1*_(Q}Jlz-kqfkILLd5!!X@!>uFnOP(xK_tmnLMm>^+nE@%gk5J}Lg&-NunRSe za0?Zm`g(Q*h*fsZSLVF1yj>r5OWOMGc2#FNa&C?xr`j0(1ig!j!d;7zJja_WS(fCM zC0o0b`573jWzUO$ucbjb`)pu)yGEHTkJsoqz8?}A1rjo*AQB3G*jvG+lo84YzyJCB zeGL?Ps0!M@@AzwALH)YJ(>wpU2z*aLun6@(<5b-fM2;fjd@BTo{CgCrpAyD@E-<7W zVw)Jdm-Y>IGH<2lxx=o%aXe*4KH6ZnP@9cG?lp9qoPSm(lY9K~2cGZqi`|xr9<$^- zLm$aU^JQ5X%Nq?*jn4Fc8dfvt`(7@Q-;jK~7;m&U8ZFo58LhKYuerItNFP_%=j-F^ z;f&!sPI2EXDuDNME?6JjO)g!__HA2nT6{})OU~h2rGe?^AS~1Mq4)kBLBCS8ldlg8 z1X79AG(KG(H?-pt$eOo1uEJlKWjd=3zj?&t!ltfcFRJM4^LXBKkI%7>sUgL2fa^q)R z-dwNN4D#djbJ-c1RO+<4Rk?4KHNY7wS`GEQk5tpYL~+i4%kbDNGL2%NpclbjiOEUe z$a#On)__azTOM{71sT(ve8_|I%Olr}V_{k6tqG>O&Ek^PtCrO?tC5t}ThCS?Mn#YL zJU-mv*{&QTvul0V`8y?_+zUu~{ftR=|#{ zWm@sRjDb6DjU{Q@TFktt^G??D*l-wGv7&R>dDgTX7x$7mS>wS|;-?z1BQ?2~wHBk@ ztCPz7JKi;>{YqI%S|KXEFW&^ir z25>oNUz61}Eco#32xRL{SNRJJyj+S@0S-lWzwaghu?;iza}t(&oJN#wKP>$|x=X;&7wxpXXS z1QTR6f2)8pF@27n6E6x94M)^i(r0 zs~PG(DU85m+fb(=`vfT3j$Y%@$MUW$&r&SCcNm>ef@J1 zcW{+q2Ln86a%Owwj%x#T`<7pyjk@8s)VaDW(ys1Tu(3yGSomE264y3mZv7c>w%fAD zlG%IAQ&=jcux}Bjx^i=MvS8MuQ24CQW@f^zS29U|j*v0pxwB`}S5ytQAXt`cMZ>%s z_FGvu+J)5hF8}Yd+7z)KHMeO zO-Mp3uGg#;uCw-`HaH=-9~iO`y%O0QaN3=nnsr~xG;P_bXQj(@it~0^Nvg>FYL?vB z*HvQGcKmb{rf~&PyH2&YuXt0r^n0Z|9>V4L_R&Yk1G?i3JDQMQsZk!=YrC-R zO*b;R+UV0UU>n7a`)wTkP82L5tC<3x)`{m2}B$RVlY0i zVM`xPIBh`6;aO3_)u!FK3qE{MhZ9mS#)^oGwu{}G6V7?JT$t1Z*a8E`f<*C>zk1vb zP?gKH>Z)P)v2XhHuUU>}9S5(z3}5WALm_r@*^fBJ2$66XETum3oy^7`=aLWdUsziv zs;X5ve9ivm!A$G97xI0tdFz?y&9FGRNw70tURIpp&bXCBFQni@pIcs{niS-Ye3b;A z1p}_S*wy(_-=;St&z0I0lfST(G|5q!M(1n~WDAHqa<%(#Wl6lUx_a_H?k5t)2&gj; zj%^lM`{SXI%E3oD4__Dt^{&?xMF+yu)UO)DxIDOtrf`{Szw`{OKia<1Y z&i?wNA1~KUcH(B^=38S2167PxytLW^=zamJ@#W7|*L*cScaw`tV-2&bLw+xta}OVr z{W^lPZ3lNqslMkE*Hnt50)?|1x{&P~j7LriQqGb#rj$z%Xc2sV;rZ-%N%3Q2bvxc* zzP$ORorVQ|n>9!?WDk`pv(>HW(Spm2jdIp!^A6Ybp-=S@7V`5kvm8_sf4Ra>{x>7l z&0g!7?;9_tqhvJCaQz+by~do+-5?xN%oWRP)hloHE&s9K_E4Eb=Z(L2>^O9TKRZ+J zTyIaf*laKh(%&+Q)0<9c7aHI)Ws>alz|uBpHSTzKjwkQ^(AmLSLF;2IvKyO*xuks$ z0cE|xAn&T8PX}#IqXv$?Zm_ue_oK5i>`(P#)Hm02W9scCZ>4%wuXg8p3ZScWXu~{? zsBA#0pM)I0;&}=8 z`D%sFRv4gNy;SCN$32GAMJXs+3>U)8=8%j@eP7s>Un~k$u;2k&CWsdNN^E;W+RH{; zDIqs`Ud8;U5GVs9?DA&s=6Z0#eO!JV;=OL&cD8B}+zYwiwd{KKwUy!ek;^yAQOl+U z(vfkDxx3;I4d(uK3L>jPLieF_YIHXJGiJJ0v-j?0+IBE!{(@E)gQl%qljX1P>gO}O zkEiN&R9%nYTT{uxJ4^#I8=4N?kZq;t@jFd;=HAA}J{oKFcZH^tbasU{LnklN_QjG) zx6g+rE@&iLe(ZM;p{zXdL)pCRj9qxa0C#nkipBm&)?I~H#n~G7EJlaN4!!9pcfL_( zsq9Vtv8`}x!^jka-h8E8egi0`vQl^7_j1}X_tm)OR+z#B52@g6mQG{R>1T2MiR`28 z7dwOskhH~I0hq+ejGe&DQ26jif7HyA1Z6?iHR9lRl2J{oIVPmA&bl|kjyCMBlGtjy z(qhck5E#!+Z9MMgj>PLL$+H2rP4qZ8SWgdsGluI;6Gc3OU8<)?8W3Da>%elGP6lkD zju3r}Os>YlQzTYX(UUPrdH(0?>1X=7y{du{iK0WznQ~)m;xEu^RgQ!8hSb3}Okz84 zR8DPryr=ZJ=#j$1#W;&j-D4l@!jVa;Quw@F*L%E2!HL3*ON5bz``UvKZVvC}yB$($ z!3)c!EL@)G^SrZ4?qihjB(A-dK<;PqrPE#=|0QtwvxBwF+c$J2t4TsSyitsA?6FBb z&ivu68RYLC26=wPHgWqKk+a6xTzgzKQJB;9%MN&20uF0_w`@2wE_Emsg8_wxv#DDr7bxAgsGb?Vr_1l*I z^?q{cg&}8}(X#68C5W}A9<0JJj|gnIVA0mA^vL zYeq);^V4n9>Tww23ks*%TC;vQXa zPWuTM;*;jdC*T3G2;eX@7@xvas2++ z_tL%wK%r>?n`Bz|K+=93W%rB3;K^LlLR$>k(Lw#DDZA>PCz?I(jK6kZt5cSY#lqdp z8CeFYI;PXpUzICyL@jn;5Q_RN9qf%+)H%g{AY=6S?(3$v!aSWi_c1$ zb`*pL(H{oXnhDK_(Z18ZK8sU*%qKcNpxwWgHTdkQtIf)C2RH@7xRDVBEj<51FclGw(262Xj5S-oeK8)F_C8g(rYS8%@cH*HD#YwZZMK`MG z8N^MX^5jMAulWA1Jd%uqy^-PZ^aBeb0f@rcvG?KXM-va`Q*9})@*a4+<0ut*1;)2$ z)WTjQ1Fgbsrwp*ZRTR<_*de)}Gl+dj3;oAS&l7?0gPTG#Y*_{X_e12>OtenkI? zDyz4DK0WJQVh;b1v4f_m`@Zev=O%Q=v*~CrDFf5l;>pzBC`VPrFw%ZKo^h=n%m`{YH+6|3#SbW#oG{BJ8+vWO$Kx-UJ^_wK!2 zVn&lq&WEG*Vuz`-?@({V?kLjgl_=k#woRSt{L9*KohR#|Qu7kLeH1&IL-mRl%SpcP zo$`#W4%~NWVl{6Kv?nOrAs4Pd1Pq#5%!ay14vXulCJS39LH*j=YsNK`H0y*jmz+CC zPsn+U|C1Wu%X)yBcs7X|`h(NP=7+4#pu++2{i*IRHZznfBDis+T{6_7TeYKdH86_J zyx~3>-{43+mCc4BudJ*Cv7pW_`0F1ZO_so8MmW$7o@Mb>92UJWO)O`}LlmB~~!6MlFaFkzGa zI`t+o-n;450f$oq--r2}ic)%KC9Pcyxg?A|v8U|bG_F?U00ACpMXybypQQ6gNjpjRAmnW5)Zv$KIc}&fy z$?qYXu@MPt?{DO?XTNRs;Strc$&?v#q0TP-CIo$twvOCSQRwcS){r{rRzUOMNDZBR zo{Gwo@cW+cjc19%!#5w1U6|WOK;%%So={0#V5p^XccDn&>8e{if`SO5@_hOJ} z?qbCaPsS7*u0oQMlwTY_8u(b}pP_eae=zY(W{=FPP6$Xo`wSsY-e26g|N5X60Gk;+ z%p)I*`EphvV^zG^e66tedJZ0+q^gF7Bj!=c4tZ-0ogosM5{knwlDjtRy-AZ^z3|$hck|uR2V<+|LY3 z*it{SC2B2LH{rra>`|^%p2y)~YSD0m2|yb9w&X9>S<&>OAi8)#G=02uW8Ko|B7%~LCO6329_tUQa#($(o}DfP2b zqosx!@n_7zOcPl>^0DX2g_lHavBPST(@h;!W2%FV`-?|nX6RO!u8HrCPvLCy7!GPP z@5gti=~T3HHy%L+>b?LvM$5EvrQQ-}1c7vv1!|>3r*Bma%TWtWM?Lj#@)u0SWK-hR z6YYE@A;_nc7FvY$(+5V${Y&DS%`2^fPMBbFd?rWv`DLVsG<#$Tv-X*M02V|1K>@RK zJ38h(zfZU4@I|?f+FoeL`qiT1(Mj}Xjn5bk=|tfd6&x%F5}Kf>FUm;U4>jyYF9S*@ z0(*;_3OC&qy_N06`3XKeLx&(Lz3kh#b&MnzA-z8!;ijz(r?P`(#CN!Iv)=Iea{Eib zPP!tqeEX5d*f85&ax{3mW>rl+d=7QYMO{_$H9Yg3a*G8#vc{4aeAALQzpygwrp1(V zzt1x1$@Pe{Q&gh7R-CZwiEB9!iFL`1E&V0}B?Rneo&xo4K#}_2+#B-Z_HH$Qsmycd zd?zzxskCVV%U=e)R<2rxgAzad;|>3$2zm$WZ(23X&x4P>LJ`Nk_DkNX zvgM2PRewu4j$XPWG{B91+y_BTHw%jcL6Zs-3{6;QZRaGyH?KxQmd-oI=#ZHh#@jMN zA38{QqfkuZ(+eY^*J4s+Q2tHVf{S{k{k#SH`r^O+hg=QE7er=4!8H31XS)v=Ste$A z!@z&J-v8}jVEk6zx43_3TR|u$GT{)F-mI{I)xKbr0nc79Il=GS(~_$R6g;iqvE^C8??m5(E%HFTSuio9*D_}b^EejdBc zXx5mHd#K6AKaqt?F_sy}lWcB~ax5lohXk=UM=9-;_ z1XoeNMEFaUf6b&XsVhzI zj0TS+-qZd|4@XAO%E6A_>!Vn8#fXbkh`wZSjHG*6^5tfpuH8<&%ge&hHMMXv*s?0N zIgGbP(k^DMrAnbM0%a5SnwUvTB}{&5XZ%}Ns@@bU{&QhfUUrnDICX=f&{l?QC8bT; zisFu=IDK!n32D8#h$_8lBZvw3Jc^}vOqwil*-VZ4wom@?7x62H7_$GEavPY^NL()9 zk13&Tz>*WLXYrnY+zOuu_{O()@f3fqE3^aK`SBkQ>XZXW(q)+j^REsAaM{gt zBy4l@WH-TUm82@klP%q0!5 z&=7Y+PvHJ*qrAtu+|s{ciQ(OsSbu|t=txk2=Nv=*GZUo$OHBa(6o~B=H9P2e#oeBK-Q

;WQ|LnP$(3+d!%vNI$& z@yB2B{+*bAIKH{oTa{2;GE>kjxHnp&&g?iOgps-w%`>nK=)!!!pd$H+(RCJa|RyR z&f|F7f)%#k%a}e4f7tR0~uT!7xcaoZ}lo)`)jtk z*(dbSP}#z2m5at4a60c%>ZArM=&5op6^ok?3phFTT2+P6I(@*xrJ$8QYJ9w3B8Tif zTM~tBI2O&MU0=$u1>2tn9XsvvrU;F|W!UtkxkvhZOVi%GWuG{{d!qWd(ex5UZO&0{ z{pet4xt|Ll=;7&{gHLrDht%$270eIX&NZ9#Oyxn@e{Tt)!b;0yJyRtmm9v00c+T=;=h@zpYn6(vPlf&!$*r3m8q7N&oaggw zQ1hGX<@s$sz}pm%91PiqpsVb}PcM5sXQ&X3ARhNpJF-`lw^4~ZD8@q;Ah>uq1mb@h zQ^V{C@Ac`DTDiKu_qh4}SH4Q{Co(CocN)GK;;_1GYaGnEE==jw0_rw(9_D$m(v@(p zn^XAVpU?)?9wKml5r-ko#U%&8bQD7J1;vzx;v5qwfgP-g7>W z!{>-r(vkYQKEwdKE12&;&c%ktkmDznieSKpf4Y8&el*H8<8{zK%9mTcJ=0LJF6Pe4?;D2{8 zeWu{Hqo6cNV(fpyM;(Y$wZu~Ic>jbue+DWL!$pIG3;%o#N{k7Th%h7$L_sJj zDXhdN`m^cAP(Xy!q8Hu&7v>a)fv`wKM4SUMM4gDBxYxYvg16~xqkUaG!~t|d1t}=N zHGa^ah{y`fKy%_afeHm8Ixf8rnCuaK=t67|;Zzl6FtPqbIm=Km=uQGH*T2z8RS3+2 zAsmbH&)Nv?05nyWSxrOt?<^R=EYw6_P}2M{ec>Z8dYy98G6L28WhoVO@I=M*PlAH} zSguY4jPBAKm+`NSia`R5fQ>hIVt$}{)|3$2aG-qk^VjT zzhZxJ0GN&{A%XC}S#+8bjLu^sDTj@p!eNufwHQVXZ@;x`;TA}_*juLV;uk6f5wTG1z9T{!}YJ3`M)9nW(`T@_x8v!dG!*LkQ&@05`Xva zPtoq)XETrG0^5Or=5CH41uuBcpKwg%KUSw0LgX4;KXHXW^F#9g2`1RVN{#l%pnig= zzziHndT{}5C-a`GKRfWu3_aqo(KY`8x?iQiQ!yuP){GI0Ft?rNXjZ-!I>L(S{Kx)d zMBvWN47Myde{Rhn54BhkEi6ToP3K_GMP#EVaK}RK(vyK|lRgFp*N1Nz zE22o*i#I(Ugea$Y;6IsnIMuB-erw0(T;Yd3l)>fbMVZ)F1p=Gte5<<{-9f?(Jd3k&?S?llE6tzEHW%Na2@mWLZN|-#hIN`%SO6{hl z_1UHxO@O~sxg1TKa^(cbED7R1D)$<rF>V=!Q{rDWL_FghzyUV!W*%do=8f_(a z8lEE>0blHjK%FG>UjEwtdAxb}LjJ6*+e1@oil&KgLcBQJ0CPU(j&0T5`&ASvb;p+8vV&xd;p3z1};<@x)froTM_Y5-h^ ztUG4tpLGX&ixNQ95VR22|2SC4T|s1<`vSru|Kg+|8UWYs-Jk!n)sfIh#lQwNxWvT$ z<1!TZap18fUc3DZR`Jd61J4TiMa}jfZIIQ3L z7|>Yjdwcrp3pp7iXY+0sBDi9LN|_Cl_7d^B&=FPpLHT6O7xDUuSgInfxeJ9%n{7W6l^ziohn2ej@3ay zNSrJt%QGwd#2?xRRq#aPhoyZ)hz#U-e{~kdQ(m(mih9&e3PpIfz(Q%`Vz^f-x_>Y^ z$b9x_+|HZEwv9h-jh1TGv{pBPHUPIFNOLL~A3NU?^#ADGTADWU)8Rl(U~_)7)$lM1)DEQ152XEDe|3}f1I3k{n!dMl z?pXpt?l7F(m>^r`hLA>;vJ^+f(#dC8bGvMtFKtcL&|$AT7+3vUgC$af8#w`G*lL{g zvlpP<%>ugnv0S^2OK>3ivPt8#%X$KWlp9a$!H@`EYIhXb%-b4J9<+0v%sTdm-KB;j zluU90S58~{jPvszTTC_MYTBvH?Uy4u_FY73in8L}|E*v;@xyeH5i$P^h*Y9Ho|qtH z2cyl-R=MMvR^?5mq;Sx6vAL2sOBlw#MD+Q5W5_rMjXZxQ*a4$m zJZ1moQGans`DW`>$?Ss9#ZGOehl;wUN#dcmihxp4LKGr_YWnqg&fq!23aC}biD;5(1`*RT*bQJm4?Be4AHYI{`vE)`Cozw3a$$WQ5!e zF<%{nEfeUDOs~y!y^U4w*jOO_V)CWsVNRHK0HbZ*UC~P*Agi|Sqz^yM3@B-wL$xI+ z{oyNgu#u&#oBO(AO$YZ4p}iTg{?Mgq2;|uFQvmuz#B>M!GIO2mo1>GgHBSrFOieu?75AKR;EmCzAWPL-k$=2~swk5*IHk{!D)07(4Xm+wtgDF=O zDi39dlna~<_lf1l&fFriek~ULF!ff;0QXi3&u$~NH1Zrk?fpPiu^LxZZ9ZT=CU&4^ zqC%ngfzxqShOFJgH&9KCjN-dMPZziIG!2s%E z9`U^_sP(;!M8TY3i5-o-*WH0}I8DM>&tagU_<+5!#%xfI5lX73PCJJ43O72E5Upz`^7Bjd z&r_oIiMI!guU%YlVa+d{o!iMs%}m`mGaj?qo=xg2+(RAT9zzVZpf=Ub?Zp~Y1_-WxuMr8+qnY~%Xg?fHEtC9$i**# zt4gojW2-2}VA&^UlPO|B$0E59J4}nnL4Oq^aNU?F@5nyh2Ni$>Ogl)5X){JW;+gO| znRBgQvjN1RONUVt(AT)EZWLs5oMVr<$l1023Yh|085`!V2WA-tOLxs5_0)CT;rWm# zAN*I!lKKuZB^xO5O%?V|e{)Enl?}*;L>*L&J8Xa>sqR?WOdU&-VVQQY_!>*8#xN1%)B8TQGH0F z?^T7{A%H*u23Y0|cbF2!?k3#-4GQSV{ln>@N>Az%<3if!d@_6ER?V8T6z5KMW+o3e zM!(#dZ0V_~R%b3Il$kmW0bo}Zsh?M(xFr^6uRF6XPk<7~I z)QTSaTUFmP8GttTHD~k`69!w&EB`Hxi0JS^H66ZoANe1A55`~sbXjfAspP{yBaV@S z40_J!=y3ZFLi_*wB^?tqC!Y!)Wj^d-XW}NFWZv-8w*-IPFcwi_Bx9c7OWa4!Bz3v` z^+o(iVl)i4lB7))ksbCMI6`0gcSVqpWgec6Da#X`=%Df|2wFk$drad> zLEwHtdQ^-Nh@yp;C-e9Z==>ju=C4X|3>8e$gWhBxaHcfTAA_(6N@bd#X;#zN?^EI< zKovc}w_=jXjfg>(@$U=xv~EmRWKUTXWR~_J5}yxH-lmF^d$WC^Ye1}pWwPMtbU!H< z!l^)N26Q_M0&k>p*k~oaM&u5JgK43Ae)q>PZL&uNriBt?4WVj9 z@Ux=Bz6apIK@!if==#kXuuFm!dJ>T+14?X3nxBWG@)zOkp2yY)fc-i>|8y7FQIH8y zKnE+);DOhRSnIV#T#l=RiH^Mq5Yeme$cS&$A&}|uy;nZq#tfLX)>9)0fFo+6Yg0EK ze*Wn=Jx`9RHi|&6T&w&O}-Ilx0bo*VT6fz5It=2n|3(q+Q3em1J$# zOjkYYJ9?kO-uMYVz|3%8O*IAj$v!e=z`ykmy%&)MSzy1}2*CFU5r^ST&3(lfV*IAjJT<`(`QB=rs=Ehu0L*0UH%T5N7 zI$|HOgMzs>Kcc~h)fGg9^clN$K4w1>J-1J{AF$9=kR<_ZCm4F?5j;^{&jKL}N*@C| z-#T%ypJV3f=7+cB2%8upfNnlL{s%y2G)HDPfg-pcNOmtzRRy5{@u^bK zQvF@#TgCvao0z=O>V5#PH0XeflIR*cIsl(f2fi{`7+Rs%-^*BG-R!!zo9Q={jppuss@M9Gm@6*j>z>aE zkA78Wfzg-fHlH{D1OW z;<;TF+5|jtzuDj)Hb?LfF)y=+lGqErH!%0>vo-IfHbe;u%51NzvFhAoWN#q(ha*ag znx?Dc8L~2;c0dwGYWwT-ozdC-UAue$eKJ|g9>RB9Kt{l{*pZU$bKWD9%92;uu5}EB z5e;&0Jx%ZWO)gjWJ*&&m-3Npdjo(M07vI5D;d4czf!@cb{DGGCM5j%o(`w{{M(zx+Qa?2v+ydHE-Di>P8BW<`esea3F zUHGlhE5v|6Ea~M`OKZXf>63PY5f^9L+ZWG&@f=Hi$y*`OJG&jPewK)JuR*^#@r|Q( zj$sAUe0k#{wj((0AojT+p*ozQB9XqG!a~gm=EM>TDS19vT5Yr^3fqvr1zR* zgfav{dvz$ogtm1*SsgiZjf}qY2C8)qgqZ!)vggsJNwf2o+7Mr4;6Q=ehK7Bz zPI>^K*sdw;xdj$!a&z%=!xmN)?ZsczR-ERR7Fg2IxGHcM7~7|pZjrM2sM`va$7%12 z-)$HH5T$D`>mA@mUQxC4Y$_iK&}Tp2QU1n3DdFC|B5uM!E!HdE6xFLDHXA?mq)_LG zB>3ry>EzJqtRqqF04E#FAw>N_Pml${XV6Wp#lr2;ubOiUu)x55NAbe|-F$bp>76u} z@3mKYh+6(&7GKNLjZ9Fttb(Otyq~Prk<3^IAk5U+2cf0-&+nx2P_moK-!mX)5|BEu z^u1TT*Msm^%iJwOS3=y@MSak|psB`5W-NeOg9?jKgV=`PNv_qJ_^OYDX z+qV=P#p#Nj?XN$fq)XS~Dt7Mi%dNpv3Fv*i(KbuZq18b)1SrH7Pv;LuT*!YGi>z>z z!$|sQYO$2M$D#}2mUA!j=*;=#?%}DEbkt%QOfOa12#ZWL;#CR%dX7X6Xs1{G~75YDt+3=;Lu|H|Q@l*Y{Py*^t7M+xN8 zG*lsYfE6U3@eP`IKGL+7WrgU2+~GfcZB^t6aj;!+jjN&KYM3%5Yq2Zxuwb*-zi?=J{1nbR?dh2#&NA`Q_CxW{83e0D$gSuxZd--b>8hEFgYYOC=ek{6-*TxI@ z2cQBxJC3w!@Y||rZ=;mNCk8PU3C$3&xTL7YaT-}4nSIatR+t>E15SIeKHH#RBC6eO zVE5G?y6oB+aY-6T2sg1akXo=Hrr4vjx^)JRdK#8J`+W8$&|P6xy6T9k-wa-TNvtI;$fdsaMVukx)n8!1KPm%ANX|G(%-p zIsUWd+1Gzh;xi(9l?rakK7^!k*A;o5?7EsIZ#I4f zR7C&b`Y^Q}H^`Ci-^GO&`JhbWoC2Y1@pWV@=+aU!13=WOfq;url$D@5NCoo!A&x(TGOX*_yzO#k0omC?#Q=I(s5`d}GJ?*=~z(eT3{ws8F(fubU8Au~@?J z6(AKp15ko1g3}H%q9RJNu1l%~$Xj~8)sy+ZaM?kienDzE$svDUR+AbC$FnhS4*Jdi z`R^lG7(8)9;5ygRP(YWMp8otHCf;e9hB*Gg3tDof2mPFbwX2+H6!@`zt8!LpJ)n#; z>}zc@S`!>hWA8(>z-B%7V5Z=yt30>6nea9c;y6HpO2z&J)6qkjHnhrg zo~W{IUY&7Zv=e+viRAwjWm;XGF5(JMY?{Z-pmkR@$M*_D@HK^kEV#0ys5bVm_RD7U zM+AWZXlLwI?JvX&-Y0YDbM}6M_ppML1lI<$?bcjyAFyU>^~n zeVGQb+vxK7L^2J;<#;~Ldx;=d;_2k; zXQ1YTFwpCIp_nTCPXwbx*_3fV92W0SfGA_TOD4eXY|@(ozFiHDvx_f5=>HMBkOLh+ zys{nC64rWAJ5^#zYzI{AnSHOgWjpIv8fce4Ca`6B(HX=jjXBzCX<*Zi{)j8TP=Jw#q@_dQ) zOBiUI{vYK9u8a}q6dflt&8Bt~63lG)Xn1xc-}(#EVtIiNOhM_VhUQIE_1=?KI`@ag z31F*BVaLzc`}6xJQ(%1#AZG}ZwFM^^#bylWKqFg` zPX6-lUGw3b;3jY~2#z=ew4yXrWJnqe-+rhyKtPNyLE$e5&TM3FgEKy>k*ovH0o_*L zAL1c;Wa16nHl6SI8@0ct2o2^d&|*Vt`}rBC?cKtg4&1igm%!!?eqqZfDfS4L;0V=ENU@>k?LWW_ zk-Ag7EJu~Ed#e$;e;O^K2cF|T3JyW(F^6N}nb7hvdAwu->yf@mczt9Aw{d~M4+*ID@Ly#c_J_#~z z4aSSEgfR;H?-&T6ErNB*0JVc1B_Y754JnBDX-z;9Thn~O-dPA6o6vQg3+YuH{h?fx zV4sd=C7{}$Z3EgwvTIFAtPOytOamZwaNbNbHrgT$zVYMl;f&?xYj?UKBo>I4 z=2!8L>d50~{eY*>e>5?i6$-fK4?Pck#=0BQZhM2%n=27yJTZQRa{cedZIM|x=bYjZ zWyH>tYv5OnbxZn3cZWD=r^6+2kXWL!@NCrcxkm0|c1jA4emI`@w40G5ky5|lVrBn; z{|Ip1{x1;|%H(n7sA#`#PBa6tp{`>apxG$pwqm(0P~_Bof%Rwc@i3BekGc*)g3t;J^Bm zzF~)`Ae`AoqZbhAl|d9tw4ZPjQ2dcq=jo}g5#Q3i613L~42i zrr%GStYZkybdk&ufF7Ut!OX(-#51Ml?dm>Z*3n9s+{#c9i;vjurewm6Om5jLHUU+~ zy;?zrQ9|ta4Z6*R&EO$PqOM3{d)le66#z)A#Uk$;v|XQDUEEbFHRJ3#6XN*I10ifK zj|famfwhZsnT?g81TG++Oce#Q<|8OD;`t7YnJU@&$S#Y0tcM~yLd&r48jG0>t8W0 z+SEHG5ZZ6dGWXAVSJXX)!nOHF0jJCpEF@k+J8W*y4-P&hGm_Ai3{a?@z5YW=32LDD z*OBk&H-U~;$|Kj=;fDYxy8L=qZ#*L9Ru+qASsS3oxZw$eLV@NfT70Dd4#Sm<*=HFA zcG*%9^*yY8uF-y}T5Yax2Y~7-$7eXf_WKqRTnDYo`j?9&Vtn3;sk?J6WK0}Ndu8bi zm2_-p5=clq>4@I~7?EQ9jeKSSvwvLi(84M_DDv1i(P2h`+GO4iy?qgJ`V;5)d-HgG zxmj}xuX4Uj)yi{lUfWPlq-6_L6ro9qH2yRT7$xh=>5!grbt}buUv5K1a&vItRVmdF z9xZl|INn#rXBBvVMG^KxC^FhRVVKumb_UWgd{DX{6&@8`ez}@%NA*_UrXcbdePB+$ z{$Ylc1;OF$EpJnBph_fnjh*2bk=Ju|Kru9p<}HCo&fog5v1U<1D;+;yC-H6+ONa$Q zHH-cavZ&I5bWY}KkD^Rkkday>7F~CWiM=s)zGIC(Hmk?7e;DEjC@IqRs{vF%1T!Sv z%OiPDD32Rq^$Z%lzCD~P{+)qNRK$~>8(kC`A)?|un&|Z~LkOsL=DdM`svmHcJ7pA` z^AE0)lx9FDJHQ}ApZ#13GjmvH#4tSRKzmb_2g5wxT$LWudd{Rg;G`Q^{MfQ!J&B>> zGuMQVy}9k1rhO=eymN`pFOyzyyhSWBCsN;D{Jmc>@6iw9;O#}KU~0G{{;EkniGO5P z6a6;{+a7WvGW+CNFnK+WbEGh;Yo+O+=O!gNM~Uo*844pb{8ge-c!8%|p4+CC zNO;)R@S(JfEjOgw`;4a@<;Q9=cZ-hi{*uQrb`w#xtNmtimiW%QwaAdPn3s!~u4GfH zM3y^sll3?awX{-AQN&IbdS&%GiWOEOQOD2ww;=;03jj@-osh!U@~f4+?VV|}vt?Rt zF4}KcdG=asn!q+?E$iFY)-Oq`exLc!WzKoBx-<1bPaqz@xo6-WsBa{Wk;lUV`7s~6 zefNhO_Z@4Ki>6gPJUms{Vc+oJh`@S3Khx%YLAv9vZO$!K_40C-eeDAD$=*grwS8xo za_qq6JbkMem!5_DNDI2pcQ_;>9p1sL@Qp(#yx!%dqUYmI^IB2zamDdIK&ww|uSYdn z{IY85#}l}G&9UI+=m=zo(`>zu+C3_3IVd4E0o>jn*oc%dA7a$)7)+y)Ji}6tn%LyK zuL|1+OzALv{&c1T8DwxXOgSA2^9)RRzEj1GTNfg-cMPx@yTl>xey*I50D%VbxIz~? ziO%aazg4k9F)ijF6k-L8Me47hnQx}UzaC*XzpPy=3~|eNM{aXQlY0=!+Q1Bs6|of4 z)hXO3!JB-I=JL%TGzYuz@_$AUV#*?Ijgn#Z80opso1fs9qd8EhScL5u*~)5^s8l(yXFnAs*n@0`xVbqsBHZaplb~XmKxY1My?; z704~=lp-M-D=G)v2(_8xto@hYKZ7yA>fE71DmMNYs`a4aC%5lRKRFZOlc)DD$s7S> zoO!=Om+i%Z6I%acWEcJ1{No??-@g)(wbHVS&M?#LXwj?QmlUwMs=)Hod{diP)Epgj z1`AnonvAAHl{${JPM0G+6xr1kaBSHopF1+<3LEHXlz+l9pC1CNZF*?%20*BE(BSMH zC{x6%U7axmtn@i8l!J0TbW@ zi#=ww7)96n4`fpY%`QJLAH-w|ASbpE+gX(FM9=3SA>{qY-q20F-`qDLia*g1sC*`$ zp(w=|?+$H&u&E<89faTuVkJDD>PEkS9~rT6O&y|=0wXAJBbSX14X{o_ z1T)D}&S~)_kjJwWwRFq-=Q+Q<4m{5jMLcI39e~x$-_F)-Z!}C!)2zy?=;|%BAITDZ z2ByoqJLICkZ7irR<4u5Bi>-IS%x=e8BUEl@4HVVm@`9d}mY(c8;zfMa&ZXrum`&%7 z>RNZk4vK;!&!)NiH7E62ggC>V$)X7VR=nY%_DfCbf71)krh9_0(2RGnoW!fqQVV&ut zKa1YVz98t7uX_0gQAR85p`qY@!iih8E;OT9_+QPvWmr{f^e?JQP(XzR5(0vB3P?yv zECG?0MnFnI8VLafDHSD@knW{)OM|q6(vlK_64HXwjb}`(eg5}+Joi3#fAHDcXYIM> zn)99G9phK?3vA3r-_JBlf(Q}@+(a&Rv59xjip0ZWgGW@8Zw@T8Cvks{GV&axt(ae& zpQSQ&FQD_B^QXM$iy^^NiWO+KqsT}WzMw3-^Uzw4n56>ezB>kr8CT`cw8#B0ae*oNKP(wmwU-j;@0~HC$)wj_yaFLyHx6Tukzyf#R5} zrWen-%y{+Hv4xe)J$wkWX`D>C?y=UFt3QQmn-nvF0l7h8G&RbQ?-YyD4`q$xoV?A9 zquJhbGW4!(lZgHy*uWax5zV_{%aAT{g@wml4#$Us8Jo$|s1#zdD&-Il@e022Io07G zVDKC2zshyd4)6_xk1TuZ0tgbXx>$bIvO+S&2+)-Jt-64{XNaPbJl1v9VgPhrdFlBB zEqVA0{oQIjz-;bu_96HQckQlYveQ^X*ZC1pnJ6TG^S1W|ItPc6foMEEgzb}USwh9k z%fN3Rtw0Lb`bp0^*9=KoLG;ztRZ5kOR`K4HsuXF zT6jX^(y%nov<;m5mRZq!vPQUf06*b97^@un#R8$ib(`iL(juCrvd%H-fDJLEc^s5p zY%&<>uiLLy_%A)~+IwT&5$+}e2^EGH-|K%^f>GaWIFs?uFU)q&@4B zm|C+EFi>1@pXy4d+UlMuC8X=p)`EnYd_3AOBG?!baZqZqCBmXVMT`ab~cQ-=VK+3tHUq6(Fm3eNK?rNtuVicg-!%) zAEtVZ+ZSJ&RMd2f7STLk6egIo?*Mz9Ah(f*|(ACKAd-;5=lRK3EjW7#8A3ubPfijM@RCf*`nk zBL)G~Xor1WwQEb%S?`1?@}!ZVMA7D;;x7nHEpfe9%_Ff5XkkEPyxSGfF=_W2k*;iW^pZ2HVfbhw(- zMjvAM-y>?e^>cR<$VT#wBssqiS?3$@E>w`z$^9mrVTfc(@OFFDXE|dI0U>1pYoy22 zF&9DJReT|ED#bE2X3A#~zIEMZlbBI9I>BYAna8UfDfY0!kQt>x4hsn56i9#+&a#Ea z6_ZyP-W~I;_*gIk$0=binetu-8b8tP4p2l?&I`0)nVkaevABPLiG}{p=m$6ZB~Iyx zi~h>!L!ww<7kFCpr&=VI8cwaJl%0EW`u6>e!|mbHfkErTPFE}S6gcmg6_-VNS|@QG z?4!R9gDJ$bpwA8>sR7938_JPAQPi=#`(4?w)##f+yhA^GC~2+=-CO>;%IiVp*3T|C zButI$y}y4+h!V<94W^ykaJ<7FWSyYd`((;J;RQhGLRwnPbL-@u;guBrk*)veG!OYG zK@qk{yW)8OqyM@bhy_&@(|RdX`0oIB2nPj#|H-EW!GEq9eCZ?|`Yzdh5fVm1Jk1B* z>#zE1O5|hi8sHyIi4jAI=5L2yq2;yCEzw?nB+Y=FB=IS>EN`ocEx6XYyzhc&tkYBy zVoT7oz*O?j=M|R8PI051laZfsp zH~3o$8dSQJn7>KqsyYh~=7*|D0FU8b{=OIl$Y=RC4so46LinFkM&n%qP~;@>)jqdA z-An|AxxNpR&wxOzpe>cW8yGB#m*zFI*p@*Q$nSpWQYR=@C|yg0Bz|fVDll5?N~kq1 zGF6x{xq~97HL>P&*1k-mtaN>B@M(el^Umd3+1im$@?aD>qFTl5Me3kP*hr&8P0aG; z((Ys1Q{nwg;t$nI@$(F2BaGe~uJPg9K3LA27R)%8b$pq+(P?tvN*_zPt1W5ae1%5! ziJqfJ?`bET-QZT+!{_iit$V%tslrwJE&Ki(?hXCbC=Lwdg*jF+RC0r}<6ZTgS@j9e zDEfvy`Z zj3}fkrUjyTAmF8S=Xn#;sNv_50p-k3#)J6w&)+*auVZ7{TZ22sm@3UvMCbXoMLodF zAvg;ZU3_z!jp|!$W%!?1R_1@+>)0H|E!}4bEIX1DCz9{3uUbE6<5qrFD}xzDv|;v`_qn ztZIX*TQ0~AD#+t&l!7$LEN*l_Ccf?&QuYn1qqh;pKAewy#pLlY7wI1 z{Idoz`0x-ZddniI5XkDKN3BeMC{uAbXV1lSWF0Lh3I|DvLN6cGA&eIKF6@GB08yB8 zSN)LUyq#}Y-8N{RWP9sFq5;`59C(}w*fLhz?gQSyCzk04FqL|fPTjfD`U&w4+Xg1rO~y$XbzeOT3Ys#U{S=gBoAx)>;b@)&{! zL|#=U+7h)0d0$qK329>Vfg*QJeiQ})fh@@_gdb^RPmjN`G2K3GxZuKu?K6DP zM0D`Tk};3ENhJGE+S?h3J_F9N8Z%<)Z)3jT~=FW4JmagUR4k^oik zFyFmIwEKu&F7uP&gJKebheFscukD^Bf|i%G-JsPVyYv2qCq6)m+%+*#4_*UkwYiI;Y@q z_4HXBNBkT}!JDc-c2)q=Kkr)v78ql}#Mw;O$YWg!{iK37uwk4hMy3REd--Z@F?t0iyEu6<46Im}Bfm%V8_9YdLD@i4|B zs3FfO*i!WE%zil17G$`iTmg!4Y%ymEhqkn*U%ZE>7KY5P2&e5yNlyA;kc73)3aMa~aAdEgOOyl@3f&n3_OEpzwrK3Ir~&mPexX0Q!f84!ISFh$=Bf z7w-zF@6Pbt|3!=A;2_lE_WFYU`6)&2iiSIY?_=FnkDYu`1K`xh9b`&4{gFr^Ph!vV z^Ahy^1h@71W~i3%sHhopjWF%0eTirQqUof@bW3fYqCRHXoqqsw_Id*-`+L$}7a)G#O+3t9K!z z1CHqwBEK10>2-dZ8s35MMaYwMGI9#Ou2}y80JK+i9k{N&s=E%2*%ZK@t%cJ^Q4BDU zN|NI8T6U-WKaa|U{jAEK9J`^*olJkfHllUtt0#B|CKZD2T0Q#JXtE+ghZI6_kgu7% z2s12*I#N~HxAifd}KcWhIkX6<)>E%t3jd)T@^xpn~6UUq#7Ar^(Gj5-b< z$1_D+^kqeWi!yf>^nIs7W;`xZEwEOFdiB}L$<>CaKE>&Y@E$_{DU!PGc9_+OlXHke2 zi%5Jm#ACWVu>;AawQEn=yHx4=-B9v;6FUwG2FoCh1^R=3h$LAj-_T|09mqY*fxr8n z!xOMH@jeb<@=r&Tp;JYq3R3_X2*IMvOSoMID%MyeEQbLYt5fZO9GaC5S#hpFR%k4; zOX7&=2&?S3JalvD8#Os7cXXbp*th8$?=BsqF{KC^h8gM>gwH-TLaW2x|+HJ^A&i#qTR?5xIAOYUzbcY&^)%WMU7B8AlW5C-;1 z*6G4(YW54ywLH<}p7W;DNYQBrVtUGsvDO8~;Fy=E)Ww|NVwEz4_`ib{gkiHCfB}z0xsk&!?K<*=ssEF39X4BVeYKogX$p# zt`h+acKCVAq(9llU60Xc>_76C%rFu{We3$_O{mmI=SWxiLpiYst$q)*jX5XWnk^-B zhYhcI*qYhhru#6Y!P}Fxo6)&9EI43%XHnB$^bdMeH}SC{^dTDDCT5R2TOVW^*&u`&AAlVAuGbwtOnDz6w7+B3cJKgQ^qZ=M0kB^sB#{fi-0bl%GR4-?KbTM5l?nOy+MmYTj>XqD-Ht*crC| z>V!cAPfmclfhrmF{A1_1=p09O5$liE=~Mv%ztZc^yFMy_uJ*!Qo(?0gRWBRTn9Qf` z4JCSRTV@YiTueuB6dRC-$?EML7!lMaQWN`;P|Mu=Lv_2H{@%jxM6_No(f`WOK)t{h z1Cq7IB&ECx(Sv|o7bwus4ezN_7QfBvW@-1L`G&qa|FR_5SHIpbkuqE%FSdR6e{v&~ z)W3?8i8io__62xinoEVjy4-6Me1f zgdR4v#ZO()H29FQ2;81^%4UYX-rJ213&+>sMxgJiQXRRNEWaTdr{zSG;tdSJuL+t> z3BwoxDqc$yqQpR;VIZ?Hl&F3OKQ8k=`Et65d`ydZnq;mqVt9{p(^8Y5Ihprs%JzAK zBj7k=zDoG9iSa0El1Wt!w>X}Xp8L;S{^nGtrqDCY6u(|`&Rc`w8{_vqUp>UIBcsv_ z**`+2XXs8Vvjjo5%bK6Ig^K{!8j?1EK4vW`!k%$GuwksNtPv6u%llzpe?bp ziG*?dH4BEw!hk|D=0LM{a2b}MA$*o3!v4q7eaPoE3+UKha>-*np=^7po5kgZi!mTUM!zRQB)9c3C!5z@Ak7E4^I+l7~$+%U-l@VJ1}`5%P2e)ldh70*sjN@Cm@#4 z3_OlA^T*x~As&5qGlihvYt5hj2Z$vrjEU0SB%lxJZfFXiTY^)07JPyIH<~yXgo)xd zK*VomFy9Q^xti`KcI?DyeB$*(>`jD>}TYB zjU}L*EPq@lMwa;J8{gv*j{ps$@`_$B2k=gCpO8irC+N#jwxutytnW*IL|`#^^T3j? zCinY1gUVP)4v2y1<|xjVhIfpO{1FNzjMS-tD4gc2G?IWXN%1Hcl1N5=Azr z=i;{aD*uW4ov!--CWj$;O|3!Ucz+AbA)H;11#kcRlus5tRwTCZv1c|5{{?O~0*hx9 z=LdhMC-olp0kL=w5PX@2h5CJ4GKX{L#-;E6slIfo9&q^V{I-zOf5C6f3vlBP5Qg9A z?;%tl%#g*~Ehqf{o&-h^YoIs2uy z7KTxI53~4TDLS0Pf-i{ zdSguBEB?R!6x68bliUZVtV)j(mA3$-Bv*Bg3GzT>l(9F&vBSweE?QjMrZh;wqc%%yBQ?%4vECQ z1)|4W5Jqs%iT}=ws#Bkq4&~NanEDPz%w2IHXb^2_8P2ELjPQQ~S?KGG zRwQ7LMug)Fo$P{(ri~uVe3v^f1WYS9wnEt%f3_^*u<&foTy-mKw~gHW=lg$#aYbI3 zz#ndb$F`94pJjaKPH;NMA|NaC@!dML=R0~3&p6@>aJB9Yn+qQ+HhUzQVlz*q;f@2_ z<2Cmfj*w)xmTa`JZ|b3-*g_WqDP=PbR(`1eg@SI){rUzSmqItR1Ee)=5Z7N%?*)Z* zwMH&norO!JOFY0>6ov+&S@zF+u8`lAi3A3_=yuzGr!urU!v-LuIU;e6s`(bM?G)ie zXeWu3X?O;B5aZw%0bUjd3A1boA0b`h8EL<;5?GD} zfEMB_<9RQFvS}MbS!)H5)CUkNj7_KB8b43SrY|)=Dd-3%(liv+EEV}_}=+w2^dJ<84wJtU)_0~x#I>e&JjRd)a{MX3|Xbr&tEkn1Hr+o zb$y!xU&@OvImjzPY=~w7$5ZCxHDtmzps#&Hz!?pEZ8WGM1dQ%jR;s^AIRcylQWfSNc$3!#L3e4< zdjh%+&BwoH4My3guU6>!SwGl8?8t1Q+StsW}W@86*#y~*E;ue9-} zfz!t1p?0blY=q|oAN>h=0%z+CHeNXRXuiP3$D^DEsFYu1#la7F3SXK)COr6`s10N` zm?fd}`1=K_<~Y$u^vrQpT>hLc+p{{*Z&=?E6HXh0otU~KLKFF&~a=@fB>+XpG zF9RD^O_f>Y*YCqbkj#LQEkYq^pic&77b^}f(|;G^0upR?lpyrnpZ)ESj};PZrkv`V zbFjG_Z0bi5>oA7?Iqttp_P;x8kl){X2zEv)HkrcT*8-MLmOotm3lZdnQonc4Lq1eU zB-jyOIr2Z345E9GF(Iv6O8ei@k05>Te|_fY3>NrReVTerzh@FSx|KI(BnkIRh8iOh z4Kol>2pw)g8ioElM-p`jjZDpvuUdn_MYNNnL2IS3pE!rK~c2t#UY7PVdRw>ZKzt4qs7R{4J@<}lG2P#a?%{iv}IxVE; zsr6o{b`S4o8_yG{yg$oYAp_5(t;uLR4?1@< zcMQe3YiI^dXE!ZF{Zbptsg4Bw1WG*iwV!3KUwGBEY!5v#x`t`tpmWCW4n=2SywRin;wVOyzow_WUw8}rY zrisNIu5kUHI9%;6G#8%}$00G}p& zCT8Ylz?>OfF+rxepxXQR`0-PEq!=kKgEt*dMJ|DW#v_oX7LK6y^zfK5F+Zs2+TaTKG@A02=E@#K! z_O_R+bNv2Qm^wZ9S^xh%m9KKgnS=sw>16Kj@7&Hbs4RtqF%?k|l%x z$vJtFQu@yl7$(IYmymL(wjX=&ZMZ_CC4!SxD@)VAQ%&UQ7Zj^XJUT&bK9%8efZSwD zMD>&5ci1|hx=mr*?gnk8@uD+XAiOwF4zd)!0!Ct7@(`paR(x;+MXDnAr{u^ zz393a&0GSZuWZnH@|#t9?3990;ELVnDt0I!p(m*Y<$!6Hs+ozGvM)aWK8S`t=`9(D zCf>eW?_HHx0eL{d%sY(4pH!D$Y#G%ofF_F(md1tYaWF<%)fhon%UIz!$)}i(yeQ6b zW6>Z6No2T$1Qe#XoSpZ5fh&L_?j{duHg`&}@ zL}&nmKV++_f+vW$`I|$|6TI_BSLrc>9=iGs`!Ew3XJ7LtVP^8$-*uep&HT0)yKzN# z6oPOa;SF10y01D1q-KLii|+W3L3HQcAgT)rkQ_i>7zO&cIBv1;@(UfGYQ1VVLOz{F zcP^mk0;RT|kb=3H+U7M2#cKipuvG;h}-Vv)65dyh3d_eK#2EwY7z&a+`dG~aWT z-iXZO4>~L2_DbUhCHiX>0T#h6;lm zu_1jt$H+?;vURwI|2D2XdqQ!z~`7`I9faXDzSFxrEGx5VR0j8o0`&aW>!{3 zNW)QLe!LcYld}{aVg71WB=ps?YbyTPv&}BloKyU&PML~fXCxnw<-B>OdTJFsUKw=} zoH}3OajqYp7w&{FajfXW`cb>C>iw8?KPT)p;+L(W(nZ;AD3n{8Hj8GrU1Dhrs)HSL z5_;|DU}4p4iR7~8Gt9YGaAMjn*Okm!- zj^kAKVi(<+r*F&K!+lp7D7QFDJ+(>~TS>zA5+3MH=7q{a?I4pPu$=ui^S{1J?n!(WHWp=1Lu1TvTIA7(`eb(aXGYe*l0xWgm~_4 z#}J9on6CoI-QZczcFE7%@Z;5}`L1Spw|$p|A+44up0kA~iEdoQurkh;FMxkOPBKAB z7!sh19$ozw!6M<5Bh6OaFK_k zsN5hSh6-l;pZgv~Sm$kj4!eku?uWLtfTgO*DT*AZh74Og9Yp>_ci&c3N_2%#P{D1lgK}OG5s{P;X*#oxcg+P3fzuR*> zxv^;auaAV#_?-3JyaZus@6Wluh_ZvAK(EMr$vlcXU3Uy@>5DisxO58}a_&X<6+mzd z3zTqP>zFx5dnVYW>yE+(Ehyt71}#7YidxVn1gK+0gjQ}($`hJ##nsF&T|=yTfx2P< zy;Cl~IG!<&%Plagt_BY+g@NA=Y=Z(HX1k^Id=k!^JkSju5Nno%)5e73UVF^ma(MIy z9-VQ19l#GfpzFF2W%E>2fWh|d* zOax)D?=-=SPw5f!ExapJO)s4{0HIb1m*mQx5S@(#>^E-1qR|D|+HR0U%gf1y zz=2yqy2b=rPk|_WPf!fW!25G+Jbs)e;*_?RUPM~HF@luPgGP;y?6kPs9lg;afu34; zEMhjUSGsU0LQwT6caovT2b|*L!!B04O?zF=! z_6hDQ!vGq!J3-UT8(b|Wo^zs}#^D$|33wXb@^WJzdj|>ODyn1jX>{_g$GrIFwqUO; zI+zZ8b)o&3l)Q>cBy8(3c^e=ki=9{fK2oXk{P}Yx)v@D?()zR|ts{zC%af;U$5z#2 z_{>Zp-ukn&ksWYHTce(wP$v4ab)rDd6gM3}&ShnHgiMWpQ6QcV$MqcztbDo?u<}jl z3h+^Q4fz;&Th)V+QC)TGoxLe}Dhy{L{6JiEs~Dz&TF1 zxLC_S>--7BDD}fxauG*QM%rb*J+`0Kd8#T3#Z_ezQU57XANiIu@GWC(1t($QfV&+s zBjdC?kUZOSQN=C`dPU9!HyUO16%->lnBf^Vu_wyb$&Um)=-64wnIIj?)k>*YEn3 zD1D`D3*)C5F_wAOE=S)52W%MLuFSEkj%=yhRH<-Vv(k>++P?47I%+*N07q;{fx&^7 zQjqW2Ayg=du#g^Ycm8=fLzR;Ky}ikjwi^TKKgeTFY(@VhwK zW}j}_KQG13omd0hr3rEutRpQ5ragk4;=Vbm&u-UmG z74&zoNr#?Wl*W+Y$)BEkfE5Otn&=S5zk`jX&(FhU9X)2z=Ok_&f`8;MtH|U@89(@6 Df2{ve diff --git a/man/staggered.Rd b/man/staggered.Rd index cd45f86..bda6504 100644 --- a/man/staggered.Rd +++ b/man/staggered.Rd @@ -18,7 +18,10 @@ staggered( use_DiD_A0 = ifelse(is.null(A_0_list), TRUE, FALSE), return_full_vcv = FALSE, return_matrix_list = FALSE, - use_last_treated_only = FALSE + use_last_treated_only = FALSE, + compute_fisher = FALSE, + num_fisher_permutations = 500, + skip_data_check = FALSE ) } \arguments{ @@ -49,6 +52,12 @@ staggered( \item{return_matrix_list}{If true, the function returns a list of the A_0_list and A_theta_list matrices along with betastar. This is used for internal recursive calls to calculate the variance-covariance matrix, and will generally not be needed by the end-user. Default is False.} \item{use_last_treated_only}{If true, then A_0_list and A_theta_list are created to only make comparisons with the last treated cohorts (as suggested by Sun and Abraham), rather than using not-yet-treated units as comparisons. If set to TRUE (and use_DiD_A0 = TRUE), then beta=1 corresponds with the Sun and Abraham estimator.} + +\item{compute_fisher}{If true, computes a Fisher Randomization Test using the studentized estimator.} + +\item{num_fisher_permutations}{The number of permutations to use in the Fisher Randomization Test (if compute_fisher = TRUE). Default is 500.} + +\item{skip_data_check}{If true, skips checks that the data is balanced and contains the colums i,t,g,y. Used in internal recursive calls to increase speed, but not recommended for end-user.} } \value{ resultsDF A data.frame containing: estimate (the point estimate), se (the standard error), and se_neyman (the Neyman standard error). If a vector-valued eventTime is provided, the data.frame contains multiple rows for each eventTime and an eventTime column. If return_full_vcv = TRUE and estimand = "eventstudy", the function returns a list containing resultsDF and the full variance covariance for the event-study estimates (vcv) as well as the Neyman version of the covariance matrix (vcv_neyman). (If return_matrix_list = TRUE, it likewise returns a list containing lists of matrices used in the vcv calculation.) diff --git a/man/staggered_cs.Rd b/man/staggered_cs.Rd index df0eb7f..3657cb7 100644 --- a/man/staggered_cs.Rd +++ b/man/staggered_cs.Rd @@ -15,7 +15,10 @@ staggered_cs( A_0_list = NULL, eventTime = 0, return_full_vcv = FALSE, - return_matrix_list = FALSE + return_matrix_list = FALSE, + compute_fisher = FALSE, + num_fisher_permutations = 500, + skip_data_check = FALSE ) } \arguments{ @@ -40,6 +43,12 @@ staggered_cs( \item{return_full_vcv}{If this is true and estimand = "eventstudy", then the function returns a list containing the full variance-covariance matrix for the event-plot estimates in addition to the usual dataframe with the estimates} \item{return_matrix_list}{If true, the function returns a list of the A_0_list and A_theta_list matrices along with betastar. This is used for internal recursive calls to calculate the variance-covariance matrix, and will generally not be needed by the end-user. Default is False.} + +\item{compute_fisher}{If true, computes a Fisher Randomization Test using the studentized estimator.} + +\item{num_fisher_permutations}{The number of permutations to use in the Fisher Randomization Test (if compute_fisher = TRUE). Default is 500.} + +\item{skip_data_check}{If true, skips checks that the data is balanced and contains the colums i,t,g,y. Used in internal recursive calls to increase speed, but not recommended for end-user.} } \value{ resultsDF A data.frame containing: estimate (the point estimate), se (the standard error), and se_neyman (the Neyman standard error). If a vector-valued eventTime is provided, the data.frame contains multiple rows for each eventTime and an eventTime column. If return_full_vcv = TRUE and estimand = "eventstudy", the function returns a list containing resultsDF and the full variance covariance for the event-study estimates (vcv) as well as the Neyman version of the covariance matrix (vcv_neyman). (If return_matrix_list = TRUE, it likewise returns a list containing lists of matrices used in the vcv calculation.) diff --git a/man/staggered_sa.Rd b/man/staggered_sa.Rd index 6fe0c93..3a0aeb2 100644 --- a/man/staggered_sa.Rd +++ b/man/staggered_sa.Rd @@ -15,7 +15,10 @@ staggered_sa( A_0_list = NULL, eventTime = 0, return_full_vcv = FALSE, - return_matrix_list = FALSE + return_matrix_list = FALSE, + compute_fisher = FALSE, + num_fisher_permutations = 500, + skip_data_check = FALSE ) } \arguments{ @@ -40,6 +43,12 @@ staggered_sa( \item{return_full_vcv}{If this is true and estimand = "eventstudy", then the function returns a list containing the full variance-covariance matrix for the event-plot estimates in addition to the usual dataframe with the estimates} \item{return_matrix_list}{If true, the function returns a list of the A_0_list and A_theta_list matrices along with betastar. This is used for internal recursive calls to calculate the variance-covariance matrix, and will generally not be needed by the end-user. Default is False.} + +\item{compute_fisher}{If true, computes a Fisher Randomization Test using the studentized estimator.} + +\item{num_fisher_permutations}{The number of permutations to use in the Fisher Randomization Test (if compute_fisher = TRUE). Default is 500.} + +\item{skip_data_check}{If true, skips checks that the data is balanced and contains the colums i,t,g,y. Used in internal recursive calls to increase speed, but not recommended for end-user.} } \value{ resultsDF A data.frame containing: estimate (the point estimate), se (the standard error), and se_neyman (the Neyman standard error). If a vector-valued eventTime is provided, the data.frame contains multiple rows for each eventTime and an eventTime column. If return_full_vcv = TRUE and estimand = "eventstudy", the function returns a list containing resultsDF and the full variance covariance for the event-study estimates (vcv) as well as the Neyman version of the covariance matrix (vcv_neyman). (If return_matrix_list = TRUE, it likewise returns a list containing lists of matrices used in the vcv calculation.) From e3e42a0c54868dcc569e560d215d4c69d628f5af Mon Sep 17 00:00:00 2001 From: jonathanroth_MSFt Date: Wed, 28 Apr 2021 10:42:27 -0700 Subject: [PATCH 04/13] Changed cov calc from Rfast to coop --- DESCRIPTION | 4 ++-- R/compute_efficient_estimator_and_se.R | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 170317f..5f21135 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -19,8 +19,8 @@ Imports: magrittr, MASS, stats, - Rfast, - tidyr + tidyr, + coop LinkingTo: Rcpp, RcppEigen diff --git a/R/compute_efficient_estimator_and_se.R b/R/compute_efficient_estimator_and_se.R index 2554a58..e115b84 100644 --- a/R/compute_efficient_estimator_and_se.R +++ b/R/compute_efficient_estimator_and_se.R @@ -42,7 +42,7 @@ compute_g_level_summaries <- function(df, is_balanced = TRUE){ #Compute means Ybar_g and covariance S_g Ybar_g <- base::colMeans(dfg) - S_g <- Rfast::cova(dfg) + S_g <- coop::covar(dfg) return(list(Ybar_g = Ybar_g, From f55dfd6a3ce6c5eee726be9be671779c3c9fb248 Mon Sep 17 00:00:00 2001 From: jonathanroth_MSFt Date: Wed, 28 Apr 2021 11:09:17 -0700 Subject: [PATCH 05/13] Added Rcpp to namespace and put dplyr:: in front of ungroup() --- NAMESPACE | 1 + R/compute_efficient_estimator_and_se.R | 3 ++- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/NAMESPACE b/NAMESPACE index 7a1d74f..566a9d7 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -3,5 +3,6 @@ export(staggered) export(staggered_cs) export(staggered_sa) +import(Rcpp) importFrom(magrittr,"%>%") useDynLib(staggered) diff --git a/R/compute_efficient_estimator_and_se.R b/R/compute_efficient_estimator_and_se.R index e115b84..442ce43 100644 --- a/R/compute_efficient_estimator_and_se.R +++ b/R/compute_efficient_estimator_and_se.R @@ -21,7 +21,7 @@ compute_g_level_summaries <- function(df, is_balanced = TRUE){ #Reshape so that Y_{it} is a column for all t df <- df %>% tidyr::pivot_wider(id_cols = c("i","g"), names_from = "t", values_from = "y") %>% - ungroup() %>% + dplyr::ungroup() %>% dplyr::select(-i) compute_Ybar_Sbar_g <- function(g){ @@ -884,6 +884,7 @@ processDF <- function(df, i, g, t, y){ #' @useDynLib staggered #' @importFrom magrittr "%>%" +#' @import Rcpp #' @title Calculate the efficient adjusted estimator in staggered rollout designs #' @description This functions calculates the efficient estimator for staggered rollout designs proposed by Roth and Sant'Anna. #' @param df A data frame containing panel data with the variables y (an outcome), i (an individual identifier), t (the period in which the outcome is observe), g (the period in which i is first treated, with Inf denoting never treated) From 236ae67a75f6608bf5e90157b867cfdecf159f3e Mon Sep 17 00:00:00 2001 From: jonathanroth_MSFt Date: Wed, 28 Apr 2021 13:56:01 -0700 Subject: [PATCH 06/13] Update betastar calculation to use an rcpp solver which makes it faster, especially for the case with full X --- R/RcppExports.R | 4 ++++ R/compute_efficient_estimator_and_se.R | 3 ++- src/RcppExports.cpp | 13 +++++++++++++ src/code.cpp | 10 ++++++++++ 4 files changed, 29 insertions(+), 1 deletion(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index 5b1ef97..f83ffe0 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -5,3 +5,7 @@ eigenMapMatMult <- function(A, B) { .Call('_staggered_eigenMapMatMult', PACKAGE = 'staggered', A, B) } +solve_least_squares_svd <- function(A, B) { + .Call('_staggered_solve_least_squares_svd', PACKAGE = 'staggered', A, B) +} + diff --git a/R/compute_efficient_estimator_and_se.R b/R/compute_efficient_estimator_and_se.R index 442ce43..9f83bb5 100644 --- a/R/compute_efficient_estimator_and_se.R +++ b/R/compute_efficient_estimator_and_se.R @@ -184,7 +184,8 @@ compute_Betastar <- function(Ybar_g_list, X_theta_cov <- base::Reduce(x = X_theta_cov_list, f = '+') #betastar <- solve(Xvar) %*% X_theta_cov - betastar <- MASS::ginv(Xvar) %*% X_theta_cov + #betastar <- MASS::ginv(Xvar) %*% X_theta_cov + betastar <- solve_least_squares_svd(Xvar,X_theta_cov) return(betastar) } diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index be7d1c2..d8018d6 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -18,9 +18,22 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// solve_least_squares_svd +Eigen::MatrixXd solve_least_squares_svd(Eigen::MatrixXd A, Eigen::MatrixXd B); +RcppExport SEXP _staggered_solve_least_squares_svd(SEXP ASEXP, SEXP BSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< Eigen::MatrixXd >::type A(ASEXP); + Rcpp::traits::input_parameter< Eigen::MatrixXd >::type B(BSEXP); + rcpp_result_gen = Rcpp::wrap(solve_least_squares_svd(A, B)); + return rcpp_result_gen; +END_RCPP +} static const R_CallMethodDef CallEntries[] = { {"_staggered_eigenMapMatMult", (DL_FUNC) &_staggered_eigenMapMatMult, 2}, + {"_staggered_solve_least_squares_svd", (DL_FUNC) &_staggered_solve_least_squares_svd, 2}, {NULL, NULL, 0} }; diff --git a/src/code.cpp b/src/code.cpp index f6c8768..88d684e 100644 --- a/src/code.cpp +++ b/src/code.cpp @@ -10,3 +10,13 @@ SEXP eigenMapMatMult(const Eigen::Map A, Eigen::Map Date: Thu, 29 Apr 2021 13:14:58 -0700 Subject: [PATCH 07/13] Updated to use a faster solver in rcpp --- R/RcppExports.R | 4 ++++ R/compute_efficient_estimator_and_se.R | 3 ++- src/RcppExports.cpp | 13 +++++++++++++ src/code.cpp | 11 +++++++++++ 4 files changed, 30 insertions(+), 1 deletion(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index f83ffe0..95a1b85 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -9,3 +9,7 @@ solve_least_squares_svd <- function(A, B) { .Call('_staggered_solve_least_squares_svd', PACKAGE = 'staggered', A, B) } +solve_least_squares_normal <- function(A, B) { + .Call('_staggered_solve_least_squares_normal', PACKAGE = 'staggered', A, B) +} + diff --git a/R/compute_efficient_estimator_and_se.R b/R/compute_efficient_estimator_and_se.R index 9f83bb5..06ac3db 100644 --- a/R/compute_efficient_estimator_and_se.R +++ b/R/compute_efficient_estimator_and_se.R @@ -185,7 +185,8 @@ compute_Betastar <- function(Ybar_g_list, #betastar <- solve(Xvar) %*% X_theta_cov #betastar <- MASS::ginv(Xvar) %*% X_theta_cov - betastar <- solve_least_squares_svd(Xvar,X_theta_cov) + #betastar <- solve_least_squares_svd(Xvar,X_theta_cov) + betastar <- solve_least_squares_normal(Xvar,X_theta_cov) #fast method of solving (Xvar)^-1 X_theta_cov return(betastar) } diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index d8018d6..d34fdfd 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -30,10 +30,23 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// solve_least_squares_normal +Eigen::MatrixXd solve_least_squares_normal(Eigen::MatrixXd A, Eigen::MatrixXd B); +RcppExport SEXP _staggered_solve_least_squares_normal(SEXP ASEXP, SEXP BSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< Eigen::MatrixXd >::type A(ASEXP); + Rcpp::traits::input_parameter< Eigen::MatrixXd >::type B(BSEXP); + rcpp_result_gen = Rcpp::wrap(solve_least_squares_normal(A, B)); + return rcpp_result_gen; +END_RCPP +} static const R_CallMethodDef CallEntries[] = { {"_staggered_eigenMapMatMult", (DL_FUNC) &_staggered_eigenMapMatMult, 2}, {"_staggered_solve_least_squares_svd", (DL_FUNC) &_staggered_solve_least_squares_svd, 2}, + {"_staggered_solve_least_squares_normal", (DL_FUNC) &_staggered_solve_least_squares_normal, 2}, {NULL, NULL, 0} }; diff --git a/src/code.cpp b/src/code.cpp index 88d684e..b0e3eb7 100644 --- a/src/code.cpp +++ b/src/code.cpp @@ -20,3 +20,14 @@ Eigen::MatrixXd solve_least_squares_svd( Eigen::MatrixXd x = A.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(B); return(x); } + + +// [[Rcpp::export]] +Eigen::MatrixXd solve_least_squares_normal( + Eigen::MatrixXd A, + Eigen::MatrixXd B +) { + Eigen::MatrixXd x = (A.transpose() * A).ldlt().solve(A.transpose() * B); + return(x); +} + From ff484f0a7cc4d366de3ff197d56b3efbe3975844 Mon Sep 17 00:00:00 2001 From: jonathanroth_MSFt Date: Thu, 20 May 2021 21:01:16 -0700 Subject: [PATCH 08/13] Removed the dplyr::select in the frt, since this appears to create occational issues on Azure --- R/compute_efficient_estimator_and_se.R | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/R/compute_efficient_estimator_and_se.R b/R/compute_efficient_estimator_and_se.R index 06ac3db..f2b3e0f 100644 --- a/R/compute_efficient_estimator_and_se.R +++ b/R/compute_efficient_estimator_and_se.R @@ -1178,7 +1178,8 @@ staggered <- function(df, i_g_table$g <- i_g_table$g[randIndex] #Merge the new treatment assignments back with the original - df <- dplyr::left_join(df %>% dplyr::select(-g), + df$g <- NULL + df <- dplyr::left_join(df, i_g_table, by = c("i")) @@ -1188,10 +1189,13 @@ staggered <- function(df, if(compute_fisher){ #Find unique pairs of (i,g). This will be used for computing the permutations - i_g_table <- df %>% - dplyr::filter(t == min(t)) %>% - dplyr::select(i,g) + # i_g_table <- df %>% + # dplyr::filter(t == min(t)) %>% + # dplyr::select(i,g) + i_g_table <- df %>% + dplyr::filter(t == min(t)) + i_g_table <- i_g_table[,c("i","g")] FRTResults <- purrr::map_dfr(.x = 1:num_fisher_permutations, From eb739ee054b899482deee90406f350718f650582 Mon Sep 17 00:00:00 2001 From: jonathanroth_MSFt Date: Fri, 21 May 2021 10:30:10 -0700 Subject: [PATCH 09/13] Fixed bug in frt that was comparing SA/CS to efficient estimator --- R/compute_efficient_estimator_and_se.R | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/R/compute_efficient_estimator_and_se.R b/R/compute_efficient_estimator_and_se.R index f2b3e0f..ad0d612 100644 --- a/R/compute_efficient_estimator_and_se.R +++ b/R/compute_efficient_estimator_and_se.R @@ -1115,6 +1115,9 @@ staggered <- function(df, ) # Xvar_list <- purrr::pmap(.l = list(A_0_list, S_g_list, N_g_list) , .f = function(A0,S,N){ return(1/N * A0 %*%S %*% base::t(A0) ) } ) + #Save the user-inputted beta (used in FRT call) + user_input_beta <- beta + if(is.null(beta)){ beta <- compute_Betastar(Ybar_g_list, A_theta_list, @@ -1200,7 +1203,8 @@ staggered <- function(df, FRTResults <- purrr::map_dfr(.x = 1:num_fisher_permutations, .f = ~ staggered::staggered(df = permuteTreatment(df, i_g_table, seed = .x), - estimand = "estimand", + estimand = NULL, + beta = user_input_beta, A_theta_list = A_theta_list, A_0_list = A_0_list, eventTime = eventTime, From c98d40e95e7fc031e103d5bd5c3254081d2bf5ee Mon Sep 17 00:00:00 2001 From: Jonathan Roth Date: Fri, 21 May 2021 14:04:26 -0700 Subject: [PATCH 10/13] Update DESCRIPTION Updated required R version --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 7e3ac08..2c53ae6 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -19,4 +19,4 @@ LinkingTo: Rcpp, RcppEigen Depends: - R (>= 3.50) + R (>= 3.5.0) From 28206e3da29f540ce588ea7db05361169d1fc76d Mon Sep 17 00:00:00 2001 From: jonathanroth_MSFt Date: Mon, 24 May 2021 11:04:30 -0700 Subject: [PATCH 11/13] Added error handling for frt draws --- R/compute_efficient_estimator_and_se.R | 37 +++++++++++++++++--------- 1 file changed, 24 insertions(+), 13 deletions(-) diff --git a/R/compute_efficient_estimator_and_se.R b/R/compute_efficient_estimator_and_se.R index ad0d612..150572a 100644 --- a/R/compute_efficient_estimator_and_se.R +++ b/R/compute_efficient_estimator_and_se.R @@ -1200,23 +1200,34 @@ staggered <- function(df, dplyr::filter(t == min(t)) i_g_table <- i_g_table[,c("i","g")] + #Now, we compute the FRT for each seed, permuting treatment for each one + #We catch any errors in the FRT simulations, and throw a warning if at least one has an error (using the remaining draws to calculate frt) FRTResults <- - purrr::map_dfr(.x = 1:num_fisher_permutations, - .f = ~ staggered::staggered(df = permuteTreatment(df, i_g_table, seed = .x), - estimand = NULL, - beta = user_input_beta, - A_theta_list = A_theta_list, - A_0_list = A_0_list, - eventTime = eventTime, - return_full_vcv = F, - return_matrix_list = F, - compute_fisher = F, - skip_data_check = T) %>% - mutate(seed = .x) - ) + purrr::map(.x = 1:num_fisher_permutations, + .f = purrr::possibly( + .f =~ staggered::staggered(df = permuteTreatment(df, i_g_table, seed = .x), + estimand = NULL, + beta = user_input_beta, + A_theta_list = A_theta_list, + A_0_list = A_0_list, + eventTime = eventTime, + return_full_vcv = F, + return_matrix_list = F, + compute_fisher = F, + skip_data_check = T) %>% mutate(seed = .x), + otherwise = NULL) + ) %>% + purrr::discard(base::is.null) %>% + purrr::reduce(.f = dplyr::bind_rows) + + successful_frt_draws <- NROW(FRTResults) + if(successful_frt_draws < num_fisher_permutations){ + warning("There was an error in at least one of the FRT simulations. Removing the problematic draws.") + } resultsDF$fisher_pval <- mean( abs(resultsDF$estimate/resultsDF$se) < abs(FRTResults$estimate/FRTResults$se) ) resultsDF$fisher_pval_se_neyman <- mean( abs(resultsDF$estimate/resultsDF$se_neyman) < abs(FRTResults$estimate/FRTResults$se_neyman) ) + resultsDF$num_fisher_permutations <- successful_frt_draws } if(!return_matrix_list){ From 8385af29e4bb7ea500edff233d3ed6d94d2fde69 Mon Sep 17 00:00:00 2001 From: jonathanroth_MSFt Date: Fri, 4 Jun 2021 12:04:33 -0700 Subject: [PATCH 12/13] Fixed bug when there is only when ATE(t,g) identified --- R/compute_efficient_estimator_and_se.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/compute_efficient_estimator_and_se.R b/R/compute_efficient_estimator_and_se.R index e9b3ca7..2ab41c9 100644 --- a/R/compute_efficient_estimator_and_se.R +++ b/R/compute_efficient_estimator_and_se.R @@ -83,7 +83,7 @@ balance_df <- function(df){ it_counts <- df %>% dplyr::group_by(i,t) %>% - dplyr::summarise(n = n()) + dplyr::summarise(n = dplyr::n()) if(max(it_counts$n) > 1 ){ stop("There are multiple observations with the same (i,t) values. The panel should have a unique outcome for each (i,t) value.") @@ -433,7 +433,7 @@ create_A0_list <- function(g_list, sum_of_lists <- function(.l){ #Given a list of lists all of the same length, this returns a list where the jth entry is the sum of the jth entry of all the lists - if(length(.l) == 1){return(.l)} + if(length(.l) == 1){return(.l[[1]])} results <- purrr::reduce(.x = .l, .f = function(l1, l2){ From 12eab654b981057d961f162e2b7d185bb5e7830a Mon Sep 17 00:00:00 2001 From: jonathanroth_MSFt Date: Thu, 24 Jun 2021 15:34:55 -0700 Subject: [PATCH 13/13] Updated staggered_stata documentation to use version 2.5.0 of rcall --- README.Rmd | 6 ++++-- README.md | 12 ++++++++++-- 2 files changed, 14 insertions(+), 4 deletions(-) diff --git a/README.Rmd b/README.Rmd index bed5704..47cc4cf 100644 --- a/README.Rmd +++ b/README.Rmd @@ -121,9 +121,11 @@ We also provide a Stata implementation (staggered_stata) via the RCall package, To install the staggered_stata package, the user first needs to install the github and RCall packages. This can be done with the following commands > net install github, from("https://haghish.github.io/github/") - > github install haghish/rcall, stable + > github install haghish/rcall, version("2.5.0") -Note that the user must have R installed before installing the rcall package. The latest version of R can be downloaded [here](https://cloud.r-project.org/). The staggered_stata package can then be installed with +**Important note:** the staggered_stata package was built under Rcall version 2.5.0, and the recent release of RCall version 3.0 has created some compatibility issues. We will try our best to fix these issues, but in the meantime it is best to install version 2.5.0, as in the command above. + +Note that the user must have R installed before installing the RCall package. The latest version of R can be downloaded [here](https://cloud.r-project.org/). The staggered_stata package can then be installed with > github install jonathandroth/staggered_stata diff --git a/README.md b/README.md index cb0ae3f..d968c29 100644 --- a/README.md +++ b/README.md @@ -46,6 +46,7 @@ for modifying and plotting the results. ``` r library(staggered) #load the staggered package library(dplyr) #load dplyr for data manipulation +#> Warning: package 'dplyr' was built under R version 4.0.5 #> #> Attaching package: 'dplyr' #> The following objects are masked from 'package:stats': @@ -55,6 +56,7 @@ library(dplyr) #load dplyr for data manipulation #> #> intersect, setdiff, setequal, union library(ggplot2) #load ggplot2 for plotting the results +#> Warning: package 'ggplot2' was built under R version 4.0.5 library(purrr) df <- staggered::pj_officer_level_balanced #load the officer data @@ -180,9 +182,15 @@ the github and RCall packages. This can be done with the following commands > net install github, from("https://haghish.github.io/github/") - > github install haghish/rcall, stable + > github install haghish/rcall, version("2.5.0") -Note that the user must have R installed before installing the rcall +**Important note:** the staggered\_stata package was built under Rcall +version 2.5.0, and the recent release of RCall version 3.0 has created +some compatibility issues. We will try our best to fix these issues, but +in the meantime it is best to install version 2.5.0, as in the command +above. + +Note that the user must have R installed before installing the RCall package. The latest version of R can be downloaded [here](https://cloud.r-project.org/). The staggered\_stata package can then be installed with