correct_gsynth.R

## ========================================================================== ##
# Project: GIRFT Elective Hubs Evaluation 

# Team: Improvement Analytics Unit (IAU) at the Health Foundation

# Script: correct_gsynth.R

# Corresponding author: Stefano Conti (e. stefano.conti@health.org.uk)

# Last modified: 02 August 2024

# Description: 
# Derive inferences on effect parameters and counterfactual outcomes from a GSynth model fitted 
# via the gsynth::gsynth() function to correct for biased results from original code implementation 
# in the presence of a staggered intervention.  N.B.: function can be applied to either balanced or 
# staggered panel data-sets.

# Inputs: 
# Data-sets        : 'data.dat'
# Fitted model     : 'gsynth_mdl'
# Tuning parameters: 'time.vec', 'alpha_ci', 'dest.dir', 'out.prf'

# Outputs:
# All outputs saved into global environment (not externally)
# inf.dat (a data-frame of inferences around ATT, counterfactual outcome and counterfactual-to-intervention 
#          outcome percent change at and across pre- and post-intervention time points and at and across 
#          intervention units)

# Dependencies: 
# None undefined externally.

# Notes:
# None in particular.
## ========================================================================== ##

correct_gsynth.fn <- function(gsynth_mdl, 
                              data.dat, 
                              time.vec=NULL, 
                              alpha_ci=.05, 
                              dest.dir=getwd(), 
                              out.prf="att.delta"
                              )
  {
  # Function correctly implements staggering of an intervention in a GSynth Interactive
  # Fixed Effects regression framework and tabulates, both within a list and optionally 
  # as a .csv dump, inferences on ATT and outcome percent change parameters 
  # (i) across key study times, as well as (ii) by individual intervention 
  # units.
  # 
  # 
  # Arguments:
  # 
  # gsynth_mdl: a GSynth model object (i.e. a model object generated by gsynth())
  # data.dat  : a data-frame for the data-set training the GSynth IFE model
  # time.vec  : a vector of positive integers denoting key post-intervention times
  # alpha_ci  : a numeric probability of Type 1 Error for the H_0) ATT = 0 null 
  #             hypothesis
  # dest.dir  : a character string denoting the destination folder path for the 
  #             tabular output
  # out.prf   : a character string denoting a prefix for tabular output
  # 
  # 
  # Output: 
  # 
  # A data-frame and ASCII (.csv) table of inferences around the ATT (i.e. 
  # point estimate, standard error, confidence interval, p-value), outcome levels 
  # (i.e. mean observed and counterfactual post-intervention values), percent outcome 
  # change (i.e. point estimate, confidence interval) parameters and model diagnostic 
  # statistics.
  
  
  # Function inputs for testing (staggered intervention, weighted outcome):
  # 
  # rm(list=ls())
  # project.bck <- ""  # Set AWS S3 project bucket
  # gsynth_mdl <- s3read_using(readRDS,
  #                            object=file.path("gsynth results/gsynth objects/pooled/new_hub_elect_rate.rds"),
  #                            bucket=project.bck
  #                            )
  # data.dat <- readRDS("~/iaelecthubs1/Data/new.no_hub.rds")
  # time.vec <- seq.int(12)
  # alpha_ci <- .05
  # dest.dir <- file.path("~/iaelecthubs1/Output/")
  # out.prf <- "att.delta"
  
  # Function inputs for testing (staggered intervention, unweighted outcome):
  # 
  # rm(list=ls())
  # data(gsynth)
  # turnout <- cbind(turnout, w=rpois(nrow(turnout), lambda=50))
  # gsynth_mdl <- gsynth(turnout ~ policy_edr + policy_mail_in + policy_motor,
  #                      data=turnout, index=c("abb", "year"), force="two-way",
  #                      CV=TRUE, r=c(0, 5), se=TRUE, inference="parametric", nboots=1e3,
  #                      parallel=TRUE, cores=4
  #                      )
  # data.dat <- turnout
  # time.vec <- seq.int(4)
  # alpha_ci <- .05
  # dest.dir <- file.path("~/iaelecthubs1/Output/")
  # out.prf <- "att.delta"
  
  # Function inputs for testing (non-staggered intervention, weighted outcome):
  # 
  # rm(list=ls())
  # project.bck <- ""  # Set AWS S3 project bucket
  # gsynth_mdl <- s3read_using(readRDS,
  #                            object=file.path("gsynth results/gsynth objects/pooled/old_hub_elect_rate.rds"),
  #                            bucket=project.bck
  #                            )
  # data.dat <- readRDS("~/iaelecthubs1/Data/old.no_hub.rds")
  # time.vec <- seq.int(12)
  # alpha_ci <- .05
  # dest.dir <- file.path("~/iaelecthubs1/Output/")
  # out.prf <- "att.delta"
  
  # Function inputs for testing (non-staggered intervention, unweighted outcome):
  # 
  # rm(list=ls())
  # data(gsynth)
  # gsynth_mdl <- gsynth(Y ~ D + X1 + X2,
  #                      data=simdata, index=c("id", "time"), force="two-way",
  #                      CV=TRUE, r=c(0, 5), se=TRUE, inference="parametric",
  #                      nboots=1e3, parallel=TRUE, cores=4
  #                      )
  # data.dat <- simdata
  # time.vec <- seq.int(4)
  # alpha_ci <- .05
  # dest.dir <- file.path("~/iaelecthubs1/Output/")
  # out.prf <- "att.delta"
  
  
  ##############
  ## Preamble ##
  ##############
  
  require("gsynth")  # Load required GSynth routines library
  
  stopifnot(class(gsynth_mdl) == "gsynth", 
            is.data.frame(data.dat) & all(unlist(gsynth_mdl[c("Y", "D", "index", "X")]) %in% names(data.dat)), 
            is.null(time.vec) | 
              (min(time.vec) >= 1 & max(time.vec) <= with(gsynth_mdl, expr=T - min(T0)) & 
                 sort(time.vec) == seq.int(min(time.vec), max(time.vec))
               ), 
            alpha_ci > 0 | alpha_ci < 1, 
            is.null(dest.dir) | file.exists(dest.dir), 
            is.character(out.prf)
            )  # Set break conditions on input incongruencies
  
  
  time.vec <- if(is.null(time.vec)) 
    seq.int(.5 * (with(gsynth_mdl, expr=T - min(T0)))) else   # Set vector of post-intervention times of interest
      sort(time.vec)  # Sort vector of post-intervention times of interest
  
  
  ##########################
  ## Derive staggered     ##
  ## effects and outcomes ##
  ##########################
  
  ## Pre-intervention 
  ## effects staggering 
  
  eff.pre_stg.vec.ls <- with(gsynth_mdl, 
                             expr=sapply(id.tr, 
                                         FUN=function(trs) 
                                           rev(eff[seq.int(T0[match(trs, table=id.tr)]), trs]
                                               ), 
                                         simplify=FALSE
                                         )
                             )  # Derive list by intervention unit of vectors by pre-intervention time of staggered effect estimates
  
  
  eff.pre_stg.arr <- sapply(eff.pre_stg.vec.ls, 
                            FUN="[", i=seq.int(with(gsynth_mdl, expr=max(T0)))
                            )  # Derive 2d-array by pre-intervention time, intervention unit of staggered effect estimates
  
  rownames(eff.pre_stg.arr) <- - with(gsynth_mdl, expr=seq.int(max(T0)))  # Rename staggered pre-intervention effect estimates 2d-array row margin
  
  
  ## Pre-intervention 
  ## outcomes staggering 
  
  out.pre_stg.arr.ls <- with(gsynth_mdl, 
                             expr=sapply(id.tr, 
                                         FUN=function(trs) 
                                           cbind(tr=rev(Y.tr[seq.int(T0[match(trs, table=id.tr)]), trs]), 
                                                 ct=rev(Y.ct[seq.int(T0[match(trs, table=id.tr)]), trs])
                                                 ), 
                                         simplify=FALSE
                                         )
                             )  # Derive list by intervention unit of 2d-arrays by pre-intervention time, unit status of outcomes
  
  
  out.pre_stg.arr <- sapply(out.pre_stg.arr.ls, 
                            FUN=function(arr) 
                              apply(arr, 
                                    MARGIN=2, 
                                    FUN="[", i=seq.int(with(gsynth_mdl, expr=max(T0)))
                                    ), 
                            simplify="array"
                            )  # Derive 3d-array by pre-intervention time, unit status, intervention unit of staggered outcomes
  
  rownames(out.pre_stg.arr) <- rownames(eff.pre_stg.arr)  # Rename staggered pre-intervention outcomes 3d-array row margin
  
  
  ## Post-intervention 
  ## effects staggering 
  
  eff.post_stg.vec.ls <- with(gsynth_mdl, 
                              expr=sapply(id.tr, 
                                          FUN=function(trs) 
                                            eff[- seq.int(T0[match(trs, table=id.tr)]), trs], 
                                          simplify=FALSE
                                          )
                              )  # Derive list by intervention unit of vectors by post-intervention time of staggered effect estimates
  
  eff.post_stg.arr <- sapply(eff.post_stg.vec.ls, 
                             FUN="[", i=seq.int(with(gsynth_mdl, expr=T - min(T0)))
                             )  # Derive 2d-array by post-intervention time, intervention unit of staggered effect estimates
  
  rownames(eff.post_stg.arr) <- seq.int(with(gsynth_mdl, expr=T - min(T0)))  # Rename staggered post-intervention effect estimates 2d-array row margin
  
  
  ## Post-intervention 
  ## outcomes staggering 
  
  out.post_stg.arr.ls <- with(gsynth_mdl, 
                              expr=sapply(id.tr, 
                                          FUN=function(trs) 
                                            cbind(tr=Y.tr[- seq.int(T0[match(trs, table=id.tr)]), trs], 
                                                  ct=Y.ct[- seq.int(T0[match(trs, table=id.tr)]), trs]
                                                  ), 
                                          simplify=FALSE
                                          )
                              )  # Derive list by intervention unit of 2d-arrays by post-intervention time, unit status of outcomes
  
  
  out.post_stg.arr <- sapply(out.post_stg.arr.ls, 
                             FUN=function(arr) 
                               apply(arr, 
                                     MARGIN=2, 
                                     FUN="[", i=seq.int(with(gsynth_mdl, expr=T - min(T0)))
                                     ), 
                             simplify="array"
                             )  # Derive 3d-array by post-intervention time, unit status, intervention unit of staggered outcomes
  
  rownames(out.post_stg.arr) <- rownames(eff.post_stg.arr)  # Rename staggered outcomes 3d-array row margin
  
  
  ## Organise staggered effects 
  ## and outcomes in lists
  
  eff_stg.arr.ls <- list(pre=eff.pre_stg.arr, post=eff.post_stg.arr)  # Set list by study period of staggered effect estimates 2d-arrays
  
  out_stg.arr.ls <- list(pre=out.pre_stg.arr, post=out.post_stg.arr)  # Set list by study period of staggered outcomes 3d-arrays
  
  
  ######################
  ## Derive staggered ##
  ## outcome weights  ##
  ######################
  
  if(! is.null(gsynth_mdl$W))
    {
    unit_index.idx <- with(gsynth_mdl, 
                           expr=which(c(any(id.tr %in% unique(data.dat[[index[1]]])), 
                                        any(id.tr %in% unique(data.dat[[index[2]]]))
                                        )
                                      )
                           )  # Derive index of 'unit' in "index" gsynth() input
         
         
         time_index.idx <- 3 - unit_index.idx  # Derive index of 'time' in "index" gsynth() input
         
         
         wgt.frm <- with(gsynth_mdl, 
                         expr=as.formula(paste(W, 
                                               paste(index[c(time_index.idx, unit_index.idx)], 
                                                     collapse=" + "
                                                     ), 
                                               sep=" ~ "
                                               )
                                         )
                         )  # Set formula for formatting outcome weights data-frame to 2d-array by unit, time
         
         
         wgt.arr <- xtabs(wgt.frm, 
                          data=data.dat, 
                          subset=with(gsynth_mdl, 
                                      expr=data.dat[[index[unit_index.idx]]] %in% id.tr
                                      )
                          )  # Derive 2d-array by time, intervention unit of outcome weights
         
         
         ## Pre-intervention 
         ## weights staggering
         
         wgt.pre_stg.vec.ls <- with(gsynth_mdl, 
                                    expr=sapply(id.tr, 
                                                FUN=function(trs) 
                                                  rev(wgt.arr[seq.int(T0[match(trs, table=id.tr)]), trs]), 
                                                simplify=FALSE
                                                )
                                    )  # Derive list by intervention unit of vectors by pre-intervention time of staggered outcome weights
         
         
         wgt.pre_stg.arr <- sapply(wgt.pre_stg.vec.ls, 
                                   FUN="[", i=seq.int(max(gsynth_mdl$T0))
                                   )  # Derive 2d-array by pre-intervention time, intervention unit of staggered outcome weights
         
         rownames(wgt.pre_stg.arr) <- rownames(eff.pre_stg.arr)  # Rename staggered pre-intervention outcome weights 2d-array row margin
         
         
         wgt.pre.rel_time_stg.arr <- sweep(wgt.pre_stg.arr, 
                                           MARGIN=1, 
                                           STATS=rowSums(wgt.pre_stg.arr, na.rm=TRUE), 
                                           FUN="/"
                                           )  # Derive 2d-array by pre-intervention time, intervention unit of staggered time-relative outcome weights
         
         
         wgt.pre.rel_unit_stg.arr <- sweep(wgt.pre_stg.arr, 
                                           MARGIN=2, 
                                           STATS=colSums(wgt.pre_stg.arr, na.rm=TRUE), 
                                           FUN="/"
                                           )  # Derive 2d-array by pre-intervention time, intervention unit of staggered unit-relative outcome weights
         
         
         ## Weighted pre-intervention 
         ## mean staggered outcomes
         
         out.pre.tr_stg_time.vec <- rowSums(wgt.pre.rel_time_stg.arr * out_stg.arr.ls$pre[, "tr", ], 
                                            na.rm=TRUE
                                            )  # Derive vector by pre-intervention time of time-weighted mean staggered intervention outcomes
         
         
         out.pre.tr_stg_unit.vec <- colSums(wgt.pre.rel_unit_stg.arr * out_stg.arr.ls$pre[, "tr", ], 
                                            na.rm=TRUE
                                            )  # Derive vector by intervention unit of unit-weighted mean staggered intervention outcomes
         
         
         ## Post-intervention 
         ## weights staggering
         
         wgt.post_stg.vec.ls <- with(gsynth_mdl, 
                                     expr=sapply(id.tr, 
                                                 FUN=function(trs) 
                                                   wgt.arr[- seq.int(T0[match(trs, table=id.tr)]), trs], 
                                                 simplify=FALSE
                                                 )
                                     )  # Derive list by intervention unit of vectors by post-intervention time of staggered outcome weights
         
         
         wgt.post_stg.arr <- sapply(wgt.post_stg.vec.ls, 
                                    FUN="[", i=with(gsynth_mdl, 
                                                    expr=seq.int(T - min(T0))
                                                    )
                                    )  # Derive 2d-array by post-intervention time, intervention unit of staggered outcome weights
         
         rownames(wgt.post_stg.arr) <- rownames(eff.post_stg.arr)  # Rename staggered post-intervention outcome weights 2d-array row margin
         
         
         wgt.post.rel_time_stg.arr <- sweep(wgt.post_stg.arr, 
                                            MARGIN=1, 
                                            STATS=rowSums(wgt.post_stg.arr, na.rm=TRUE), 
                                            FUN="/"
                                            )  # Derive 2d-array by post-intervention time, intervention unit of staggered time-relative outcome weights
         
         
         wgt.post.rel_unit_stg.arr <- sweep(wgt.post_stg.arr, 
                                            MARGIN=2, 
                                            STATS=colSums(wgt.post_stg.arr, na.rm=TRUE), 
                                            FUN="/"
                                            )  # Derive 2d-array by post-intervention time, intervention unit of staggered unit-relative outcome weights
         
         
         ## Weighted post-intervention 
         ## mean staggered outcomes
         
         out.post.tr_stg_time.vec <- rowSums(wgt.post.rel_time_stg.arr * out_stg.arr.ls$post[, "tr", ], 
                                             na.rm=TRUE
                                             )  # Derive vector by post-intervention time of time-weighted mean staggered intervention outcomes
         
         
         out.post.tr_stg_unit.vec <- colSums(wgt.post.rel_unit_stg.arr * out_stg.arr.ls$post[, "tr", ], 
                                             na.rm=TRUE
                                             )  # Derive vector by intervention unit of unit-weighted mean staggered intervention outcomes
         } else 
           {
             ## Unweighted pre-intervention 
             ## mean staggered outcomes
             
             out.pre.tr_stg_time.vec <- rowMeans(out.pre_stg.arr[, "tr", ], na.rm=TRUE)  # Derive vector by pre-intervention time of mean staggered intervention outcomes
             
             
             out.pre.tr_stg_unit.vec <- colMeans(out.pre_stg.arr[, "tr", ], na.rm=TRUE)  # Derive vector by pre-intervention time of mean staggered intervention outcomes
             
             
             ## Unweighted post-intervention 
             ## mean staggered outcomes
             
             out.post.tr_stg_time.vec <- rowMeans(out.post_stg.arr[, "tr", ], na.rm=TRUE)  # Derive vector by post-intervention time of mean staggered intervention outcomes
             
             
             out.post.tr_stg_unit.vec <- colMeans(out.post_stg.arr[, "tr", ], na.rm=TRUE)  # Derive vector by post-intervention time of mean staggered intervention outcomes
             }
  
  
  ## Organise un- / weighted 
  ## staggered ATT and outcomes 
  ## estimates in lists
  
  att_stg_time.vec.ls <- sapply(eff_stg.arr.ls, 
                                FUN=rowMeans, na.rm=TRUE, simplify=FALSE
                                )  # Derive list by study period of vectors by study time of ATT estimates
  
  
  att_stg_unit.vec.ls <- sapply(eff_stg.arr.ls, 
                                FUN=colMeans, na.rm=TRUE, simplify=FALSE
                                )  # Derive list by study period of vectors by intervention unit of ATT estimates
  
  
  att_stg.vec.ls.ls <- list(pre=list(time=att_stg_time.vec.ls$pre, 
                                     unit=att_stg_unit.vec.ls$pre
                                     ), 
                            post=list(time=att_stg_time.vec.ls$post, 
                                      unit=att_stg_unit.vec.ls$post
                                      )
                            )  # Set list by study period of list by effect margin of ATT estimates vectors
  
  
  att_stg_mean.arr <- sapply(att_stg.vec.ls.ls, 
                             FUN=function(lst) 
                               sapply(lst, FUN=mean)
                             )  # Derive 2d-array by effect margin, study period of average ATT estimates
  
  
  out.tr_stg.vec.ls.ls <- list(pre=list(time=out.pre.tr_stg_time.vec, 
                                        unit=out.pre.tr_stg_unit.vec
                                        ), 
                               post=list(time=out.post.tr_stg_time.vec, 
                                         unit=out.post.tr_stg_unit.vec
                                         )
                               )  # List by study period of lists by effect margin of mean staggered outcome vectors
  
  
  out_stg_mean.arr.ls.ls <- list(pre=list(time=cbind(tr=out.tr_stg.vec.ls.ls$pre$time, 
                                                     ct=out.tr_stg.vec.ls.ls$pre$time - att_stg.vec.ls.ls$pre$time
                                                     ),   # Derive 2d-array by pre-intervention time, comparator of time-weighted mean staggered outcomes
                                          unit=cbind(tr=out.tr_stg.vec.ls.ls$pre$unit, 
                                                     ct=out.tr_stg.vec.ls.ls$pre$unit - att_stg.vec.ls.ls$pre$unit
                                                     )  # Derive 2d-array by pre-intervention time, comparator of unit-weighted mean staggered outcomes
                                          ),   # Derive list by effect margin of 2d-arrays by pre-intervention time, comparator of weighted mean staggered outcomes
                                 post=list(time=cbind(tr=out.tr_stg.vec.ls.ls$post$time, 
                                                      ct=out.tr_stg.vec.ls.ls$post$time - att_stg.vec.ls.ls$post$time
                                                      ),   # Derive 2d-array by post-intervention time, comparator of time-weighted mean staggered outcomes
                                           unit=cbind(tr=out.tr_stg.vec.ls.ls$post$unit, 
                                                      ct=out.tr_stg.vec.ls.ls$post$unit - att_stg.vec.ls.ls$post$unit
                                                      )  # Derive 2d-array by post-intervention time, comparator of unit-weighted mean staggered outcomes
                                           )  # Derive list by effect margin of 2d-arrays by post-intervention time, comparator of weighted mean staggered outcomes
                                 )  # Derive list by study period of lists by effect margin of 2d-arrays by study time, comparator of weighted mean staggered outcomes
  
  
  ###########################
  ## Derive staggered      ##
  ## effects and outcomes  ##
  ## from bootstrap sample ##
  ###########################
  
  eff_stg_boot.arr.ls.ls <- with(gsynth_mdl, 
                                 expr=list(pre=sapply(id.tr, 
                                                      FUN=function(trs) 
                                                        apply(eff.boot[seq.int(T0[match(trs, table=id.tr)]), trs, , drop=FALSE], 
                                                              MARGIN=2:3, 
                                                              FUN=rev
                                                              ), 
                                                      simplify=FALSE
                                                      ),   # Derive list by intervention unit of degenerate 3d-arrays by pre-intervention time, bootstrap draw of staggered effect estimates
                                           post=sapply(id.tr, 
                                                       FUN=function(trs) 
                                                         eff.boot[- seq.int(T0[match(trs, table=id.tr)]), trs, , drop=FALSE], 
                                                       simplify=FALSE
                                                       )  # Derive list by intervention unit of degenerate 3d-arrays by post-intervention time, bootstrap draw of staggered effect estimates
                                           ), 
                                 simplify=FALSE
                                 )  # Derive list by study time of lists by intervention unit of degenerate 3d-arrays by study time, bootstrap draw of staggered effect estimates
  
  
  eff_stg_boot.arr.ls <- list(pre=sapply(eff_stg_boot.arr.ls.ls$pre, 
                                         FUN=function(arr) 
                                           apply(arr, 
                                                 MARGIN=3, 
                                                 FUN="[", i=seq.int(with(gsynth_mdl, expr=max(T0)))
                                                 ), 
                                         simplify="array"
                                         ),   # Derive 3d-array by pre-intervention time, bootstrap draw, intervention unit of staggered effect estimates
                              post=sapply(eff_stg_boot.arr.ls.ls$post, 
                                          FUN=function(arr) 
                                            apply(arr, 
                                                  MARGIN=3, 
                                                  FUN="[", i=seq.int(with(gsynth_mdl, expr=T - min(T0)))
                                                  ), 
                                          simplify="array"
                                          )  # Derive 3d-array by post-intervention time, bootstrap draw, intervention unit of staggered effect estimates
                              )  # Derive list by study period of 3d-arrays by study time, bootstrap draw, intervention unit of staggered effect estimates
  
  rownames(eff_stg_boot.arr.ls$pre) <- rownames(eff_stg.arr.ls$pre)  # Rename staggered pre-intervention effect estimates 3d-array row margin
  
  rownames(eff_stg_boot.arr.ls$post) <- rownames(eff_stg.arr.ls$post)  # Rename staggered post-intervention effect estimates 3d-array row margin
  
  
  eff_stg_mean_boot.arr.ls.ls <- list(pre=list(time=apply(eff_stg_boot.arr.ls$pre, 
                                                          MARGIN=seq.int(2), 
                                                          FUN=mean, na.rm=TRUE
                                                          ),   # Derive 2d-array by pre-intervention time, bootstrap draw of mean time-staggered pre-intervention effect estimates
                                               unit=apply(eff_stg_boot.arr.ls$pre, 
                                                          MARGIN=3:2, 
                                                          FUN=mean, na.rm=TRUE
                                                          )  # Derive 2d-array by intervention trust, bootstrap draw of mean unit-staggered pre-intervention effect estimates
                                               ),   # Derive list by effect margin of 2d-arrays by intervention trust, bootstrap draw of mean staggered pre-intervention effect estimates
                                      post=list(time=apply(eff_stg_boot.arr.ls$post, 
                                                           MARGIN=seq.int(2), 
                                                           FUN=mean, na.rm=TRUE
                                                           ),   # Derive 2d-array by post-intervention time, bootstrap draw of mean time-staggered post-intervention effect estimates
                                                unit=apply(eff_stg_boot.arr.ls$post, 
                                                           MARGIN=3:2, 
                                                           FUN=mean, na.rm=TRUE
                                                           )  # Derive 2d-array by intervention trust, bootstrap draw of mean unit-staggered post-intervention effect estimates
                                                )  # Derive list by effect margin of 2d-arrays by intervention trust, bootstrap draw of mean staggered post-intervention effect estimates
                                      )  # Derive list by study period of lists by effect margin of 2d-arrays by intervention trust, bootstrap draw of mean staggered effect estimates
  
  
  out.ct_stg_boot.arr.ls <- sapply(names(eff_stg_boot.arr.ls), 
                                   FUN=function(whn) 
                                     sweep(-eff_stg_boot.arr.ls[[whn]], 
                                           MARGIN=c(1, 3), 
                                           STATS=out_stg.arr.ls[[whn]][, "tr", ], 
                                           FUN="+"
                                           ), 
                                   simplify=FALSE
                                   )  # Derive list by study period of 3d-arrays by study time, bootstrap draw, intervention unit of staggered counterfactual outcomes
  
  
  out.ct_stg_mean_boot.arr.ls.ls <- list(pre=list(time=apply(out.ct_stg_boot.arr.ls$pre, 
                                                             MARGIN=seq.int(2), 
                                                             FUN=mean, na.rm=TRUE
                                                             ),   # Derive 2d-array by pre-intervention time, bootstrap draw of mean pre-intervention time-staggered counterfactual outcomes
                                                  unit=apply(out.ct_stg_boot.arr.ls$pre, 
                                                             MARGIN=3:2, 
                                                             FUN=mean, na.rm=TRUE
                                                             )  # Derive 2d-array by intervention unit, bootstrap draw of mean pre-intervention unit-staggered counterfactual outcomes
                                                  ),   # Derive list by effect margin of 2d-arrays by margin, bootstrap draw of mean pre-intervention staggered counterfactual outcomes
                                         post=list(time=apply(out.ct_stg_boot.arr.ls$post, 
                                                             MARGIN=seq.int(2), 
                                                             FUN=mean, na.rm=TRUE
                                                             ),   # Derive 2d-array by post-intervention time, bootstrap draw of mean post-intervention time-staggered counterfactual outcomes
                                                   unit=apply(out.ct_stg_boot.arr.ls$post, 
                                                              MARGIN=3:2, 
                                                              FUN=mean, na.rm=TRUE
                                                              )  # Derive 2d-array by intervention unit, bootstrap draw of mean intervention unit-staggered counterfactual outcomes
                                                   )  # Derive list by effect margin of 2d-arrays by margin, bootstrap draw of mean post-intervention staggered counterfactual outcomes
                                         )  # Derive list by study period of lists by effect margin of 2d-arrays by margin, bootstrap draw of mean staggered counterfactual outcomes
  
  
  ############################
  ## Derive inferences on   ##
  ## ATT and outcome change ##
  ## across units and by    ##
  ## study period           ##
  ############################
  
  att.out_mean.dat <- with(gsynth_mdl, 
                           expr=data.frame(comparator=paste(ifelse(min(T0) == max(T0), "old", "new"), "hub", sep="_"), 
                                           outcome=Y,   # Set outcome name
                                           time=c(paste(-1, -max(T0), sep=" - "), 
                                                  paste(1, T - min(T0), sep=" - ")
                                                  ),   # Set intervention times of interest
                                           unit="Pooled",   # Set intervention units of interest
                                           att_mean=c(att_stg_mean.arr["time", "pre"], 
                                                      att_stg_mean.arr["time", "post"]
                                                      ),   # derive vector by study period of average ATT estimate across intervention units
                                           att_se=c(sd(colMeans(eff_stg_mean_boot.arr.ls.ls$pre$time)), 
                                                    sd(colMeans(eff_stg_mean_boot.arr.ls.ls$post$time))
                                                    )  # Derive vector by study period of ATT standard error across intervention units
                                           )
                           )  # Derive data-frame of ATT inferences across intervention units, post-intervention times
  
  att.out_mean.dat <- within(att.out_mean.dat, 
                             expr=
                               {
                                 att_ci.lower=att_mean + qnorm(.5 * alpha_ci) * att_se  # Derive vector by study period of ATT lower confidence bound across intervention units
                                 
                                 att_ci.upper=att_mean + qnorm(1 - .5 * alpha_ci) * att_se  # Derive vector by study period of ATT upper confidence bound across intervention units
                                 
                                 pval=2 * pnorm(abs(att_mean) / att_se, 
                                                lower.tail=FALSE
                                                )  # Derive vector by study period of ATT p-values across intervention units
                                 
                                 out_tr=c(mean(out_stg_mean.arr.ls.ls$pre$time[, "tr"]), 
                                          mean(out_stg_mean.arr.ls.ls$post$time[, "tr"])
                                          )  # Derive vector by study period of mean staggered intervention outcomes across intervention units
                                 
                                 out_ct=c(mean(out_stg_mean.arr.ls.ls$pre$time[, "ct"]), 
                                          mean(out_stg_mean.arr.ls.ls$post$time[, "ct"])
                                          )  # Derive vector by study period of mean staggered counterfactual outcome across intervention units
                                 
                                 out_ct_se=c(sd(colMeans(out.ct_stg_mean_boot.arr.ls.ls$pre$time)), 
                                             sd(colMeans(out.ct_stg_mean_boot.arr.ls.ls$post$time))
                                             )  # Derive vector by study period of staggered counterfactual outcome standard errors across intervention units
                                 
                                 out_ct_ci.lower=out_ct + qnorm(.5 * alpha_ci) * out_ct_se  # Derive vector by study period of staggered counterfactual outcome lower confidence bound across intervention units
                                 
                                 out_ct_ci.upper=out_ct + qnorm(1 - .5 * alpha_ci) * out_ct_se  # Derive vector by study period of staggered counterfactual outcome upper confidence bound across intervention units
                                 
                                 delta_mean=1e2 * att_mean / out_ct  # Derive vector by study period of mean staggered outcome change across intervention units
                                 
                                 delta_ci.lower=1e2 * att_ci.lower / out_ct  # Derive vector by study period of staggered outcome change approximate lower confidence bound across intervention units
                                 
                                 delta_ci.upper=1e2 * att_ci.upper / out_ct  # Derive vector by study period of staggered outcome change approximate upper confidence bound across intervention units
                                 }
                             )  # Augment average ATT inferences data-frame with second-order and outcome change inferences across intervention units, intervention times
  
  att.out_mean.dat <- with(gsynth_mdl, 
                           expr=within(att.out_mean.dat, 
                                       expr=
                                         {
                                           sigma2=sigma2  # Derive residuals variance
                                           
                                           size_fct=unname(r.cv)  # Derive cross-validated latent factors number
                                           
                                           size_period=c(max(T0), T - min(T0))  # Derive vector by study period of study period lengths across intervention units
                                           
                                           size_tr=c(0, Ntr)  # Derive average intervention units sample size across intervention units, post-intervention times
                                           
                                           size_ct=N - size_tr  # Derive average control units sample size across intervention units, post-intervention times
                                           }
                                       )
                           )  # Augment average ATT inferences data-frame with model diagnostics across intervention units, post-intervention times
  
  
  ############################
  ## Derive inferences on   ##
  ## ATT and outcome change ##
  ## across units and by    ##
  ## study time             ##
  ############################
  
  att.out_time.dat <- with(gsynth_mdl, 
                           expr=data.frame(comparator=paste(ifelse(min(T0) == max(T0), "old", "new"), "hub", sep="_"), 
                                           outcome=Y,   # Set outcome name
                                           time=c(- seq.int(max(T0)), seq.int(T - min(T0))),   # Set intervention times of interest
                                           unit="Pooled",   # Set intervention units of interest
                                           att_mean=c(att_stg.vec.ls.ls$pre$time, 
                                                      att_stg.vec.ls.ls$post$time
                                                      ),   # Derive vector by study times of average ATT estimates across intervention units
                                           att_se=c(apply(eff_stg_mean_boot.arr.ls.ls$pre$time, 
                                                          MARGIN=1, 
                                                          FUN=sd
                                                          ), 
                                                    apply(eff_stg_mean_boot.arr.ls.ls$post$time, 
                                                          MARGIN=1, 
                                                          FUN=sd
                                                          )
                                                    )  # Derive vector by study times of ATT standard errors across intervention units
                                           )
                           )  # Derive data-frame of ATT inferences by study time across intervention units
  
  att.out_time.dat <- within(att.out_time.dat, 
                             expr=
                               {
                                 att_ci.lower=att_mean + qnorm(.5 * alpha_ci) * att_se  # Derive vector by study time of ATT lower confidence bounds across intervention units
                                 
                                 att_ci.upper=att_mean + qnorm(1 - .5 * alpha_ci) * att_se  # Derive vector by study time of ATT upper confidence bounds across intervention units
                                 
                                 pval=2 * pnorm(abs(att_mean) / att_se, 
                                                lower.tail=FALSE
                                                )  # Derive vector by study time of ATT p-values across intervention units
                                 
                                 out_tr=c(out_stg_mean.arr.ls.ls$pre$time[, "tr"], 
                                          out_stg_mean.arr.ls.ls$post$time[, "tr"]
                                          )  # Derive vector by study time of mean staggered intervention outcomes across intervention units
                                 
                                 out_ct=c(out_stg_mean.arr.ls.ls$pre$time[, "ct"], 
                                          out_stg_mean.arr.ls.ls$post$time[, "ct"]
                                          )  # Derive vector by study time of mean staggered counterfactual outcomes across intervention units
                                 
                                 out_ct_se=apply(rbind(out.ct_stg_mean_boot.arr.ls.ls$pre$time, 
                                                       out.ct_stg_mean_boot.arr.ls.ls$post$time
                                                       ), 
                                                 MARGIN=1, 
                                                 FUN=sd
                                                 )  # Derive vector by study time of staggered counterfactual outcomes standard errors across intervention units
                                 
                                 out_ct_ci.lower=out_ct + qnorm(.5 * alpha_ci) * out_ct_se  # Derive vector by post-intervention time of staggered counterfactual outcome lower confidence bound across intervention units
                                 
                                 out_ct_ci.upper=out_ct + qnorm(1 - .5 * alpha_ci) * out_ct_se  # Derive vector by post-intervention time of staggered counterfactual outcome upper confidence bound across intervention units
                                 
                                 delta_mean=1e2 * att_mean / out_ct  # Derive vector by post-intervention time of mean staggered outcome change across intervention units
                                 
                                 delta_ci.lower=1e2 * att_ci.lower / out_ct  # Derive vector by post-intervention time of staggered outcome change approximate lower confidence bound across intervention units
                                 
                                 delta_ci.upper=1e2 * att_ci.upper / out_ct  # Derive vector by post-intervention time of staggered outcome change approximate upper confidence bound across intervention units
                                 }
                             )  # Augment ATT inferences data-frame with second-order and outcome change inferences by post-intervention time across intervention units
  
  att.out_time.dat <- with(gsynth_mdl, 
                           expr=within(att.out_time.dat, 
                                       expr=
                                         {
                                           sigma2=sigma2  # Derive residuals variance
                                           
                                           size_fct=unname(r.cv)  # Derive cross-validated latent factors number
                                           
                                           size_period=c(rep.int(max(T0), times=max(T0)), 
                                                  rep.int(T - min(T0), times=T - min(T0))
                                                  )  # Derive vector by study time of study period lengths across intervention units
                                           
                                           size_tr=c(rep.int(0, times=max(T0)), 
                                                     rev(rowSums(D.tr)[-seq.int(min(T0))])
                                                     )  # Derive vector by study time of intervention units sample size
                                           
                                           size_ct=N - size_tr  # Derive vector by study time of control units sample size across intervention units
                                           }
                                       )
                           )  # Augment ATT inferences data-frame with model diagnostics by study time across intervention units
  
  
  #############################
  ## Derive inferences on    ##
  ## ATT and outcome change  ##
  ## across units and across ##
  ## key study times         ##
  #############################
  
  att.out_time_mean.dat <- with(gsynth_mdl, 
                                expr=data.frame(comparator=paste(ifelse(min(T0) == max(T0), "old", "new"), "hub", sep="_"), 
                                                outcome=Y,   # Set outcome name
                                                time=c(paste(-time.vec[c(1, length(time.vec))], collapse=" - "), 
                                                       paste(time.vec[c(1, length(time.vec))], collapse=" - ")
                                                       ),   # Set intervention times of interest
                                                unit="Pooled",   # Set intervention units of interest
                                                att_mean=c(mean(att_stg.vec.ls.ls$pre$time[as.character(-time.vec)]), 
                                                           mean(att_stg.vec.ls.ls$post$time[as.character(time.vec)])
                                                           ),   # Derive vector by study periods of average ATT estimates across intervention units, key study times
                                                att_se=c(sd(colMeans(eff_stg_mean_boot.arr.ls.ls$pre$time[as.character(-time.vec), ])), 
                                                         sd(colMeans(eff_stg_mean_boot.arr.ls.ls$post$time[as.character(time.vec), ]))
                                                         )  # Derive vector by study period of ATT standard errors across intervention units, key study times
                                                )
                                )  # Derive data-frame of mean ATT inferences across intervention units, key study times
  
  att.out_time_mean.dat <- within(att.out_time_mean.dat, 
                                  expr=
                                    {
                                      att_ci.lower=att_mean + qnorm(.5 * alpha_ci) * att_se  # Derive vector by study period of ATT lower confidence bounds across intervention units, key study times
                                      
                                      att_ci.upper=att_mean + qnorm(1 - .5 * alpha_ci) * att_se  # Derive vector by study period of ATT upper confidence bounds across intervention units, key study times
                                      
                                      pval=2 * pnorm(abs(att_mean) / att_se, 
                                                     lower.tail=FALSE
                                                     )  # Derive vector by study period of ATT p-values across intervention units, key study times
                                      
                                      out_tr=c(mean(out_stg_mean.arr.ls.ls$pre$time[as.character(-time.vec), "tr"]), 
                                               mean(out_stg_mean.arr.ls.ls$post$time[as.character(time.vec), "tr"])
                                               )  # Derive vector by study period of mean staggered intervention outcome across intervention units, key study times
                                      
                                      out_ct=c(mean(out_stg_mean.arr.ls.ls$pre$time[as.character(-time.vec), "ct"]), 
                                               mean(out_stg_mean.arr.ls.ls$post$time[as.character(time.vec), "ct"])
                                               )  # Derive vector by study period of mean staggered counterfactual outcomes across intervention units, key study times
                                      
                                      out_ct_se=c(sd(colMeans(out.ct_stg_mean_boot.arr.ls.ls$pre$time[as.character(-time.vec), ])), 
                                                  sd(colMeans(out.ct_stg_mean_boot.arr.ls.ls$post$time[as.character(time.vec), ]))
                                                  )  # Derive vector by study period of staggered counterfactual outcomes standard errors across intervention units, key study times
                                      
                                      out_ct_ci.lower=out_ct + qnorm(.5 * alpha_ci) * out_ct_se  # Derive vector by study period of staggered counterfactual outcome lower confidence bounds across intervention units, key study times
                                      
                                      out_ct_ci.upper=out_ct + qnorm(1 - .5 * alpha_ci) * out_ct_se  # Derive vector by study period of staggered counterfactual outcome upper confidence bounds across intervention units, key study times
                                      
                                      delta_mean=1e2 * att_mean / out_ct  # Derive vector by study period of mean staggered outcome change across intervention units, key study times
                                      
                                      delta_ci.lower=1e2 * att_ci.lower / out_ct  # Derive vector by study period of staggered outcome change approximate lower confidence bounds across intervention units, key study times
                                      
                                      delta_ci.upper=1e2 * att_ci.upper / out_ct  # Derive vector by study period of staggered outcome change approximate upper confidence bounds across intervention units, key study times
                                      }
                                  )  # Augment average ATT inferences data-frame with second-order and outcome change inferences across intervention units, key study times
  
  att.out_time_mean.dat <- with(gsynth_mdl, 
                                expr=within(att.out_time_mean.dat, 
                                            expr=
                                              {
                                                sigma2=sigma2  # Derive residuals variance
                                                
                                                size_fct=unname(r.cv)  # Derive cross-validated latent factors number
                                                
                                                size_period=pmin(c(max(T0), T - min(T0)), length(time.vec))  # Derive study period length across intervention units
                                                
                                                size_tr=c(0, rev(rowSums(D.tr[-seq.int(min(T0)), ]))[max(time.vec)])  # Derive vector by study period of intervention units sample size
                                                
                                                size_ct=N - size_tr  # Derive vector by study period of control units sample size across intervention units
                                                }
                                            )
                                )  # Augment average ATT inferences data-frame with model diagnostics across intervention units, post-intervention times
  
  
  ############################
  ## Derive inferences on   ##
  ## ATT and outcome change ##
  ## by unit and across     ##
  ## study times            ##
  ############################
  
  att.out_unit.dat <- with(gsynth_mdl, 
                           expr=data.frame(comparator=paste(ifelse(min(T0) == max(T0), "old", "new"), "hub", sep="_"), 
                                           outcome=Y,   # Set outcome name
                                           time=rep(c(paste(-1, -max(T0), sep=" - "), 
                                                      paste(1, T - min(T0), sep=" - ")
                                                      ), 
                                                    each=Ntr
                                                    ),   # Set intervention times of interest
                                           unit=id.tr,   # Set intervention units of interest
                                           att_mean=c(att_stg.vec.ls.ls$pre$unit, 
                                                      att_stg.vec.ls.ls$post$unit
                                                      ),   # Derive vector by intervention unit of average ATT estimates across study periods
                                           att_se=c(apply(eff_stg_mean_boot.arr.ls.ls$pre$unit, 
                                                          MARGIN=1, 
                                                          FUN=sd
                                                          ), 
                                                    apply(eff_stg_mean_boot.arr.ls.ls$post$unit, 
                                                          MARGIN=1, 
                                                          FUN=sd
                                                          )
                                                    )  # Derive vector by intervention unit of ATT standard errors across study periods
                                           )
                           )  # Derive data-frame of ATT inferences by intervention unit across study periods
  
  att.out_unit.dat <- within(att.out_unit.dat, 
                             expr=
                               {
                                 att_ci.lower=att_mean + qnorm(.5 * alpha_ci) * att_se  # Derive vector by intervention unit of ATT lower confidence bounds across study periods
                                 
                                 att_ci.upper=att_mean + qnorm(1 - .5 * alpha_ci) * att_se  # Derive vector by intervention unit of ATT upper confidence bounds across study periods
                                 
                                 pval=2 * pnorm(abs(att_mean) / att_se, 
                                                lower.tail=FALSE
                                                )  # Derive vector by intervention unit of ATT p-values across study periods
                                 
                                 out_tr=c(out_stg_mean.arr.ls.ls$pre$unit[, "tr"], 
                                          out_stg_mean.arr.ls.ls$post$unit[, "tr"]
                                          )  # Derive vector by intervention unit of mean staggered intervention outcomes across study periods
                                 
                                 out_ct=c(out_stg_mean.arr.ls.ls$pre$unit[, "ct"], 
                                          out_stg_mean.arr.ls.ls$post$unit[, "ct"]
                                          )  # Derive vector by intervention unit of mean staggered counterfactual outcomes across study periods
                                 
                                 out_ct_se=c(apply(out.ct_stg_mean_boot.arr.ls.ls$pre$unit, 
                                                   MARGIN=1, 
                                                   FUN=sd
                                                   ), 
                                             apply(out.ct_stg_mean_boot.arr.ls.ls$post$unit, 
                                                   MARGIN=1, 
                                                   FUN=sd
                                                   )
                                             )  # Derive vector by intervention unit of staggered counterfactual outcomes standard errors across study periods
                                 
                                 out_ct_ci.lower=out_ct + qnorm(.5 * alpha_ci) * out_ct_se  # Derive vector by intervention unit of staggered counterfactual outcome lower confidence bounds across study periods
                                 
                                 out_ct_ci.upper=out_ct + qnorm(1 - .5 * alpha_ci) * out_ct_se  # Derive vector by intervention unit of staggered counterfactual outcome upper confidence bounds across study periods
                                 
                                 delta_mean=1e2 * att_mean / out_ct  # Derive vector by intervention unit of mean staggered outcome change across study periods
                                 
                                 delta_ci.lower=1e2 * att_ci.lower / out_ct  # Derive vector by intervention unit of staggered outcome change approximate lower confidence bounds across study periods
                                 
                                 delta_ci.upper=1e2 * att_ci.upper / out_ct  # Derive vector by intervention unit of staggered outcome change approximate upper confidence bounds across study periods
                                 }
                             )  # Augment ATT inferences data-frame with second-order and outcome change inferences by intervention unit across study periods
  
  att.out_unit.dat <- with(gsynth_mdl, 
                           expr=within(att.out_unit.dat, 
                                       expr=
                                         {
                                           sigma2=sigma2  # Derive residuals variance
                                           
                                           size_fct=unname(r.cv)  # Derive cross-validated latent factors number
                                           
                                           size_period=c(T0, T - T0)  # Derive vector by intervention unit of study period lengths across study periods
                                           
                                           size_tr=rep(c(0, Ntr), each=Ntr)  # Derive vector by intervention unit of intervention units sample size across study periods
                                           
                                           size_ct=N - size_tr  # Derive vector by intervention unit of control units sample size across study periods
                                           }
                                       )
                           )  # Augment ATT inferences data-frame with model diagnostics by intervention unit across study periods
  
  
  ############################
  ## Derive inferences on   ##
  ## ATT and outcome change ##
  ## by unit and across key ##
  ## study times            ##
  ############################
  
  att.out_unit_mean.dat <- with(gsynth_mdl, 
                                expr=data.frame(comparator=paste(ifelse(min(T0) == max(T0), "old", "new"), "hub", sep="_"), 
                                                outcome=Y,   # Set outcome name
                                                time=rep(c(paste(-time.vec[c(1, length(time.vec))], collapse=" - "), 
                                                           paste(time.vec[c(1, length(time.vec))], collapse=" - ")
                                                           ), 
                                                         each=Ntr
                                                         ), # Set intervention times of interest
                                                unit=id.tr,   # Set intervention units of interest
                                                att_mean=c(colMeans(eff_stg.arr.ls$pre[as.character(-time.vec), ], na.rm=TRUE), 
                                                           colMeans(eff_stg.arr.ls$post[as.character(time.vec), ], na.rm=TRUE)
                                                           ),   # Derive vector by intervention unit of average ATT estimates across key study times
                                                att_se=c(apply(apply(eff_stg_boot.arr.ls$pre[as.character(-time.vec), , ], 
                                                                     MARGIN=3:2, 
                                                                     FUN=mean, na.rm=TRUE
                                                                     ), 
                                                               MARGIN=1, 
                                                               FUN=sd
                                                               ), 
                                                         apply(apply(eff_stg_boot.arr.ls$post[as.character(time.vec), , ], 
                                                                     MARGIN=3:2, 
                                                                     FUN=mean, na.rm=TRUE
                                                                     ), 
                                                               MARGIN=1, 
                                                               FUN=sd
                                                               )
                                                         )  # Derive vector by intervention unit of ATT standard errors across key study times
                                                )
                                )  # Derive data-frame of ATT inferences by intervention unit across key study times
  
  att.out_unit_mean.dat <- within(att.out_unit_mean.dat, 
                                  expr={
                                    att_ci.lower=att_mean + qnorm(.5 * alpha_ci) * att_se  # Derive vector by intervention unit of ATT lower confidence bounds across key study times
                                    
                                    att_ci.upper=att_mean + qnorm(1 - .5 * alpha_ci) * att_se  # Derive vector by intervention unit of ATT upper confidence bounds across key study times
                                    
                                    pval=2 * pnorm(abs(att_mean) / att_se, 
                                                   lower.tail=FALSE
                                                   )  # Derive vector by intervention unit of ATT p-values across key study times
                                    
                                    out_tr=c(colMeans(out_stg.arr.ls$pre[as.character(-time.vec), "tr", ], na.rm=TRUE), 
                                             colMeans(out_stg.arr.ls$post[as.character(time.vec), "tr", ], na.rm=TRUE)
                                             )  # Derive vector by intervention unit of mean staggered intervention outcomes across key study times
                                    
                                    out_ct=c(colMeans(out_stg.arr.ls$pre[as.character(-time.vec), "ct", ], na.rm=TRUE), 
                                             colMeans(out_stg.arr.ls$post[as.character(time.vec), "ct", ], na.rm=TRUE)
                                             )  # Derive vector by intervention unit of mean staggered counterfactual outcomes across key study times
                                    
                                    out_ct_se=c(apply(apply(out.ct_stg_boot.arr.ls$pre[as.character(-time.vec), , ], 
                                                            MARGIN=3:2, 
                                                            FUN=mean, na.rm=TRUE
                                                            ), 
                                                      MARGIN=1, 
                                                      FUN=sd
                                                      ), 
                                                apply(apply(out.ct_stg_boot.arr.ls$post[as.character(time.vec), , ], 
                                                            MARGIN=3:2, 
                                                            FUN=mean, na.rm=TRUE
                                                            ), 
                                                      MARGIN=1, 
                                                      FUN=sd
                                                      )
                                                )  # Derive vector by intervention unit of staggered counterfactual outcomes standard errors across post-intervention times
                                    
                                    out_ct_ci.lower=out_ct + qnorm(.5 * alpha_ci) * out_ct_se  # Derive vector by intervention unit of staggered counterfactual outcome lower confidence bound by key intervention time
                                    
                                    out_ct_ci.upper=out_ct + qnorm(1 - .5 * alpha_ci) * out_ct_se  # Derive vector by intervention unit of staggered counterfactual outcome upper confidence bound by key intervention time
                                    
                                    delta_mean=1e2 * att_mean / out_ct  # Derive vector by intervention unit of mean staggered outcome change across key study times
                                    
                                    delta_ci.lower=1e2 * att_ci.lower / out_ct  # Derive vector by intervention unit of staggered outcome change approximate lower confidence bound across key study times
                                    
                                    delta_ci.upper=1e2 * att_ci.upper / out_ct  # Derive vector by intervention unit of staggered outcome change approximate upper confidence bound across key study times
                                    }
                                  )  # Augment ATT inferences data-frame with second-order and outcome change inferences by intervention unit across key study times
  
  att.out_unit_mean.dat <- with(gsynth_mdl, 
                                expr=within(att.out_unit_mean.dat, 
                                            expr=
                                              {
                                                sigma2=sigma2  # Derive residuals variance
                                                size_fct=unname(r.cv)  # Derive cross-validated latent factors number
                                                
                                                size_period=pmin(c(T0, T - T0), length(time.vec))  # Derive study period length across intervention units
                                                
                                                size_tr=rep(c(0, rowSums(D.tr)[min(T - min(T0) + 1, T - length(time.vec))]), 
                                                            each=Ntr
                                                            )  # Derive vector by study period of intervention units sample size
                                                
                                                size_ct=N - size_tr  # Derive vector by study period of control units sample
                                                }
                                            )
                                )  # Augment ATT inferences data-frame with model diagnostics by intervention unit across key study times
  
  
  ###############################
  ## Derive inferences on      ##
  ## ATT and outcome change    ##
  ## by unit and by study time ##
  ###############################
  
  att.out.dat <- with(gsynth_mdl, 
                      expr=data.frame(comparator=paste(ifelse(min(T0) == max(T0), "old", "new"), "hub", sep="_"), 
                                      outcome=Y,   # Set outcome name
                                      time=c(apply(rbind(T0, T - T0), 
                                                   MARGIN=2, 
                                                   FUN=function(vec) c(-seq.int(vec[1]), seq.int(vec[2]))
                                                   )
                                             ),   # Set intervention times of interest
                                      unit=rep(id.tr, each=T),   # Set intervention units of interest
                                      att_mean=c(sapply(id.tr, 
                                                        FUN=function(trs) 
                                                          c(eff_stg.arr.ls$pre[seq.int(T0[match(trs, table=id.tr)]), trs], 
                                                            eff_stg.arr.ls$post[seq.int(T - T0[match(trs, table=id.tr)]), trs]
                                                            )
                                                        )
                                                 ),   # Derive vector by study time, intervention unit of ATT estimates
                                      att_se=c(sapply(id.tr, 
                                                      FUN=function(trs) 
                                                        c(apply(eff_stg_boot.arr.ls$pre, MARGIN=c(1, 3), FUN=sd)[seq.int(T0[match(trs, table=id.tr)]), trs], 
                                                          apply(eff_stg_boot.arr.ls$post, MARGIN=c(1, 3), FUN=sd)[seq.int(T - T0[match(trs, table=id.tr)]), trs]
                                                          )
                                                      )
                                               )  # Derive vector by study time, intervention unit of ATT standard errors
                                      )
                      )  # Derive data-frame of ATT inferences by study time, intervention unit
  
  att.out.dat <- with(gsynth_mdl, 
                      expr=within(att.out.dat, 
                                  expr=
                                    {
                                      att_ci.lower=att_mean + qnorm(.5 * alpha_ci) * att_se  # Derive vector by study time, intervention unit of ATT lower confidence bounds
                                      
                                      att_ci.upper=att_mean + qnorm(1 - .5 * alpha_ci) * att_se  # Derive vector by study time, intervention unit of ATT upper confidence bounds
                                      
                                      pval=2 * pnorm(abs(att_mean) / att_se, 
                                                     lower.tail=FALSE
                                                     )  # Derive vector by study time, intervention unit of ATT p-values
                                      
                                      out_tr=c(sapply(id.tr, 
                                                      FUN=function(trs) 
                                                        c(out_stg.arr.ls$pre[seq.int(T0[match(trs, table=id.tr)]), "tr", trs], 
                                                          out_stg.arr.ls$post[seq.int(T - T0[match(trs, table=id.tr)]), "tr", trs]
                                                          )
                                                      )
                                               )  # Derive vector by study time, intervention unit of mean staggered intervention outcomes
                                      
                                      out_ct=c(sapply(id.tr, 
                                                      FUN=function(trs) 
                                                        c(out_stg.arr.ls$pre[seq.int(T0[match(trs, table=id.tr)]), "ct", trs], 
                                                          out_stg.arr.ls$post[seq.int(T - T0[match(trs, table=id.tr)]), "ct", trs]
                                                          )
                                                      )
                                               )  # Derive vector by study time, intervention unit of mean staggered counterfactual outcomes
                                      
                                      out_ct_se=c(sapply(id.tr, 
                                                         FUN=function(trs) 
                                                           c(apply(out.ct_stg_boot.arr.ls$pre, MARGIN=c(1, 3), FUN=sd)[seq.int(T0[match(trs, table=id.tr)]), trs], 
                                                             apply(out.ct_stg_boot.arr.ls$post, MARGIN=c(1, 3), FUN=sd)[seq.int(T - T0[match(trs, table=id.tr)]), trs]
                                                             )
                                                         )
                                                  )  # Derive vector by post-intervention time, intervention unit of staggered counterfactual outcomes standard errors
                                      
                                      out_ct_ci.lower=out_ct + qnorm(.5 * alpha_ci) * out_ct_se  # Derive vector by study time, intervention unit of staggered counterfactual outcome lower confidence bound
                                      
                                      out_ct_ci.upper=out_ct + qnorm(1 - .5 * alpha_ci) * out_ct_se  # Derive vector by study time, intervention unit of staggered counterfactual outcome upper confidence bound
                                      
                                      delta_mean=1e2 * att_mean / out_ct  # Derive vector by study time, intervention unit of mean staggered outcome change
                                      
                                      delta_ci.lower=1e2 * att_ci.lower / out_ct  # Derive vector by study time, intervention unit of staggered outcome change approximate lower confidence bound
                                      
                                      delta_ci.upper=1e2 * att_ci.upper / out_ct  # Derive vector by study time, intervention unit of staggered outcome change approximate upper confidence bound
                                      }
                                  )
                      )  # Augment ATT inferences data-frame with second-order and outcome change inferences by study time, intervention unit
  
  att.out.dat <- with(gsynth_mdl, 
                      expr=within(att.out.dat, 
                                  expr=
                                    {
                                      sigma2=sigma2  # Derive residuals variance
                                      
                                      size_fct=unname(r.cv)  # Derive cross-validated latent factors number
                                      
                                      size_period=c(apply(rbind(T0, T - T0), 
                                                   MARGIN=2, 
                                                   FUN=function(vec) rep.int(vec, times=vec))
                                             )  # Derive vector by study time, intervention unit of study period lengths
                                      
                                      size_tr=rep.int(c(rowSums(D.tr)[seq.int(min(T0))], 
                                                        rev(rowSums(D.tr)[- seq.int(min(T0))])
                                                        ), 
                                                      times=Ntr
                                                      )  # Derive vector by study time, intervention unit of intervention units sample sizes
                                      
                                      size_ct=N - size_tr  # Derive vector by study time, intervention unit of control units sample sizes
                                      }
                                  )
                      )  # Augment ATT inferences data-frame with model diagnostics by study time, intervention unit
  
  
  #############################
  ## Tabulate ATT and change ##
  ## inferences as output    ##
  #############################
  
  inf.dat <- rbind(att.out_mean.dat, att.out_time.dat, att.out_time_mean.dat, 
                   att.out_unit.dat, att.out_unit_mean.dat, att.out.dat, 
                   make.row.names=FALSE
                   )  # Bind ATT and change inferences data-frames
  
  inf.dat <- inf.dat[, c("comparator", "outcome", "time", "unit", 
                         paste("att", 
                               c("mean", "se", "ci.lower", "ci.upper"), 
                               sep="_"
                               ), 
                         "pval", 
                         paste("out", c("tr", "ct"), sep="_"), 
                         paste("out_ct", c("se", "ci.lower", "ci.upper"), 
                               sep="_"
                               ), 
                         paste("delta", c("mean", "ci.lower", "ci.upper"), sep="_"), 
                         "sigma2", paste("size", c("fct", "period", "tr", "ct"), sep="_")
                         )]  # Re-order inferences data-frame columns
  
  
  ###########################
  ## Export ATT and change ##
  ## inferences to disc    ##
  ###########################
  
  if(file.exists(dest.dir)) 
    write.table(inf.dat, 
                file=file.path(dest.dir, 
                               paste(out.prf, "att.delta.csv", sep="_")
                               ), 
                quote=TRUE, sep=",", row.names=FALSE, col.names=TRUE
                )  # Export ATT and change inferences data-frame as .csv file
  
  
  return(inf.dat)  # Return as output ATT and change inferences data-frame
  }