<jupyter_start><jupyter_text>Multi-investment Attribution: Distinguish the Effects of Multiple Outreach EffortsA startup that sells software would like to know whether its multiple outreach efforts were successful in attracting new customers or boosting consumption among existing customers. They would also like to distinguish the effects of several incentives on different kinds of customers. In other words, they would like to learn the **heterogeneous treatment effect** of each investment on customers' software usage. In an ideal world, the startup would run several randomized experiments where each customer would receive a random assortment of investments. However, this can be logistically prohibitive or strategically unsound: the startup might not have the resources to design such experiments or they might not want to risk losing out on big opportunities due to lack of incentives.In this customer scenario walkthrough, we show how tools from the [EconML](https://aka.ms/econml) library can use historical investment data to learn the effects of multiple investments.  Summary1. [Background](Background)2. [Data](Data)3. [Get Causal Effects with EconML](Get-Causal-Effects-with-EconML)4. [Understand Treatment Effects with EconML](Understand-Treatment-Effects-with-EconML)5. [Make Policy Decisions with EconML](Make-Policy-Decisions-with-EconML)6. [Conclusions](Conclusions)  BackgroundIn this scenario, a startup that sells software provides two types of incentives to its customers: technical support and discounts. A customer might be given one, both or none of these incentives. The startup has historical data on these two investments for 2,000 customers, as well as how much revenue these customers generated in the year after the investments were made. They would like to use this data to learn the optimal incentive policy for each existing or new customer in order to maximize the return on investment (ROI).The startup faces two challenges: 1) the dataset is biased because historically the larger customers received the most incentives and 2) the observed outcome combines effects from two different investments. Thus, they need a causal model that can accommodate multiple concurrent interventions. **Solution:** EconML’s `Doubly Robust Learner` model jointly estimates the effects of multiple discrete treatments. The model uses flexible functions of observed customer features to filter out spurious correlations in existing data and deliver the causal effect of each intervention on revenue.<jupyter_code># Some imports to get us started
import warnings
warnings.simplefilter('ignore')
# Utilities
import os
import urllib.request
import numpy as np
import pandas as pd

# Generic ML imports
from xgboost import XGBRegressor, XGBClassifier

# EconML imports
from econml.dr import LinearDRLearner

import matplotlib.pyplot as plt
import seaborn as sns

%matplotlib inline<jupyter_output><empty_output><jupyter_text>DataThe data* contains ~2,000 customers and is comprised of:* Customer features: details about the industry, size, revenue, and technology profile of each customer.* Interventions: information about which incentive was given to a customer.* Outcome: the amount of product the customer bought in the year after the incentives were given.Feature Name | Type | Details :--- |:--- |:--- **Global Flag** | W | whether the customer has global offices**Major Flag** | W | whether the customer is a large consumer in their industry (as opposed to SMC - Small Medium Corporation - or SMB - Small Medium Business)**SMC Flag** | W | whether the customer is a Small Medium Corporation (SMC, as opposed to major and SMB)**Commercial Flag** | W | whether the customer's business is commercial (as opposed to public secor)**IT Spend** | W | \\$ spent on IT-related purchases **Employee Count** | W | number of employees**PC Count** | W | number of PCs used by the customer**Size** | X | customer's size given by their yearly total revenue **Tech Support** | T | whether the customer received tech support (binary)**Discount** | T | whether the customer was given a discount (binary)**Revenue** | Y | \\$ Revenue from customer given by the amount of software purchased**To protect the privacy of the startup's customers, the data used in this scenario is synthetically generated and the feature distributions don't correspond to real distributions. However, the feature names have preserved their names and meaning.*<jupyter_code># Import the sample multi-attribution data
file_url = "https://msalicedatapublic.blob.core.windows.net/datasets/ROI/multi_attribution_sample.csv"
multi_data = pd.read_csv(file_url)
# Data sample
multi_data.head()
# Define estimator inputs
T_bin = multi_data[
    ["Tech Support", "Discount"]
]  # multiple interventions, or treatments
Y = multi_data["Revenue"]  # amount of product purchased, or outcome
X = multi_data[["Size"]]  # heterogeneity feature
W = multi_data.drop(
    columns=["Tech Support", "Discount", "Revenue", "Size"]
)  # controls<jupyter_output><empty_output><jupyter_text>We investigate below whether the number of investments given is correlated with the size of the customer. We note that the average customer size is larger for more incentives given.<jupyter_code># Average customer size per incentive combination
multi_data[["Size", "Tech Support", "Discount"]].groupby(
    by=["Tech Support", "Discount"], as_index=False
).mean().astype(int)<jupyter_output><empty_output><jupyter_text>The data was generated using the following underlying treatment effect function:$$\text{treatment_effect(Size)} = (5,000 + 2\% \cdot \text{Size}) \cdot I_\text{Tech Support} + (5\% \cdot \text{Size}) \cdot I_\text{Discount}$$Therefore, the treatment effect depends on the customer's size as follows: tech support provides an consumption boost of \$5,000 + 2\% Size and a discount provides an consumption boost of 5\% Size.**This is the relationship we seek to learn from the data.**<jupyter_code># Define underlying treatment effect function
TE_fn = lambda X: np.hstack([5000 + 2 / 100 * X, 5 / 100 * X])
true_TE = TE_fn(X)

# Define true coefficients for the three treatments
# The third coefficient is just the sum of the first two since we assume an additive effect
true_coefs = [2 / 100, 5 / 100, 7 / 100]
true_intercepts = [5000, 0, 5000]
treatment_names = ["Tech Support", "Discount", "Tech Support & Discount"]<jupyter_output><empty_output><jupyter_text>Get Causal Effects with EconML To get causal effects, we use EconML's `LinearDRLearner`* estimator. This estimator requires a set of discrete treatments $T$ that corresponds to different types of interventions. Thus, we first map the binary interventions tech support and discount into one categorical variable:Tech support| Discount| Treatment encoding| Details:--- |:--- |:--- |:---0 | 0 | 0 | no incentive1 | 0 | 1 | tech support only0 | 1 | 2 | discount only1 | 1 | 3 | both incentivesThe estimator takes as input the outcome of interest $Y$ (amount of product purchased), a discrete treatment $T$ (interventions given), heterogeneity features $X$ (here, customer's size) and controls $W$ (all other customer features).The LinearDRLearner also requires two auxiliary models to model the relationships $T\sim (W, X)$ (`model_propensity`) and $Y \sim (W, X)$(`model_regression`). These can be generic, flexible classification and regression models, respectively.  **This estimator assumes a linear relationship between the treatment effect and a transformation of the features $X$ (e.g. a polynomial basis expansion). For more generic forms of the treatment effect, see the `DRLearner` estimator.*<jupyter_code># Transform T to one-dimensional array with consecutive integer encoding
def treat_map(t):
    return np.dot(t, 2 ** np.arange(t.shape[0]))


T = np.apply_along_axis(treat_map, 1, T_bin).astype(int)
# Train EconML model with generic helper models
model = LinearDRLearner(
    model_regression=XGBRegressor(learning_rate=0.1, max_depth=3),
    model_propensity=XGBClassifier(learning_rate=0.1, max_depth=3, objective="multi:softmax"),
    random_state=1,
)
# Specify final stage inference type and fit model
model.fit(Y=Y, T=T, X=X, W=W, inference="statsmodels")<jupyter_output>[06:14:07] WARNING: C:/Users/Administrator/workspace/xgboost-win64_release_1.3.0/src/learner.cc:1061: Starting in XGBoost 1.3.0, the default evaluation metric used with the objective 'multi:softprob' was changed from 'merror' to 'mlogloss'. Explicitly set eval_metric if you'd like to restore the old behavior.
[06:14:08] WARNING: C:/Users/Administrator/workspace/xgboost-win64_release_1.3.0/src/learner.cc:1061: Starting in XGBoost 1.3.0, the default evaluation metric used with the objective 'multi:softprob' was changed from 'merror' to 'mlogloss'. Explicitly set eval_metric if you'd like to restore the old behavior.<jupyter_text>Understand Treatment Effects with EconMLWe can obtain a summary of the coefficient values as well as confidence intervals by calling the `summary` function on the fitted model for each treatment.<jupyter_code>for i in range(model._d_t[0]):
    print(f"Investment: {treatment_names[i]}")
    print(f"True treatment effect: {true_intercepts[i]} + {true_coefs[i]}*Size")
    display(model.summary(T=i + 1))<jupyter_output>Investment: Tech Support
True treatment effect: 5000 + 0.02*Size<jupyter_text>From the summary panels, we see that the learned coefficients/intercepts are close to the true coefficients/intercepts and the p-values are small for most of these. We further use the `coef_, coef__interval` and the `intercept_, intercept__interval` methods to obtain the learned coefficient values and build confidence intervals. We compare the true and the learned coefficients through the plots below.<jupyter_code># Compare learned coefficients with true model coefficients
# Aggregate data
coef_indices = np.arange(model._d_t[0])
coefs = np.hstack([model.coef_(T=i) for i in 1 + coef_indices])
intercepts = np.hstack([model.intercept_(T=i) for i in 1 + coef_indices])

# Calculate coefficient error bars for 95% confidence interval
coef_error = np.hstack([model.coef__interval(T=i) for i in 1 + coef_indices])
coef_error[0, :] = coefs - coef_error[0, :]
coef_error[1, :] = coef_error[1, :] - coefs

# Calculate intercept error bars for 95% confidence interval
intercept_error = np.vstack(
    [model.intercept__interval(T=i) for i in 1 + coef_indices]
).T
intercept_error[0, :] = intercepts - intercept_error[0, :]
intercept_error[1, :] = intercept_error[1, :] - intercepts
# Plot coefficients
plt.figure(figsize=(6, 5))
ax1 = plt.subplot(2, 1, 1)
plt.errorbar(
    coef_indices,
    coefs,
    coef_error,
    fmt="o",
    label="Learned values\nand 95% confidence interval",
)
plt.scatter(coef_indices, true_coefs, color="C1", label="True values", zorder=3)
plt.xticks(coef_indices, treatment_names)
plt.setp(ax1.get_xticklabels(), visible=False)
plt.title("Coefficients")
plt.legend(loc=(1.05, 0.65))
plt.grid()

# Plot intercepts
plt.subplot(2, 1, 2)
plt.errorbar(coef_indices, intercepts, intercept_error, fmt="o")
plt.scatter(coef_indices, true_intercepts, color="C1", zorder=3)
plt.xticks(coef_indices, treatment_names)
plt.title("Intercepts")
plt.grid()
plt.show()<jupyter_output><empty_output><jupyter_text>Make Policy Decisions with EconMLInvestments such as tech support and discounts come with an associated cost. Thus, we would like to know what incentives to give to each customer to maximize the profit from their increased engagement. This is the **treatment policy**.In this scenario, we define a cost function as follows:* The cost of `tech support` scales with the number of PCs a customer has. You can imagine that if the software product needs tech support to be installed on each machine, there is a cost (\\$100 here) per machine.* The cost of `discount` is a fixed \\$7,000. Think of this as giving the customer the first \\$7,000 worth of product for free.* The cost of `tech support` and `discount` is the sum of the cost of each of these. Note that this might not be the case in every business application: it is possible that managing multiple incentive programs can add overhead.<jupyter_code># Define cost function
def cost_fn(multi_data):
    t1_cost = multi_data[["PC Count"]].values * 100
    t2_cost = np.ones((multi_data.shape[0], 1)) * 7000
    return np.hstack([t1_cost, t2_cost, t1_cost + t2_cost])<jupyter_output><empty_output><jupyter_text>We use the model's `const_marginal_effect` method to find the counterfactual treatment effect for each possible treatment. We then subtract the treatment cost and choose the treatment which the highest return. That is the recommended policy.<jupyter_code># Get roi for each customer and possible treatment
potential_roi = model.const_marginal_effect(X=X.values) - cost_fn(multi_data)
# Add a column of 0s for no treatment
potential_roi = np.hstack([np.zeros((X.shape[0], 1)), potential_roi])
all_treatments = np.array(["None"] + treatment_names)
recommended_T = np.argmax(potential_roi, axis=1)
ax1 = sns.scatterplot(
    x=X.iloc[:, 0].values.flatten(),
    y=multi_data["PC Count"].values,
    hue=all_treatments[recommended_T],
    hue_order=all_treatments,
    cmap="Dark2",
    s=40,
)
plt.legend(title="Investment Policy")
plt.setp(
    ax1,
    xlabel="Customer Size",
    ylabel="PC Count",
    title="Optimal Investment Policy by Customer",
)
plt.show()<jupyter_output><empty_output><jupyter_text>We compare different policies: the optimal policy we learned, the current policy, and the policy under which each customer is given all incentives. We note that the optimal policy has a much higher ROI than the alternatives.<jupyter_code>roi_current = potential_roi[np.arange(X.shape[0]), T].sum()
roi_optimal = potential_roi[np.arange(X.shape[0]), recommended_T].sum()
roi_bothT = potential_roi[:, -1].sum()
all_rois = np.array([roi_optimal, roi_current, roi_bothT])
Y_baseline = (Y - model.effect(X=X.values, T1=T)).sum()
pd.DataFrame(
    {
        "Policy": ["Optimal", "Current", "All Investments"],
        "ROI ($)": all_rois,
        "ROI (% of baseline Y)": np.round(all_rois / Y_baseline * 100, 1),
    }
)<jupyter_output><empty_output><jupyter_text>Performance of policies based on ground truth<jupyter_code>true_roi = np.zeros((X.shape[0], 4))
true_roi[:, 1:3] = TE_fn(X.iloc[:, [0]])
true_roi[:, 3] = np.sum(true_roi[:, 1:3], axis=1)
true_roi[:, 1:] -= cost_fn(multi_data)
roi_current = true_roi[np.arange(X.shape[0]), T].sum()
roi_optimal = true_roi[np.arange(X.shape[0]), recommended_T].sum()
roi_bothT = true_roi[:, -1].sum()
all_rois = np.array([roi_optimal, roi_current, roi_bothT])
Y_baseline = (Y - model.effect(X=X.values, T1=T)).sum()
pd.DataFrame(
    {
        "Policy": ["Optimal", "Current", "All Investments"],
        "ROI ($)": all_rois,
        "ROI (% of baseline Y)": np.round(all_rois / Y_baseline * 100, 1),
    }
)<jupyter_output><empty_output><jupyter_text>Policy Decisions Using Doubly Robust Policy Learning<jupyter_code>from econml.policy import DRPolicyForest

X = multi_data[['Size', 'PC Count']]
est = DRPolicyForest(random_state=1,
                     min_propensity=1e-3,
                     min_impurity_decrease=0.1,
                     min_samples_leaf=40,
                     max_samples=.6,
                     honest=True,
                     max_depth=4)
costs = np.hstack([np.zeros((X.shape[0], 1)), cost_fn(multi_data)])
est.fit(Y - costs[np.arange(X.shape[0]), T], T, X=X)
recommended_T = est.predict(X.values)
ax1 = sns.scatterplot(
    x=X.iloc[:, 0].values.flatten(),
    y=multi_data["PC Count"].values,
    hue=all_treatments[recommended_T],
    hue_order=all_treatments,
    cmap="Dark2",
    s=40,
)
plt.legend(title="Investment Policy")
plt.setp(
    ax1,
    xlabel="Customer Size",
    ylabel="PC Count",
    title="Optimal Investment Policy by Customer",
)
plt.show()
plt.figure(figsize=(15, 5))
est.plot(0, max_depth=2, treatment_names=['None'] + treatment_names)<jupyter_output><empty_output><jupyter_text>Performance based on ground truth<jupyter_code>roi_current = true_roi[np.arange(X.shape[0]), T].sum()
roi_optimal = true_roi[np.arange(X.shape[0]), recommended_T].sum()
roi_bothT = true_roi[:, -1].sum()
all_rois = np.array([roi_optimal, roi_current, roi_bothT])
Y_baseline = (Y - model.effect(X=X.values[:, [0]], T1=T)).sum()
pd.DataFrame(
    {
        "Policy": ["Optimal", "Current", "All Investments"],
        "ROI ($)": all_rois,
        "ROI (% of baseline Y)": np.round(all_rois / Y_baseline * 100, 1),
    }
)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1031390134529148, 'fraction_numerical': 0.015605381165919283, 'mean_word_length': 4.8708318708318705, 'pattern_counts': {'":': 9, '<': 51, '<?xml version=': 0, '>': 51, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '57793', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 5312, 'n_tokens_neox': 4970, 'n_words': 1783, 'path': 'notebooks/CustomerScenarios/Case Study - Multi-investment Attribution at A Software Company.ipynb', 'repo_name': 'GildasTiwangNg/EconML', 'size': '16725'}