knitr::opts_chunk$set(echo = TRUE)
Package Invoking
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.2 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.2 ✔ tibble 3.2.1
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(dplyr)
library(zoo)
##
## Attaching package: 'zoo'
##
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(data.table)
##
## Attaching package: 'data.table'
##
## The following objects are masked from 'package:lubridate':
##
## hour, isoweek, mday, minute, month, quarter, second, wday, week,
## yday, year
##
## The following objects are masked from 'package:dplyr':
##
## between, first, last
##
## The following object is masked from 'package:purrr':
##
## transpose
library(panelView)
## ## See bit.ly/panelview4r for more info.
## ## Report bugs -> yiqingxu@stanford.edu.
library(bpCausal)
library(gridExtra)
##
## Attaching package: 'gridExtra'
##
## The following object is masked from 'package:dplyr':
##
## combine
Data Loading
TSCS_Y_filter <- read.csv("/Users/apple/Quantitative\ Marketing\ Research/Statistical\ Modeling\ IX/Statistical\ Modeling\ IX\ Data/TSCS_Y_filter.csv")
TSCS_Y <- read.csv("/Users/apple/Quantitative\ Marketing\ Research/Statistical\ Modeling\ IX/Statistical\ Modeling\ IX\ Data/TSCS_Y.csv")
Formulation of
Hypothetical Treated Units in Complete Data Matrix
TSCS_Y <- as.data.table(TSCS_Y)
TSCS_Y_filter <- as.data.table(TSCS_Y_filter)
TSCS_Y_D0 <- TSCS_Y[D == 0]
D0_id_count <- TSCS_Y_D0[, .(count = .N), by = id]
subset_ids_gt53 <- D0_id_count[count >= 53, .(id)] # 95 IDs
TSCS_Y_filter_gt53 <- TSCS_Y_filter[id %in% subset_ids_gt53$id]
panelview(D = "D",
data = TSCS_Y_filter_gt53,
index = c("id", "T"),
display.all = TRUE, type = "treat", by.timing = TRUE,
xlab = "Time (Biweek)", ylab = "Pair ID",
axis.lab.gap = c(2,10), main = "Treatment Status (Subset of Active Periods >= 53)")
TSCS_Y_filter_gt53_hypo <- TSCS_Y_filter_gt53[!(T >= 54 & T <= 55)]
TSCS_Y_filter_gt53_removed <- TSCS_Y_filter_gt53[T >= 54 & T <= 55]
set.seed(123)
hypo_synth_trt <- TSCS_Y_filter_gt53_hypo
unique_ids <- unique(hypo_synth_trt$id)
n_ids <- length(unique_ids)
n_30p <- round(0.30 * n_ids)
n_20p <- round(0.20 * n_ids)
n_10p <- round(0.10 * n_ids)
n_5p <- round(0.05 * n_ids)
shuffled_ids <- sample(unique_ids)
ids_30p <- shuffled_ids[1:n_30p]
ids_20p <- shuffled_ids[(n_30p + 1):(n_30p + n_20p)]
ids_10p <- shuffled_ids[(n_30p + n_20p + 1):(n_30p + n_20p + n_10p)]
ids_5p <- shuffled_ids[(n_30p + n_20p + n_10p + 1):(n_30p + n_20p + n_10p + n_5p)]
hypo_synth_trt[id %in% ids_30p & T >= 36, D := 1]
hypo_synth_trt[id %in% ids_20p & T >= 31, D := 1]
hypo_synth_trt[id %in% ids_10p & T >= 26, D := 1]
hypo_synth_trt[id %in% ids_5p & T >= 21, D := 1]
head(hypo_synth_trt, 10)
## id T D F_trans M_trans R_trans F_redem M_redem R_redem F_charg M_charg
## 1: 2779 1 0 1 350 733 0 0.00 5 1 250
## 2: 2779 2 0 1 350 733 0 0.00 5 1 250
## 3: 2779 3 0 1 350 733 0 0.00 5 1 250
## 4: 2779 4 0 1 350 733 0 0.00 5 1 250
## 5: 2779 5 0 1 350 733 1 10.29 674 1 250
## 6: 2779 6 0 1 350 733 3 82.82 661 1 250
## 7: 2779 7 0 1 350 733 3 82.82 661 1 250
## 8: 2779 8 0 1 350 733 3 82.82 661 1 250
## 9: 2779 9 0 1 350 733 3 82.82 661 1 250
## 10: 2779 10 0 1 350 733 3 82.82 661 1 250
## R_charg F_tip M_tip user proj
## 1: 733 0 0.00 30 10
## 2: 733 0 0.00 30 10
## 3: 733 0 0.00 30 10
## 4: 733 0 0.00 30 10
## 5: 733 1 2.26 30 10
## 6: 733 3 18.22 30 10
## 7: 733 3 18.22 30 10
## 8: 733 3 18.22 30 10
## 9: 733 3 18.22 30 10
## 10: 733 3 18.22 30 10
panelview(D = "D",
data = hypo_synth_trt,
index = c("id", "T"),
display.all = TRUE, type = "treat", by.timing = TRUE,
xlab = "Time (Biweek)", ylab = "Pair ID",
axis.lab.gap = c(2,10), main = "Treatment Status with Hypothetical Treated Units")
Trend of Unit-Level
Time-Series Visualization
Transaction
grid.arrange(transM, transF, transR, ncol = 1)
Redemption
grid.arrange(redemM, redemF, redemR, ncol = 1)
Bayesian Hierarchical
Dynamic Multi-Outcome (BHDM) Latent Factor Model
Implementation of
External C++ Codes
Rcpp::sourceCpp("/Users/apple/Quantitative\ Marketing\ Research/Statistical\ Modeling\ VI/Statistical\ Modeling\ VI\ Data/bpCausal-main/src/blasso.cpp")
Implementation of
External R Scripts
source("/Users/apple/Quantitative\ Marketing\ Research/Statistical\ Modeling\ VI/Statistical\ Modeling\ VI\ Data/bpCausal-main/R/blasso_default.R")
source("/Users/apple/Quantitative\ Marketing\ Research/Statistical\ Modeling\ VI/Statistical\ Modeling\ VI\ Data/bpCausal-main/R/blasso_core.R")
Implementation of
BHDM
BHDM <- function(data, index, Yname_vector, Dname, Xname, Zname, Aname,
re, ar1, r, niter = 15000, burn = 5000,
xlasso = 1, zlasso = 1, alasso = 1, flasso = 1,
a1 = 0.001, a2 = 0.001, b1 = 0.001, b2 = 0.001,
c1 = 0.001, c2 = 0.001, p1 = 0.001, p2 = 0.001) {
out <- lapply(Yname_vector, function(Yname_single) {
bpCausal(data = data,
index = index,
Yname = Yname_single,
Dname = Dname,
Xname = Xname,
Zname = Zname,
Aname = Aname,
re = re,
ar1 = ar1,
r = r,
niter = niter,
burn = burn,
xlasso = xlasso,
zlasso = zlasso,
alasso = alasso,
flasso = flasso,
a1 = a1, a2 = a2,
b1 = b1, b2 = b2,
c1 = c1, c2 = c2,
p1 = p1, p2 = p2)
})
return(out)
}
Posteriors
OUT <-
BHDM(data = hypo_synth_trt, index = c("id", "T"),
Yname_vector = c("F_trans", "M_trans", "R_trans", "F_redem", "M_redem", "R_redem"),
Dname = "D", Xname = c(), Zname = c(), Aname = c(),
re = "both", ar1 = TRUE, r = 10)
## Simulated times: 100
## Simulated times: 200
## Simulated times: 300
## Simulated times: 400
## Simulated times: 500
## Simulated times: 600
## Simulated times: 700
## Simulated times: 800
## Simulated times: 900
## Simulated times: 1000
## Simulated times: 1100
## Simulated times: 1200
## Simulated times: 1300
## Simulated times: 1400
## Simulated times: 1500
## Simulated times: 1600
## Simulated times: 1700
## Simulated times: 1800
## Simulated times: 1900
## Simulated times: 2000
## Simulated times: 2100
## Simulated times: 2200
## Simulated times: 2300
## Simulated times: 2400
## Simulated times: 2500
## Simulated times: 2600
## Simulated times: 2700
## Simulated times: 2800
## Simulated times: 2900
## Simulated times: 3000
## Simulated times: 3100
## Simulated times: 3200
## Simulated times: 3300
## Simulated times: 3400
## Simulated times: 3500
## Simulated times: 3600
## Simulated times: 3700
## Simulated times: 3800
## Simulated times: 3900
## Simulated times: 4000
## Simulated times: 4100
## Simulated times: 4200
## Simulated times: 4300
## Simulated times: 4400
## Simulated times: 4500
## Simulated times: 4600
## Simulated times: 4700
## Simulated times: 4800
## Simulated times: 4900
## Simulated times: 5000
## Simulated times: 5100
## Simulated times: 5200
## Simulated times: 5300
## Simulated times: 5400
## Simulated times: 5500
## Simulated times: 5600
## Simulated times: 5700
## Simulated times: 5800
## Simulated times: 5900
## Simulated times: 6000
## Simulated times: 6100
## Simulated times: 6200
## Simulated times: 6300
## Simulated times: 6400
## Simulated times: 6500
## Simulated times: 6600
## Simulated times: 6700
## Simulated times: 6800
## Simulated times: 6900
## Simulated times: 7000
## Simulated times: 7100
## Simulated times: 7200
## Simulated times: 7300
## Simulated times: 7400
## Simulated times: 7500
## Simulated times: 7600
## Simulated times: 7700
## Simulated times: 7800
## Simulated times: 7900
## Simulated times: 8000
## Simulated times: 8100
## Simulated times: 8200
## Simulated times: 8300
## Simulated times: 8400
## Simulated times: 8500
## Simulated times: 8600
## Simulated times: 8700
## Simulated times: 8800
## Simulated times: 8900
## Simulated times: 9000
## Simulated times: 9100
## Simulated times: 9200
## Simulated times: 9300
## Simulated times: 9400
## Simulated times: 9500
## Simulated times: 9600
## Simulated times: 9700
## Simulated times: 9800
## Simulated times: 9900
## Simulated times: 10000
## Simulated times: 10100
## Simulated times: 10200
## Simulated times: 10300
## Simulated times: 10400
## Simulated times: 10500
## Simulated times: 10600
## Simulated times: 10700
## Simulated times: 10800
## Simulated times: 10900
## Simulated times: 11000
## Simulated times: 11100
## Simulated times: 11200
## Simulated times: 11300
## Simulated times: 11400
## Simulated times: 11500
## Simulated times: 11600
## Simulated times: 11700
## Simulated times: 11800
## Simulated times: 11900
## Simulated times: 12000
## Simulated times: 12100
## Simulated times: 12200
## Simulated times: 12300
## Simulated times: 12400
## Simulated times: 12500
## Simulated times: 12600
## Simulated times: 12700
## Simulated times: 12800
## Simulated times: 12900
## Simulated times: 13000
## Simulated times: 13100
## Simulated times: 13200
## Simulated times: 13300
## Simulated times: 13400
## Simulated times: 13500
## Simulated times: 13600
## Simulated times: 13700
## Simulated times: 13800
## Simulated times: 13900
## Simulated times: 14000
## Simulated times: 14100
## Simulated times: 14200
## Simulated times: 14300
## Simulated times: 14400
## Simulated times: 14500
## Simulated times: 14600
## Simulated times: 14700
## Simulated times: 14800
## Simulated times: 14900
## Simulated times: 15000
## Simulated times: 100
## Simulated times: 200
## Simulated times: 300
## Simulated times: 400
## Simulated times: 500
## Simulated times: 600
## Simulated times: 700
## Simulated times: 800
## Simulated times: 900
## Simulated times: 1000
## Simulated times: 1100
## Simulated times: 1200
## Simulated times: 1300
## Simulated times: 1400
## Simulated times: 1500
## Simulated times: 1600
## Simulated times: 1700
## Simulated times: 1800
## Simulated times: 1900
## Simulated times: 2000
## Simulated times: 2100
## Simulated times: 2200
## Simulated times: 2300
## Simulated times: 2400
## Simulated times: 2500
## Simulated times: 2600
## Simulated times: 2700
## Simulated times: 2800
## Simulated times: 2900
## Simulated times: 3000
## Simulated times: 3100
## Simulated times: 3200
## Simulated times: 3300
## Simulated times: 3400
## Simulated times: 3500
## Simulated times: 3600
## Simulated times: 3700
## Simulated times: 3800
## Simulated times: 3900
## Simulated times: 4000
## Simulated times: 4100
## Simulated times: 4200
## Simulated times: 4300
## Simulated times: 4400
## Simulated times: 4500
## Simulated times: 4600
## Simulated times: 4700
## Simulated times: 4800
## Simulated times: 4900
## Simulated times: 5000
## Simulated times: 5100
## Simulated times: 5200
## Simulated times: 5300
## Simulated times: 5400
## Simulated times: 5500
## Simulated times: 5600
## Simulated times: 5700
## Simulated times: 5800
## Simulated times: 5900
## Simulated times: 6000
## Simulated times: 6100
## Simulated times: 6200
## Simulated times: 6300
## Simulated times: 6400
## Simulated times: 6500
## Simulated times: 6600
## Simulated times: 6700
## Simulated times: 6800
## Simulated times: 6900
## Simulated times: 7000
## Simulated times: 7100
## Simulated times: 7200
## Simulated times: 7300
## Simulated times: 7400
## Simulated times: 7500
## Simulated times: 7600
## Simulated times: 7700
## Simulated times: 7800
## Simulated times: 7900
## Simulated times: 8000
## Simulated times: 8100
## Simulated times: 8200
## Simulated times: 8300
## Simulated times: 8400
## Simulated times: 8500
## Simulated times: 8600
## Simulated times: 8700
## Simulated times: 8800
## Simulated times: 8900
## Simulated times: 9000
## Simulated times: 9100
## Simulated times: 9200
## Simulated times: 9300
## Simulated times: 9400
## Simulated times: 9500
## Simulated times: 9600
## Simulated times: 9700
## Simulated times: 9800
## Simulated times: 9900
## Simulated times: 10000
## Simulated times: 10100
## Simulated times: 10200
## Simulated times: 10300
## Simulated times: 10400
## Simulated times: 10500
## Simulated times: 10600
## Simulated times: 10700
## Simulated times: 10800
## Simulated times: 10900
## Simulated times: 11000
## Simulated times: 11100
## Simulated times: 11200
## Simulated times: 11300
## Simulated times: 11400
## Simulated times: 11500
## Simulated times: 11600
## Simulated times: 11700
## Simulated times: 11800
## Simulated times: 11900
## Simulated times: 12000
## Simulated times: 12100
## Simulated times: 12200
## Simulated times: 12300
## Simulated times: 12400
## Simulated times: 12500
## Simulated times: 12600
## Simulated times: 12700
## Simulated times: 12800
## Simulated times: 12900
## Simulated times: 13000
## Simulated times: 13100
## Simulated times: 13200
## Simulated times: 13300
## Simulated times: 13400
## Simulated times: 13500
## Simulated times: 13600
## Simulated times: 13700
## Simulated times: 13800
## Simulated times: 13900
## Simulated times: 14000
## Simulated times: 14100
## Simulated times: 14200
## Simulated times: 14300
## Simulated times: 14400
## Simulated times: 14500
## Simulated times: 14600
## Simulated times: 14700
## Simulated times: 14800
## Simulated times: 14900
## Simulated times: 15000
## Simulated times: 100
## Simulated times: 200
## Simulated times: 300
## Simulated times: 400
## Simulated times: 500
## Simulated times: 600
## Simulated times: 700
## Simulated times: 800
## Simulated times: 900
## Simulated times: 1000
## Simulated times: 1100
## Simulated times: 1200
## Simulated times: 1300
## Simulated times: 1400
## Simulated times: 1500
## Simulated times: 1600
## Simulated times: 1700
## Simulated times: 1800
## Simulated times: 1900
## Simulated times: 2000
## Simulated times: 2100
## Simulated times: 2200
## Simulated times: 2300
## Simulated times: 2400
## Simulated times: 2500
## Simulated times: 2600
## Simulated times: 2700
## Simulated times: 2800
## Simulated times: 2900
## Simulated times: 3000
## Simulated times: 3100
## Simulated times: 3200
## Simulated times: 3300
## Simulated times: 3400
## Simulated times: 3500
## Simulated times: 3600
## Simulated times: 3700
## Simulated times: 3800
## Simulated times: 3900
## Simulated times: 4000
## Simulated times: 4100
## Simulated times: 4200
## Simulated times: 4300
## Simulated times: 4400
## Simulated times: 4500
## Simulated times: 4600
## Simulated times: 4700
## Simulated times: 4800
## Simulated times: 4900
## Simulated times: 5000
## Simulated times: 5100
## Simulated times: 5200
## Simulated times: 5300
## Simulated times: 5400
## Simulated times: 5500
## Simulated times: 5600
## Simulated times: 5700
## Simulated times: 5800
## Simulated times: 5900
## Simulated times: 6000
## Simulated times: 6100
## Simulated times: 6200
## Simulated times: 6300
## Simulated times: 6400
## Simulated times: 6500
## Simulated times: 6600
## Simulated times: 6700
## Simulated times: 6800
## Simulated times: 6900
## Simulated times: 7000
## Simulated times: 7100
## Simulated times: 7200
## Simulated times: 7300
## Simulated times: 7400
## Simulated times: 7500
## Simulated times: 7600
## Simulated times: 7700
## Simulated times: 7800
## Simulated times: 7900
## Simulated times: 8000
## Simulated times: 8100
## Simulated times: 8200
## Simulated times: 8300
## Simulated times: 8400
## Simulated times: 8500
## Simulated times: 8600
## Simulated times: 8700
## Simulated times: 8800
## Simulated times: 8900
## Simulated times: 9000
## Simulated times: 9100
## Simulated times: 9200
## Simulated times: 9300
## Simulated times: 9400
## Simulated times: 9500
## Simulated times: 9600
## Simulated times: 9700
## Simulated times: 9800
## Simulated times: 9900
## Simulated times: 10000
## Simulated times: 10100
## Simulated times: 10200
## Simulated times: 10300
## Simulated times: 10400
## Simulated times: 10500
## Simulated times: 10600
## Simulated times: 10700
## Simulated times: 10800
## Simulated times: 10900
## Simulated times: 11000
## Simulated times: 11100
## Simulated times: 11200
## Simulated times: 11300
## Simulated times: 11400
## Simulated times: 11500
## Simulated times: 11600
## Simulated times: 11700
## Simulated times: 11800
## Simulated times: 11900
## Simulated times: 12000
## Simulated times: 12100
## Simulated times: 12200
## Simulated times: 12300
## Simulated times: 12400
## Simulated times: 12500
## Simulated times: 12600
## Simulated times: 12700
## Simulated times: 12800
## Simulated times: 12900
## Simulated times: 13000
## Simulated times: 13100
## Simulated times: 13200
## Simulated times: 13300
## Simulated times: 13400
## Simulated times: 13500
## Simulated times: 13600
## Simulated times: 13700
## Simulated times: 13800
## Simulated times: 13900
## Simulated times: 14000
## Simulated times: 14100
## Simulated times: 14200
## Simulated times: 14300
## Simulated times: 14400
## Simulated times: 14500
## Simulated times: 14600
## Simulated times: 14700
## Simulated times: 14800
## Simulated times: 14900
## Simulated times: 15000
## Simulated times: 100
## Simulated times: 200
## Simulated times: 300
## Simulated times: 400
## Simulated times: 500
## Simulated times: 600
## Simulated times: 700
## Simulated times: 800
## Simulated times: 900
## Simulated times: 1000
## Simulated times: 1100
## Simulated times: 1200
## Simulated times: 1300
## Simulated times: 1400
## Simulated times: 1500
## Simulated times: 1600
## Simulated times: 1700
## Simulated times: 1800
## Simulated times: 1900
## Simulated times: 2000
## Simulated times: 2100
## Simulated times: 2200
## Simulated times: 2300
## Simulated times: 2400
## Simulated times: 2500
## Simulated times: 2600
## Simulated times: 2700
## Simulated times: 2800
## Simulated times: 2900
## Simulated times: 3000
## Simulated times: 3100
## Simulated times: 3200
## Simulated times: 3300
## Simulated times: 3400
## Simulated times: 3500
## Simulated times: 3600
## Simulated times: 3700
## Simulated times: 3800
## Simulated times: 3900
## Simulated times: 4000
## Simulated times: 4100
## Simulated times: 4200
## Simulated times: 4300
## Simulated times: 4400
## Simulated times: 4500
## Simulated times: 4600
## Simulated times: 4700
## Simulated times: 4800
## Simulated times: 4900
## Simulated times: 5000
## Simulated times: 5100
## Simulated times: 5200
## Simulated times: 5300
## Simulated times: 5400
## Simulated times: 5500
## Simulated times: 5600
## Simulated times: 5700
## Simulated times: 5800
## Simulated times: 5900
## Simulated times: 6000
## Simulated times: 6100
## Simulated times: 6200
## Simulated times: 6300
## Simulated times: 6400
## Simulated times: 6500
## Simulated times: 6600
## Simulated times: 6700
## Simulated times: 6800
## Simulated times: 6900
## Simulated times: 7000
## Simulated times: 7100
## Simulated times: 7200
## Simulated times: 7300
## Simulated times: 7400
## Simulated times: 7500
## Simulated times: 7600
## Simulated times: 7700
## Simulated times: 7800
## Simulated times: 7900
## Simulated times: 8000
## Simulated times: 8100
## Simulated times: 8200
## Simulated times: 8300
## Simulated times: 8400
## Simulated times: 8500
## Simulated times: 8600
## Simulated times: 8700
## Simulated times: 8800
## Simulated times: 8900
## Simulated times: 9000
## Simulated times: 9100
## Simulated times: 9200
## Simulated times: 9300
## Simulated times: 9400
## Simulated times: 9500
## Simulated times: 9600
## Simulated times: 9700
## Simulated times: 9800
## Simulated times: 9900
## Simulated times: 10000
## Simulated times: 10100
## Simulated times: 10200
## Simulated times: 10300
## Simulated times: 10400
## Simulated times: 10500
## Simulated times: 10600
## Simulated times: 10700
## Simulated times: 10800
## Simulated times: 10900
## Simulated times: 11000
## Simulated times: 11100
## Simulated times: 11200
## Simulated times: 11300
## Simulated times: 11400
## Simulated times: 11500
## Simulated times: 11600
## Simulated times: 11700
## Simulated times: 11800
## Simulated times: 11900
## Simulated times: 12000
## Simulated times: 12100
## Simulated times: 12200
## Simulated times: 12300
## Simulated times: 12400
## Simulated times: 12500
## Simulated times: 12600
## Simulated times: 12700
## Simulated times: 12800
## Simulated times: 12900
## Simulated times: 13000
## Simulated times: 13100
## Simulated times: 13200
## Simulated times: 13300
## Simulated times: 13400
## Simulated times: 13500
## Simulated times: 13600
## Simulated times: 13700
## Simulated times: 13800
## Simulated times: 13900
## Simulated times: 14000
## Simulated times: 14100
## Simulated times: 14200
## Simulated times: 14300
## Simulated times: 14400
## Simulated times: 14500
## Simulated times: 14600
## Simulated times: 14700
## Simulated times: 14800
## Simulated times: 14900
## Simulated times: 15000
## Simulated times: 100
## Simulated times: 200
## Simulated times: 300
## Simulated times: 400
## Simulated times: 500
## Simulated times: 600
## Simulated times: 700
## Simulated times: 800
## Simulated times: 900
## Simulated times: 1000
## Simulated times: 1100
## Simulated times: 1200
## Simulated times: 1300
## Simulated times: 1400
## Simulated times: 1500
## Simulated times: 1600
## Simulated times: 1700
## Simulated times: 1800
## Simulated times: 1900
## Simulated times: 2000
## Simulated times: 2100
## Simulated times: 2200
## Simulated times: 2300
## Simulated times: 2400
## Simulated times: 2500
## Simulated times: 2600
## Simulated times: 2700
## Simulated times: 2800
## Simulated times: 2900
## Simulated times: 3000
## Simulated times: 3100
## Simulated times: 3200
## Simulated times: 3300
## Simulated times: 3400
## Simulated times: 3500
## Simulated times: 3600
## Simulated times: 3700
## Simulated times: 3800
## Simulated times: 3900
## Simulated times: 4000
## Simulated times: 4100
## Simulated times: 4200
## Simulated times: 4300
## Simulated times: 4400
## Simulated times: 4500
## Simulated times: 4600
## Simulated times: 4700
## Simulated times: 4800
## Simulated times: 4900
## Simulated times: 5000
## Simulated times: 5100
## Simulated times: 5200
## Simulated times: 5300
## Simulated times: 5400
## Simulated times: 5500
## Simulated times: 5600
## Simulated times: 5700
## Simulated times: 5800
## Simulated times: 5900
## Simulated times: 6000
## Simulated times: 6100
## Simulated times: 6200
## Simulated times: 6300
## Simulated times: 6400
## Simulated times: 6500
## Simulated times: 6600
## Simulated times: 6700
## Simulated times: 6800
## Simulated times: 6900
## Simulated times: 7000
## Simulated times: 7100
## Simulated times: 7200
## Simulated times: 7300
## Simulated times: 7400
## Simulated times: 7500
## Simulated times: 7600
## Simulated times: 7700
## Simulated times: 7800
## Simulated times: 7900
## Simulated times: 8000
## Simulated times: 8100
## Simulated times: 8200
## Simulated times: 8300
## Simulated times: 8400
## Simulated times: 8500
## Simulated times: 8600
## Simulated times: 8700
## Simulated times: 8800
## Simulated times: 8900
## Simulated times: 9000
## Simulated times: 9100
## Simulated times: 9200
## Simulated times: 9300
## Simulated times: 9400
## Simulated times: 9500
## Simulated times: 9600
## Simulated times: 9700
## Simulated times: 9800
## Simulated times: 9900
## Simulated times: 10000
## Simulated times: 10100
## Simulated times: 10200
## Simulated times: 10300
## Simulated times: 10400
## Simulated times: 10500
## Simulated times: 10600
## Simulated times: 10700
## Simulated times: 10800
## Simulated times: 10900
## Simulated times: 11000
## Simulated times: 11100
## Simulated times: 11200
## Simulated times: 11300
## Simulated times: 11400
## Simulated times: 11500
## Simulated times: 11600
## Simulated times: 11700
## Simulated times: 11800
## Simulated times: 11900
## Simulated times: 12000
## Simulated times: 12100
## Simulated times: 12200
## Simulated times: 12300
## Simulated times: 12400
## Simulated times: 12500
## Simulated times: 12600
## Simulated times: 12700
## Simulated times: 12800
## Simulated times: 12900
## Simulated times: 13000
## Simulated times: 13100
## Simulated times: 13200
## Simulated times: 13300
## Simulated times: 13400
## Simulated times: 13500
## Simulated times: 13600
## Simulated times: 13700
## Simulated times: 13800
## Simulated times: 13900
## Simulated times: 14000
## Simulated times: 14100
## Simulated times: 14200
## Simulated times: 14300
## Simulated times: 14400
## Simulated times: 14500
## Simulated times: 14600
## Simulated times: 14700
## Simulated times: 14800
## Simulated times: 14900
## Simulated times: 15000
## Simulated times: 100
## Simulated times: 200
## Simulated times: 300
## Simulated times: 400
## Simulated times: 500
## Simulated times: 600
## Simulated times: 700
## Simulated times: 800
## Simulated times: 900
## Simulated times: 1000
## Simulated times: 1100
## Simulated times: 1200
## Simulated times: 1300
## Simulated times: 1400
## Simulated times: 1500
## Simulated times: 1600
## Simulated times: 1700
## Simulated times: 1800
## Simulated times: 1900
## Simulated times: 2000
## Simulated times: 2100
## Simulated times: 2200
## Simulated times: 2300
## Simulated times: 2400
## Simulated times: 2500
## Simulated times: 2600
## Simulated times: 2700
## Simulated times: 2800
## Simulated times: 2900
## Simulated times: 3000
## Simulated times: 3100
## Simulated times: 3200
## Simulated times: 3300
## Simulated times: 3400
## Simulated times: 3500
## Simulated times: 3600
## Simulated times: 3700
## Simulated times: 3800
## Simulated times: 3900
## Simulated times: 4000
## Simulated times: 4100
## Simulated times: 4200
## Simulated times: 4300
## Simulated times: 4400
## Simulated times: 4500
## Simulated times: 4600
## Simulated times: 4700
## Simulated times: 4800
## Simulated times: 4900
## Simulated times: 5000
## Simulated times: 5100
## Simulated times: 5200
## Simulated times: 5300
## Simulated times: 5400
## Simulated times: 5500
## Simulated times: 5600
## Simulated times: 5700
## Simulated times: 5800
## Simulated times: 5900
## Simulated times: 6000
## Simulated times: 6100
## Simulated times: 6200
## Simulated times: 6300
## Simulated times: 6400
## Simulated times: 6500
## Simulated times: 6600
## Simulated times: 6700
## Simulated times: 6800
## Simulated times: 6900
## Simulated times: 7000
## Simulated times: 7100
## Simulated times: 7200
## Simulated times: 7300
## Simulated times: 7400
## Simulated times: 7500
## Simulated times: 7600
## Simulated times: 7700
## Simulated times: 7800
## Simulated times: 7900
## Simulated times: 8000
## Simulated times: 8100
## Simulated times: 8200
## Simulated times: 8300
## Simulated times: 8400
## Simulated times: 8500
## Simulated times: 8600
## Simulated times: 8700
## Simulated times: 8800
## Simulated times: 8900
## Simulated times: 9000
## Simulated times: 9100
## Simulated times: 9200
## Simulated times: 9300
## Simulated times: 9400
## Simulated times: 9500
## Simulated times: 9600
## Simulated times: 9700
## Simulated times: 9800
## Simulated times: 9900
## Simulated times: 10000
## Simulated times: 10100
## Simulated times: 10200
## Simulated times: 10300
## Simulated times: 10400
## Simulated times: 10500
## Simulated times: 10600
## Simulated times: 10700
## Simulated times: 10800
## Simulated times: 10900
## Simulated times: 11000
## Simulated times: 11100
## Simulated times: 11200
## Simulated times: 11300
## Simulated times: 11400
## Simulated times: 11500
## Simulated times: 11600
## Simulated times: 11700
## Simulated times: 11800
## Simulated times: 11900
## Simulated times: 12000
## Simulated times: 12100
## Simulated times: 12200
## Simulated times: 12300
## Simulated times: 12400
## Simulated times: 12500
## Simulated times: 12600
## Simulated times: 12700
## Simulated times: 12800
## Simulated times: 12900
## Simulated times: 13000
## Simulated times: 13100
## Simulated times: 13200
## Simulated times: 13300
## Simulated times: 13400
## Simulated times: 13500
## Simulated times: 13600
## Simulated times: 13700
## Simulated times: 13800
## Simulated times: 13900
## Simulated times: 14000
## Simulated times: 14100
## Simulated times: 14200
## Simulated times: 14300
## Simulated times: 14400
## Simulated times: 14500
## Simulated times: 14600
## Simulated times: 14700
## Simulated times: 14800
## Simulated times: 14900
## Simulated times: 15000
Counterfactual
Estimate
counterfactual_est <- function(x) {
niter <- dim(x$sigma2)[2]
yct_i <- x$yct
yct_i <- matrix(c(yct_i[, (1):niter]), nrow(yct_i), niter)
yo_t <- x$yo_t
id_tr <- x$raw.id.tr
time_tr = x$time.tr
m_yct_mean <- apply(yct_i, 1, mean)
m_yct_ci_l <- apply(yct_i, 1, quantile, 0.025)
m_yct_ci_u <- apply(yct_i, 1, quantile, 0.975)
result_x <- data.frame(
id = id_tr,
T = as.integer(time_tr),
original_outcome = yo_t,
counterfactual_estimate = m_yct_mean,
ci_lower = m_yct_ci_l,
ci_upper = m_yct_ci_u
)
return(result_x)
}
counterfactual_results <- lapply(OUT, counterfactual_est)
Posterior
Distributions
hypo_synth_trt_post <- hypo_synth_trt |> group_by(id) |>
mutate(is_treat = ifelse(any(D == 1), 1, 0)) |> ungroup()
hypo_synth_trt_post <- hypo_synth_trt_post |> arrange(id, T) |> group_by(id) |>
mutate(T0 = case_when(any(D == 1) ~ min(T[D == 1]) - 1, TRUE ~ max(T))) |> ungroup()
## Warning: There were 33 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `T0 = case_when(any(D == 1) ~ min(T[D == 1]) - 1, TRUE ~
## max(T))`.
## ℹ In group 1: `id = 2779`.
## Caused by warning in `min()`:
## ! no non-missing arguments to min; returning Inf
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 32 remaining warnings.
hypo_synth_trt_post <- hypo_synth_trt_post[, c(1:3, 17, 18, 15:16, 4:14)]
hypo_synth_trt_post$is_treat <- as.integer(hypo_synth_trt_post$is_treat)
hypo_synth_trt_post$T0 <- as.integer(hypo_synth_trt_post$T0)
F_T <- hypo_synth_trt_post[, c(1:7, 8)]
M_T <- hypo_synth_trt_post[, c(1:7, 9)]
R_T <- hypo_synth_trt_post[, c(1:7, 10)]
F_R <- hypo_synth_trt_post[, c(1:7, 11)]
M_R <- hypo_synth_trt_post[, c(1:7, 12)]
R_R <- hypo_synth_trt_post[, c(1:7, 13)]
FT <- full_join(F_T, counterfactual_results[[1]], by = c("id", "T", "F_trans" = "original_outcome"))
colnames(FT)[c(9:11)] <- c("F_trans_ct", "F_trans_ct_lwr", "F_trans_ct_upr")
MT <- full_join(M_T, counterfactual_results[[2]], by = c("id", "T", "M_trans" = "original_outcome"))
colnames(MT)[c(9:11)] <- c("M_trans_ct", "M_trans_ct_lwr", "M_trans_ct_upr")
RT <- full_join(R_T, counterfactual_results[[3]], by = c("id", "T", "R_trans" = "original_outcome"))
colnames(RT)[c(9:11)] <- c("R_trans_ct", "R_trans_ct_lwr", "R_trans_ct_upr")
FR <- full_join(F_R, counterfactual_results[[4]], by = c("id", "T", "F_redem" = "original_outcome"))
colnames(FR)[c(9:11)] <- c("F_redem_ct", "F_redem_ct_lwr", "F_redem_ct_upr")
MR <- full_join(M_R, counterfactual_results[[5]], by = c("id", "T", "M_redem" = "original_outcome"))
colnames(MR)[c(9:11)] <- c("M_redem_ct", "M_redem_ct_lwr", "M_redem_ct_upr")
RR <- full_join(R_R, counterfactual_results[[6]], by = c("id", "T", "R_redem" = "original_outcome"))
colnames(RR)[c(9:11)] <- c("R_redem_ct", "R_redem_ct_lwr", "R_redem_ct_upr")
step1 <- left_join(FT, MT, by = c("id", "T", "D", "is_treat", "T0", "user", "proj"))
step2 <- left_join(step1, RT, by = c("id", "T", "D", "is_treat", "T0", "user", "proj"))
step3 <- left_join(step2, FR, by = c("id", "T", "D", "is_treat", "T0", "user", "proj"))
step4 <- left_join(step3, MR, by = c("id", "T", "D", "is_treat", "T0", "user", "proj"))
FINAL <- left_join(step4, RR, by = c("id", "T", "D", "is_treat", "T0", "user", "proj"))
Counterfactual Trend
Analysis
prefixes <- c("F_trans", "M_trans", "R_trans", "F_redem", "M_redem", "R_redem")
trend_plot <- function(prefix) {
plot <-
ggplot(FINAL, aes(x = T)) +
geom_ribbon(data = FINAL |> filter(is_treat == 1),
aes_string(ymin = paste0(prefix, "_ct_lwr"), ymax = paste0(prefix, "_ct_upr"), group = "id"),
alpha = 0.2, fill = "pink") +
geom_line(data = FINAL |> filter(is_treat == 0),
aes_string(y = prefix, group = "id"), color = "grey") +
geom_line(data = FINAL |> filter(is_treat == 1),
aes_string(y = paste0(prefix, "_ct"), group = "id", color = "factor(D)")) +
scale_color_manual(values = c("0" = "orange", "1" = "red")) +
scale_y_continuous(name = paste0(prefix)) +
theme_minimal() +
labs(title = paste0("Plot for ", prefix))
return(plot)
}
lapply(prefixes, trend_plot)
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
Average Treatment
Effect on the Treated (ATT)
unique(FINAL$T0)
## [1] 53 20 30 25 35
FINAL_SA20 <- filter(FINAL, T0 == 20)
SA20_id <- unique(FINAL_SA20$id)
FINAL_SA25 <- filter(FINAL, T0 == 25)
SA25_id <- unique(FINAL_SA25$id)
FINAL_SA30 <- filter(FINAL, T0 == 30)
SA30_id <- unique(FINAL_SA30$id)
FINAL_SA35 <- filter(FINAL, T0 == 35)
SA35_id <- unique(FINAL_SA35$id)
titles <- c("Transaction (Frequency)", "Transaction (Monetary)", "Transaction (Recency)",
"Redemption (Frequency)", "Redemption (Monetary)", "Redemption (Recency)")
Staggered Adoption
(T0 = 20)
OUT_ATT_SA20 <- lapply(1:6, function(i) {
effSummary(OUT[[i]], usr.id = SA20_id, cumu = FALSE, rela.period = TRUE)
})
ATT_view_SA20 <- lapply(1:6, function(i) {
ggplot(OUT_ATT_SA20[[i]]$est.eff, aes(x = time)) +
geom_hline(yintercept = 0, linetype = "dashed", color = "darkgrey") +
geom_vline(xintercept = 0, linetype = "dashed", color = "steelblue") +
geom_ribbon(aes(ymin = estimated_ATT_ci_l, ymax = estimated_ATT_ci_u),
alpha = 0.3, fill = "pink") +
geom_line(aes(y = estimated_ATT), color = "black") +
ggtitle(titles[i]) + ylab("ATT") + xlab("Relative Time") +
theme_minimal()
})
grid.arrange(grobs = ATT_view_SA20, ncol = 3, nrow = 2)
Staggered Adoption
(T0 = 25)
OUT_ATT_SA25 <- lapply(1:6, function(i) {
effSummary(OUT[[i]], usr.id = SA25_id, cumu = FALSE, rela.period = TRUE)
})
ATT_view_SA25 <- lapply(1:6, function(i) {
ggplot(OUT_ATT_SA25[[i]]$est.eff, aes(x = time)) +
geom_hline(yintercept = 0, linetype = "dashed", color = "darkgrey") +
geom_vline(xintercept = 0, linetype = "dashed", color = "steelblue") +
geom_ribbon(aes(ymin = estimated_ATT_ci_l, ymax = estimated_ATT_ci_u),
alpha = 0.3, fill = "pink") +
geom_line(aes(y = estimated_ATT), color = "black") +
ggtitle(titles[i]) + ylab("ATT") + xlab("Relative Time") +
theme_minimal()
})
grid.arrange(grobs = ATT_view_SA25, ncol = 3, nrow = 2)
Staggered Adoption
(T0 = 30)
OUT_ATT_SA30 <- lapply(1:6, function(i) {
effSummary(OUT[[i]], usr.id = SA30_id, cumu = FALSE, rela.period = TRUE)
})
ATT_view_SA30 <- lapply(1:6, function(i) {
ggplot(OUT_ATT_SA30[[i]]$est.eff, aes(x = time)) +
geom_hline(yintercept = 0, linetype = "dashed", color = "darkgrey") +
geom_vline(xintercept = 0, linetype = "dashed", color = "steelblue") +
geom_ribbon(aes(ymin = estimated_ATT_ci_l, ymax = estimated_ATT_ci_u),
alpha = 0.3, fill = "pink") +
geom_line(aes(y = estimated_ATT), color = "black") +
ggtitle(titles[i]) + ylab("ATT") + xlab("Relative Time") +
theme_minimal()
})
grid.arrange(grobs = ATT_view_SA30, ncol = 3, nrow = 2)
Staggered Adoption
(T0 = 35)
OUT_ATT_SA35 <- lapply(1:6, function(i) {
effSummary(OUT[[i]], usr.id = SA35_id, cumu = FALSE, rela.period = TRUE)
})
ATT_view_SA35 <- lapply(1:6, function(i) {
ggplot(OUT_ATT_SA35[[i]]$est.eff, aes(x = time)) +
geom_hline(yintercept = 0, linetype = "dashed", color = "darkgrey") +
geom_vline(xintercept = 0, linetype = "dashed", color = "steelblue") +
geom_ribbon(aes(ymin = estimated_ATT_ci_l, ymax = estimated_ATT_ci_u),
alpha = 0.3, fill = "pink") +
geom_line(aes(y = estimated_ATT), color = "black") +
ggtitle(titles[i]) + ylab("ATT") + xlab("Relative Time") +
theme_minimal()
})
grid.arrange(grobs = ATT_view_SA35, ncol = 3, nrow = 2)
ATT Without Staggered
Adoption
OUT_ATT <- lapply(1:6, function(i) {
effSummary(OUT[[i]], usr.id = NULL, cumu = FALSE, rela.period = TRUE)
})
ATT_view <- lapply(1:6, function(i) {
ggplot(OUT_ATT[[i]]$est.eff, aes(x = time)) +
geom_hline(yintercept = 0, linetype = "dashed", color = "darkgrey") +
geom_vline(xintercept = 0, linetype = "dashed", color = "steelblue") +
geom_ribbon(aes(ymin = estimated_ATT_ci_l, ymax = estimated_ATT_ci_u),
alpha = 0.3, fill = "pink") +
geom_line(aes(y = estimated_ATT), color = "black") +
ggtitle(titles[i]) + ylab("ATT") + xlab("Relative Time") +
theme_minimal()
})
grid.arrange(grobs = ATT_view, ncol = 3, nrow = 2)
