Japan Sotyktu CE3 Lift Analysis
Kyle Rechard
2024-08-08
Causal Inference with a Matching and Difference-in-Difference Regression
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
This R markdown document implements a four(4) step procedure for post-hoc test/control difference-in-difference analysis using “coarsened exact matching”. The code implements a simplified version of the method outlined in https://gking.harvard.edu/publications/cem-software-coarsened-exact-matchings.
The R package “MatchIt”, https://cran.r-project.org/web//packages//MatchIt/vignettes/MatchIt.html, also has an option to implement the full version of coarsened exact match, which outputs observation weights on the control arm for use in a regression. This is less transparent, however than the sampling based on the multinomial distribution.
Japan Data Load
library(magrittr)# Data Load/cleaning
Data_initial <- read.csv('C:/Users/rechardk/OneDrive - Bristol Myers Squibb/Documents/CE3/Sotyktu_cluster_data_202407_2.csv') %>%
dplyr::select(clst_cd,month, sales, hcp_veeva, calls, ce3_calls, suggested) %>%
dplyr::mutate(Account.id = clst_cd)
BIA_Sotyktu <- read.csv('C:/Users/rechardk/OneDrive - Bristol Myers Squibb/Documents/CE3/BMS-ZS_Encluster_Sotyktu_Otezla_Sales_v1.0.csv') %>%
dplyr::mutate(acc = clst_cd) %>%
dplyr::mutate(date = as.Date(paste0(cald_dt_mth_id,"01"),tryFormats = c("%Y/%m/%d", "%Y%m%d"))) %>%
dplyr::group_by(date) %>%
dplyr::summarise(sales = sum(sales_new, na.rm = TRUE)) %>%
dplyr::ungroup()
head(Data_initial %>% dplyr::select(-Account.id))## clst_cd month sales hcp_veeva calls ce3_calls suggested
## 1 1.02011e+11 2023-01-01 83127 0016F00003N4eMlQAJ 0 0 0
## 2 1.02011e+11 2023-01-01 83127 0016F00003N5CnXQAV 0 0 0
## 3 1.02011e+11 2023-01-01 83127 0016F00003N5x3OQAR 0 0 0
## 4 1.02011e+11 2023-01-01 83127 0016F00003N5BP5QAN 0 0 0
## 5 1.02011e+11 2023-01-01 83127 0016F00003LxeO1QAJ 0 0 0
## 6 1.02011e+11 2023-01-01 83127 0016F00003N5dzdQAB 0 0 0CE3 launched in Japan for Sotyktu in November 2023. We consider the 7 months before November 2023 to be the pre-CE3 period, which we cluster level sales to the 7 months after the launch of CE3 (“post-CE3). Each cluster (account) will thus have two rows (pre- and post-CE3), identified by the dummy variable”Time”, with “0” for pre-CE3 and “1” for post-CE3.
The delineation between the “test” and “control” accounts will be based on the percentage of calls that are categorized as executed CE3 calls in the post-CE3 period, so “ce3_percentage” variable is created as well.
This cleaning is shown below:
Data <- Data_initial %>%
dplyr::distinct(Account.id, month, sales, .keep_all = TRUE) %>%
dplyr::filter(month>="2023-04-01") %>%
dplyr::filter(month<="2024-05-01") %>%
dplyr::mutate(Time = ifelse(month<"2023-11-01", 0, 1)) %>%
dplyr::group_by(Account.id, Time) %>%
dplyr::summarise(Sales = sum(sales, na.rm = TRUE),
calls = sum(calls, na.rm = TRUE),
ce3_calls = sum(ce3_calls, na.rm = TRUE),
suggested = sum(suggested, na.rm =TRUE)) %>%
dplyr::mutate(ce3_percentage = ifelse(calls>0,ce3_calls/calls, 0)) %>%
dplyr::ungroup()## `summarise()` has grouped output by 'Account.id'. You can override using the
## `.groups` argument.head(Data)## # A tibble: 6 × 7
## Account.id Time Sales calls ce3_calls suggested ce3_percentage
## <dbl> <dbl> <dbl> <int> <int> <int> <dbl>
## 1 102011001791 0 831270 2 0 0 0
## 2 102011001791 1 1828794 2 1 4 0.5
## 3 102011001793 1 83127 0 0 0 0
## 4 102012000021 0 47501. 3 0 0 0
## 5 102012000021 1 35626. 2 0 0 0
## 6 102021001774 1 665016 16 8 18 0.5Step 1: Outliering
Percentile thresholds for eliminating accounts larger than the XXth percentile (I chose the 95th percentile, but the 99th would also be a good default choice) in both volume and sales growth.
Outlier Parameters
#Parameters
pre_group_outlier_percentile = 0.95
post_group_outlier_percentile = 0.95
sales_growth_percentile_cutoff = 0.95Outliers_pre <- Data %>%
dplyr::filter(Time==0) %>%
dplyr::mutate(Sales_quantile_control = quantile(Sales, probs = pre_group_outlier_percentile, na.rm =TRUE)) %>%
dplyr::filter(Sales>Sales_quantile_control) %>%
dplyr::select(Account.id)
Outliers_post <- Data %>%
dplyr::filter(Time==1) %>%
dplyr::mutate(Sales_quantile_test = quantile(Sales, probs = post_group_outlier_percentile, na.rm =TRUE)) %>%
dplyr::filter(Sales>Sales_quantile_test) %>%
dplyr::select(Account.id)
Data1a <- Data %>%
dplyr::filter(!(Account.id %in% c(Outliers_pre$Account.id))) %>%
dplyr::filter(!(Account.id %in% c(Outliers_post$Account.id)))
####
# Remove Sales Growth Outliers
sales_growth <- Data1a %>%
dplyr::group_by(Account.id) %>%
dplyr::mutate(sales_total = sum(Sales, na.rm = TRUE),
sales_growth = ifelse(Sales>0,abs(((sales_total - Sales)/Sales-1)*100), 0)) %>%
dplyr::filter(Time == 0) %>%
dplyr::ungroup() %>%
dplyr::mutate(sales_growth_percentile = quantile(sales_growth, probs = sales_growth_percentile_cutoff, na.rm = TRUE)
) %>%
dplyr::select(Account.id, sales_growth,sales_growth_percentile) %>%
dplyr::ungroup()
Data1 <- Data1a %>%
dplyr::left_join(sales_growth, by = "Account.id") %>%
dplyr::filter(sales_growth<sales_growth_percentile)
head(Data1)## # A tibble: 6 × 9
## Account.id Time Sales calls ce3_calls suggested ce3_percentage sales_growth
## <dbl> <dbl> <dbl> <int> <int> <int> <dbl> <dbl>
## 1 1.02e11 0 8.31e5 2 0 0 0 120
## 2 1.02e11 1 1.83e6 2 1 4 0.5 120
## 3 1.02e11 0 4.75e4 3 0 0 0 25
## 4 1.02e11 1 3.56e4 2 0 0 0 25
## 5 1.02e11 0 1.58e6 18 0 0 0 121.
## 6 1.02e11 1 3.49e6 24 16 30 0.667 121.
## # ℹ 1 more variable: sales_growth_percentile <dbl>Step 1(a)
# Account Sales Tiering
HML_designation = Data1 %>%
dplyr::filter(Time ==0) %>%
dplyr::mutate(
sales_decile = dplyr::ntile(Sales,10),
activity_decile = dplyr::ntile(calls,10),
sales_tier = dplyr::case_when(
sales_decile > 7 ~ "High",
sales_decile <= 7 & sales_decile > 4 ~ "Medium",
sales_decile <= 4 ~ "Low",
TRUE ~ "other"
),
F2F_tier = dplyr::case_when(
activity_decile > 7 ~ "High",
activity_decile <= 7 & activity_decile > 4 ~ "Medium",
activity_decile <= 4 ~ "Low",
TRUE ~ "other"
)
) %>%
dplyr::select(Account.id, sales_decile,activity_decile,sales_tier, F2F_tier)
##################################################
##################################################
# Join Data
Data2 = Data1 %>%
#dplyr::select(-ce3_percentage) %>%
dplyr::left_join(HML_designation,by = 'Account.id') %>% # Drop Sales tier zero
dplyr::filter(sales_tier!="other")
table(Data2$sales_tier,Data2$F2F_tier)##
## High Low Medium
## High 96 28 68
## Low 30 135 58
## Medium 65 64 64Step 2: Create Test/Control Groups
Test and control accounts are delineated based on the composition of CE3 calls in the post-CE3 period. We examine the distribution and to determine appropriate. It is “ok” to iterate different choices for the cutoffs here in order to create a somewhat bimodal distribution, but DO NOT choose these cutoffs to minimize the p-value on the lift coefficient in the diff-in-diff regression later on.
Based on the histogram, we choose the following cutoffs:
ce3_percentile_cutoff_lower = 0.22
ce3_percentile_cutoff_upper = 0.75
quantile(Data2 %>% dplyr::select(ce3_percentage), probs = ce3_percentile_cutoff_lower, na.rm = TRUE)## 22%
## 0quantile(Data2 %>% dplyr::select(ce3_percentage), probs = ce3_percentile_cutoff_upper, na.rm = TRUE)## 75%
## 0.35All accounts below the 22nd percentile of the CE3 composition distribution (0% CE3 call composition in the post-CE3 period) were considered the “control” group. The 25% of accounts greater than the 75th percentile were considered the “Test” group.
Test_control <- Data2 %>%
dplyr::filter(Time==1) %>%
dplyr::mutate(ce3_percentage_percentile_lower = quantile(ce3_percentage, probs = ce3_percentile_cutoff_lower, na.rm = TRUE),
ce3_percentage_percentile_upper = quantile(ce3_percentage, probs = ce3_percentile_cutoff_upper, na.rm = TRUE)) %>%
dplyr::mutate(Test = ifelse(ce3_percentage <= ce3_percentage_percentile_lower,0,1),
Test = ifelse(ce3_percentage > ce3_percentage_percentile_lower & ce3_percentage <= ce3_percentage_percentile_upper,NA,Test)
) %>%
dplyr::select(Account.id,Test, ce3_percentage)
Test_control## # A tibble: 281 × 3
## Account.id Test ce3_percentage
## <dbl> <dbl> <dbl>
## 1 102011001791 1 0.5
## 2 102012000021 0 0
## 3 102031001821 1 0.667
## 4 102051001807 NA 0.432
## 5 102091001817 1 0.545
## 6 102101001779 NA 0.357
## 7 111051001550 NA 0.312
## 8 111061001592 NA 0.4
## 9 112011001655 1 0.472
## 10 112031001563 0 0
## # ℹ 271 more rows#join Test/Control
Data2 = Data2 %>%
dplyr::select(-ce3_percentage) %>%
dplyr::left_join(Test_control,by="Account.id")Step 3: Matching
We do “Coarsened Exact Matching” algorithm in the R MatchIt package as a way to remove bias from unobserved confounders in of Diff-in-Diff. It is necessary to always match the outcome variable (sales) in the test and control group, but also a parsimonious list of other variables that are the most likely to can affect the outcome (F2F interactions are generally the most influential marketing channel, so they are included).
Because the quality of the matches will degrade the more variables the algorithm is forced to match on, it is best to match on the outcome and at most one or two other influential channels in the pre-period. Total F2F calls are generally the best choice, besides our outcome variable Sales, to use to find matching Test and Control accounts in the pre-CE3 period. Other possible confounding market channels (RTE, Emails, etc.) can be included as additional control variables in the diff-in-diff regression.
First, we examine the balance of sales and F2F calls in the pre-CE3 period (without doing any matching):
#Matchit
# Match on Pre-CE3 Sales and total F2F calls
Data_matchit <- Data2 %>%
dplyr::filter(!is.na(Test)) %>%
dplyr::filter(Time==0)
# Pre matching
m.out0 <- MatchIt::matchit(Test ~ Sales +F2F_tier, data = Data_matchit,
method = NULL, distance = "glm")
summary(m.out0)##
## Call:
## MatchIt::matchit(formula = Test ~ Sales + F2F_tier, data = Data_matchit,
## method = NULL, distance = "glm")
##
## Summary of Balance for All Data:
## Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## distance 0.5491 0.5009 0.4135 1.3819 0.1246
## Sales 412178.4436 298196.9066 0.2878 1.3236 0.0876
## F2F_tierHigh 0.2429 0.2540 -0.0259 . 0.0111
## F2F_tierLow 0.4143 0.5397 -0.2546 . 0.1254
## F2F_tierMedium 0.3429 0.2063 0.2876 . 0.1365
## eCDF Max
## distance 0.1984
## Sales 0.1921
## F2F_tierHigh 0.0111
## F2F_tierLow 0.1254
## F2F_tierMedium 0.1365
##
## Sample Sizes:
## Control Treated
## All 63 70
## Matched 63 70
## Unmatched 0 0
## Discarded 0 0We want the variance ratio of the test and control groups for sales to be as close to 1 as possible. As you can see, the Test and Control groups were actually pretty well balanced prior to doing any matching (this is good luck, but will not always happen).
Now, using the R MatchIt package, we do “Coarsened Exact Matching” - after splitting the program weights the control observations in order to match the first two moments of the test group distribution of (in this case) in the pre-period.
m.out1 <- MatchIt::matchit(Test ~ Sales +F2F_tier, data = Data_matchit, method = "cem", replace =FALSE)
m.out1## A matchit object
## - method: Coarsened exact matching
## - number of obs.: 133 (original), 113 (matched)
## - target estimand: ATT
## - covariates: Sales, F2F_tiersummary(m.out1)##
## Call:
## MatchIt::matchit(formula = Test ~ Sales + F2F_tier, data = Data_matchit,
## method = "cem", replace = FALSE)
##
## Summary of Balance for All Data:
## Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## Sales 412178.4436 298196.9066 0.2878 1.3236 0.0876
## F2F_tierHigh 0.2429 0.2540 -0.0259 . 0.0111
## F2F_tierLow 0.4143 0.5397 -0.2546 . 0.1254
## F2F_tierMedium 0.3429 0.2063 0.2876 . 0.1365
## eCDF Max
## Sales 0.1921
## F2F_tierHigh 0.0111
## F2F_tierLow 0.1254
## F2F_tierMedium 0.1365
##
## Summary of Balance for Matched Data:
## Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## Sales 300554.1062 283024.2946 0.0443 0.8988 0.0312
## F2F_tierHigh 0.2632 0.2632 0.0000 . 0.0000
## F2F_tierLow 0.4561 0.4561 -0.0000 . 0.0000
## F2F_tierMedium 0.2807 0.2807 -0.0000 . 0.0000
## eCDF Max Std. Pair Dist.
## Sales 0.1178 0.1413
## F2F_tierHigh 0.0000 0.0000
## F2F_tierLow 0.0000 0.0000
## F2F_tierMedium 0.0000 0.0000
##
## Sample Sizes:
## Control Treated
## All 63. 70
## Matched (ESS) 40.28 57
## Matched 56. 57
## Unmatched 7. 13
## Discarded 0. 0You can see the distributions of
plot(m.out1, type = "density", interactive = FALSE,
which.xs = ~Sales +F2F_tier)
m.data <- MatchIt::match.data(m.out1) %>%
dplyr::select(Account.id, weights)
head(m.data)## # A tibble: 6 × 2
## Account.id weights
## <dbl> <dbl>
## 1 102011001791 1
## 2 102012000021 0.632
## 3 102091001817 1
## 4 112011001655 1
## 5 112031001563 1.11
## 6 112211001596 1.31############################################
# Final Analysis Data
Data_matched = Data2 %>%
dplyr::left_join(m.data,by= 'Account.id')
#Data3 = test_control %>%
# dplyr::left_join(Data2,by= 'Account.id')Step 4: Lift Analysis with Difference-in-Difference Regression
Perform a weighted least squares diff-in-diff regression (the weights on the control observations come from the matching algorithm in the pre-CE3 period).
#regression on matched data
#summary(lm(Sales~test_1 + Time + test_1:Time ,data = Data3))
summary(lm(Sales~Test*Time ,data = Data_matched, weights=weights))##
## Call:
## lm(formula = Sales ~ Test * Time, data = Data_matched, weights = weights)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -957660 -300554 -134300 281335 2949309
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 283024 75930 3.727 0.000246 ***
## Test 17530 106909 0.164 0.869905
## Time 331562 107381 3.088 0.002274 **
## Test:Time 242418 151193 1.603 0.110275
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 568200 on 222 degrees of freedom
## (382 observations deleted due to missingness)
## Multiple R-squared: 0.159, Adjusted R-squared: 0.1477
## F-statistic: 13.99 on 3 and 222 DF, p-value: 2.174e-08# Average Treatment effect of the treated
ATT = lm(Sales~Test + Time + Test:Time ,data = Data_matched)$coefficients[4]# Manually compute DiD estimator (check)
total_1_1 = Data_matched %>% dplyr::filter(Test ==1 & Time ==1) %>% dplyr::select(Sales) %>% sum()
total_1_0 = Data_matched %>% dplyr::filter(Test ==1 & Time ==0) %>% dplyr::select(Sales) %>% sum()
total_0_1 = Data_matched %>% dplyr::filter(Test ==0 & Time ==1) %>% dplyr::select(Sales) %>% sum()
total_0_0 = Data_matched %>% dplyr::filter(Test ==0 & Time ==0) %>% dplyr::select(Sales) %>% sum()
n_11 = Data_matched %>% dplyr::filter(Test ==1 & Time ==1) %>% dplyr::select(Account.id) %>% dplyr::distinct()
n_10 = Data_matched %>% dplyr::filter(Test ==1 & Time ==0) %>% dplyr::select(Account.id) %>% dplyr::distinct()
n_01 = Data_matched %>% dplyr::filter(Test ==0 & Time ==1) %>% dplyr::select(Account.id) %>% dplyr::distinct()
n_00 = Data_matched %>% dplyr::filter(Test ==0 & Time ==0) %>% dplyr::select(Account.id) %>% dplyr::distinct()
ATT_manual = (total_0_0/dim(n_00)[1] - total_0_1/dim(n_01)[1]) - (total_1_0/dim(n_10)[1] - total_1_1/dim(n_11)[1])
ATT_manual## [1] 240687.5Compute Average Treatment effect of the Treated
#sales lift percentage of counterfactual test group
counterfactual_per_account_sales = (total_1_0/dim(n_10)[1] + (total_0_1/dim(n_01)[1] - total_0_0/dim(n_00)[1]))
# Test over counterfactual control percentage
as.numeric(ATT/counterfactual_per_account_sales)*100## [1] 31.52158Compute two other sales lift percentages - Treatment impact of the treated as percentage of sa sales in the accounts with mapped HCPs (“aligned sales”) and Treatment impact of the treated as as a percentage of total sales.
total_mapped =Data %>% dplyr::filter(Time==1) %>% dplyr::select(Sales) %>% sum()
total_sales =BIA_Sotyktu %>%
dplyr::filter(date>="2023-11-01" & date<="2024-05-01" ) %>%
dplyr::select(sales) %>%
sum()
sum(m.out1$treat) # initial treatment sample size ## [1] 70summary(m.out1)$nn[9] # effective treatment sample size## [1] 57#Lift as a percentage of sales in the 502 clusters
as.numeric(ATT*summary(m.out1)$nn[9]/total_mapped)*100## [1] 2.510341#Lift as a percentage of total sales
as.numeric(ATT*summary(m.out1)$nn[9]/total_sales)*100## [1] 1.059708