作者簡(jiǎn)介

Yiwen，攜程數(shù)據(jù)分析師，專(zhuān)注用戶(hù)增長(zhǎng)、因果推斷、數(shù)據(jù)科學(xué)等領(lǐng)域。

一、背景

如何科學(xué)地推斷某個(gè)產(chǎn)品策略對(duì)觀測(cè)指標(biāo)產(chǎn)生的效應(yīng)非常重要，這能夠幫助產(chǎn)品和運(yùn)營(yíng)更精準(zhǔn)地得到該策略的價(jià)值，從而進(jìn)行后續(xù)方向的迭代及調(diào)整。

在因果推斷框架下，效果評(píng)估的黃金準(zhǔn)則一定是“AB實(shí)驗(yàn)”，因?yàn)閷?shí)驗(yàn)的分流被認(rèn)為是完全隨機(jī)且均勻的，在此基礎(chǔ)上對(duì)比實(shí)驗(yàn)組與對(duì)照組的指標(biāo)差異就可以體現(xiàn)某個(gè)干預(yù)帶來(lái)的增量值。但是很多場(chǎng)景下，我們較難進(jìn)行嚴(yán)格的AB實(shí)驗(yàn)，例如對(duì)于酒店的定價(jià)；現(xiàn)金獎(jiǎng)勵(lì)的發(fā)放等等，不適宜向不同人群展現(xiàn)不同的內(nèi)容。對(duì)于這些問(wèn)題，我們會(huì)采取因果推斷的方法來(lái)進(jìn)行策略的效果評(píng)估。

本文主要介紹BSTS模型原理以及CausalImpact對(duì)模型的代碼實(shí)現(xiàn)，旨在面對(duì)一些具有特定周期性特點(diǎn)的數(shù)據(jù)時(shí)，更精準(zhǔn)科學(xué)地進(jìn)行因果效應(yīng)值的估計(jì)。下文將首先對(duì)模型原理進(jìn)行簡(jiǎn)要闡釋?zhuān)浑S后利用模擬數(shù)據(jù)展示代碼邏輯，最后在具體的業(yè)務(wù)場(chǎng)景中進(jìn)行實(shí)踐。

二、現(xiàn)有方法及潛在問(wèn)題

大部分運(yùn)營(yíng)和產(chǎn)品在評(píng)估一些全量上線的策略效果時(shí)，最常用的方式就是看上線前后的效果差異。但這種方法最大的問(wèn)題在于其假設(shè)前提：假設(shè)上線的功能是唯一影響效果的變量（即沒(méi)有任何其他干預(yù)和混淆變量），但這個(gè)假設(shè)現(xiàn)實(shí)中往往很難實(shí)現(xiàn)。

于是我們嘗試使用更多因果推斷的方法，例如PSM（傾向分匹配法），在所有非實(shí)驗(yàn)組的用戶(hù)群中，找到與實(shí)驗(yàn)組用戶(hù)的特征非常相似的一群人，將他們的指標(biāo)數(shù)據(jù)（例如下單率，訂單收益等等）與實(shí)驗(yàn)組的用戶(hù)進(jìn)行對(duì)比，從而體現(xiàn)出干預(yù)帶來(lái)的影響。但這個(gè)方法較為依賴(lài)選取的用戶(hù)特征與最后的匹配效果。

再例如SCM（合成控制方法），利用一些未受干預(yù)的地區(qū)合成一個(gè)“類(lèi)似的虛擬地區(qū)”來(lái)與“上線策略的地區(qū)”進(jìn)行整體的對(duì)比。但這也需要一個(gè)關(guān)鍵假設(shè)：可以找到長(zhǎng)期變化趨勢(shì)高度同步的地區(qū)來(lái)進(jìn)行對(duì)照，而這個(gè)條件往往也很難實(shí)現(xiàn)。

進(jìn)而在傳統(tǒng)SCM的基礎(chǔ)上，我們企圖通過(guò)類(lèi)似集成學(xué)習(xí)的方法，將多個(gè)未干預(yù)的對(duì)照組作為輸入值，再結(jié)合實(shí)驗(yàn)組自身長(zhǎng)期的時(shí)間序列波動(dòng)情況，擬合出一個(gè)未受干預(yù)的虛擬對(duì)照組，從而將“對(duì)照組與實(shí)驗(yàn)組高度同步”的強(qiáng)假設(shè)降為弱假設(shè)。本文介紹的BSTS模型就是用來(lái)刻畫(huà)某種“長(zhǎng)期的時(shí)間序列波動(dòng)”的數(shù)據(jù)模型，CausalImpact是用來(lái)針對(duì)這樣的數(shù)據(jù)進(jìn)行因果效應(yīng)值的估計(jì)。下文中我們將詳細(xì)介紹這兩個(gè)工具。

三、模型介紹

BSTS模型 （Bayesian Structured Time Series）稱(chēng)為“貝葉斯結(jié)構(gòu)化時(shí)間序列”，正如其名，它的主要特點(diǎn)體現(xiàn)在：

• 適用于有結(jié)構(gòu)特征的時(shí)間序列數(shù)據(jù)
• 利用貝葉斯的思想來(lái)進(jìn)行參數(shù)估計(jì)

結(jié)構(gòu)化的時(shí)間序列數(shù)據(jù)在日常生活中不少見(jiàn)，尤其像攜程這樣的OTA行業(yè)，平臺(tái)的訂單情況其實(shí)是有一定時(shí)間規(guī)律的，例如周末和節(jié)假日是訂單高峰期；周中是訂單平峰期等。另一方面，貝葉斯的思想是指在得到樣本數(shù)據(jù)之前，即對(duì)要估計(jì)的參數(shù)有一些“先天認(rèn)知”），隨后基于這樣的認(rèn)知，結(jié)合樣本數(shù)據(jù)再得到后驗(yàn)分布（如下方公式展示）

故BSTS模型主要即對(duì)結(jié)構(gòu)化時(shí)序數(shù)據(jù)進(jìn)行模型擬合及預(yù)測(cè)，在擬合的過(guò)程中使用到了貝葉斯的先驗(yàn)思想。其好處是能夠給出預(yù)測(cè)值的置信區(qū)間，使得預(yù)測(cè)結(jié)果更科學(xué)可信。下文將對(duì)這幾種思想逐一進(jìn)行介紹。

3.1 狀態(tài)空間模型

結(jié)構(gòu)化的時(shí)間序列數(shù)據(jù)是指某一觀測(cè)數(shù)據(jù)的背后其實(shí)隱藏著隨時(shí)間變化而變化的不同狀態(tài)，其中觀測(cè)值與狀態(tài)值之間有對(duì)應(yīng)關(guān)系；不同時(shí)刻的狀態(tài)之間也有轉(zhuǎn)換關(guān)系。我們一般用以下?tīng)顟B(tài)空間模型來(lái)刻畫(huà)這兩種映射邏輯：

(1) 稱(chēng)為觀測(cè)方程，反映觀測(cè)值與其背后隱藏狀態(tài)的關(guān)系；(2) 稱(chēng)為狀態(tài)方程，反映隨時(shí)間推移各個(gè)狀態(tài)之間的轉(zhuǎn)換。；都是不同變量之間的“關(guān)系映射矩陣”； 是獨(dú)立于其他變量且服從正態(tài)分布的噪聲。所謂數(shù)據(jù)的“結(jié)構(gòu)化”，主要包括：

• Linear Local Trend（局部趨勢(shì)）：一定時(shí)間內(nèi)的單調(diào)性（單調(diào)上升或下降）
• Seasonality（季節(jié)性因子）：固定長(zhǎng)度的變化，類(lèi)似于一年四季的溫度變化
• Cyclical（周期性）：類(lèi)似季節(jié)性但波動(dòng)時(shí)間不固定，波動(dòng)頻率也不固定的變化

圖3-1：觀測(cè)數(shù)據(jù)及其結(jié)構(gòu)化元素。第一張圖體現(xiàn)原數(shù)據(jù)的波動(dòng)情況；第二張?bào)w現(xiàn)季節(jié)性因子的情況；第三張圖體現(xiàn)局部趨勢(shì)的情況。

如果希望在映射關(guān)系中加入?yún)f(xié)變量X，可以將(1)拓展為：

其中表示協(xié)變量X與觀測(cè)數(shù)據(jù)之間的關(guān)系，如果協(xié)變量項(xiàng)表現(xiàn)很好（如有顯著影響）的話，那對(duì)應(yīng)的local trend就會(huì)相對(duì)較弱。上述三個(gè)方程中的參數(shù)將在后文中展示估計(jì)方式。

3.2 貝葉斯及MCMC（馬爾可夫蒙特卡洛方法）

假設(shè)狀態(tài)方程(2)中各個(gè)時(shí)刻的狀態(tài)序列為表示模型中所有的參數(shù)。我們現(xiàn)在希望對(duì)θ進(jìn)行估計(jì)，核心步驟如下：

• 對(duì)θ設(shè)置先驗(yàn)分布以及初始狀態(tài)的分布 
• 構(gòu)造馬爾科夫鏈，用MCMC方法得到
• 通過(guò)貝葉斯公式計(jì)算得到參數(shù)的后驗(yàn)分布

下面對(duì)于各個(gè)步驟中用到的方法進(jìn)行簡(jiǎn)要說(shuō)明：

1）貝葉斯估計(jì)：BSTS模型的一大特點(diǎn)就是在參數(shù)估計(jì)上使用了貝葉斯估計(jì)，即在估計(jì)之前先給出參數(shù)設(shè)置先驗(yàn)分布，隨后再結(jié)合樣本數(shù)據(jù)給出參數(shù)的后驗(yàn)分布。不同類(lèi)型的參數(shù)一般有一些常用的先驗(yàn)分布，例如均值一般使用正態(tài)分布，，方差使用inverse-Gamma分布，協(xié)方差矩陣可以使用IW分布等等。值得注意的是，先驗(yàn)分布的設(shè)置一定程度上會(huì)影響后續(xù)MCMC收斂的情況以及后驗(yàn)分布的準(zhǔn)確性，因此并不能太過(guò)隨意地設(shè)置先驗(yàn)分布，應(yīng)盡可能多地根據(jù)實(shí)際數(shù)據(jù)推導(dǎo)出最合適的先驗(yàn)分布，或是比較各先驗(yàn)分布下后驗(yàn)分布和似然函數(shù)的值來(lái)進(jìn)行選擇。

2）MCMC方法：我們嘗試構(gòu)造一條馬爾可夫鏈（一種特殊的序列，當(dāng)前時(shí)刻的狀態(tài)值僅與前一時(shí)刻的狀態(tài)值有關(guān)，最終序列會(huì)收斂到某個(gè)穩(wěn)定的分布），使得其最終收斂的穩(wěn)態(tài)分布就是參數(shù)的后驗(yàn)分布。這一過(guò)程我們可以通過(guò)Gibbs采樣實(shí)現(xiàn)：設(shè)置先驗(yàn)分布之后，從初始狀態(tài)出發(fā)，每次固定α采樣θ；再固定θ采樣α，逐漸一次次更新兩組參數(shù)，最終形成一條服從馬爾可夫性質(zhì)的鏈路，可以證明其穩(wěn)態(tài)收斂的分布就是，其中代表所有的觀測(cè)數(shù)據(jù)。

3）預(yù)測(cè)值估計(jì)：得到之后，我們從該分布中對(duì)（α，θ）進(jìn)行采樣，再代入狀態(tài)空間方程(1)中對(duì)y進(jìn)行預(yù)測(cè)，得到，其中表示時(shí)間點(diǎn)n之后y的預(yù)測(cè)值。       

圖3-2：展示某結(jié)構(gòu)化時(shí)序數(shù)據(jù)及其背后各個(gè)狀態(tài)轉(zhuǎn)換的過(guò)程。狀態(tài)α 包含 Local trend:（局部趨勢(shì)); local level:(局部趨勢(shì)的均值) 以及協(xié)變量x,表示所有的觀測(cè)數(shù)據(jù);表示根據(jù)狀態(tài)模型得到的預(yù)測(cè)數(shù)據(jù)。分別表示的標(biāo)準(zhǔn)差這些參數(shù)均通過(guò)MCMC的方式得到估計(jì)。

四、模型應(yīng)用與代碼實(shí)現(xiàn)

以上我們給出了BSTS模型及MCMC方法的簡(jiǎn)要理論推導(dǎo)及結(jié)果輸出，核心目的就是對(duì)觀測(cè)值y做出預(yù)測(cè)。接下來(lái)我們將介紹如何在因果推斷場(chǎng)景中應(yīng)用BSTS模型。

在對(duì)政策的效果評(píng)估上，我們核心想要的是觀測(cè)對(duì)象“反事實(shí)值”，例如“如果沒(méi)有這個(gè)廣告投放，用戶(hù)的瀏覽情況會(huì)怎樣？”相較于傳統(tǒng)的PSM或SCM方法，BSTS勝在其能夠?qū)τ跁r(shí)間序列數(shù)據(jù)進(jìn)行效果評(píng)估；同時(shí)利用貝葉斯估計(jì)輸出反事實(shí)值y的預(yù)測(cè)，并給出預(yù)測(cè)值的置信區(qū)間，能一定程度上降低反事實(shí)值預(yù)測(cè)的波動(dòng)性，提升效應(yīng)評(píng)估的準(zhǔn)確性與穩(wěn)定性。

在實(shí)踐應(yīng)用上，可以通過(guò)谷歌開(kāi)源的CausalImpact包來(lái)實(shí)現(xiàn)BSTS模型，在Python和R中均可調(diào)用，具體代碼實(shí)現(xiàn)詳見(jiàn)參考文獻(xiàn)[7][8]。

圖4-1：展示執(zhí)行代碼時(shí)的三個(gè)步驟：訓(xùn)練BSTS模型；反事實(shí)值預(yù)測(cè)；計(jì)算因果效應(yīng)值，包括效應(yīng)值的點(diǎn)估計(jì)及置信區(qū)間。

4.1 代碼實(shí)現(xiàn)

下面通過(guò)模擬數(shù)據(jù)展示代碼的具體命令

import tensorflow as tf
import tensorflow_probability as tfp
from causalimpact import CausalImpact
# 模型初始化 - 自定義時(shí)間序列數(shù)據(jù)：
def plot_time_series_components(ci):
     component_dists = tfp.sts.decompose_by_component(ci.model, ci.observed_time_series, ci.model_samples)
       num_components = len(component_dists)
mu, sig = ci.mu_sig if ci.mu_sig is not None else 0.0, 1.0
for i, (component, component_dist) in enumerate(component_dists.items()):
         component_mean = component_dist.mean().numpy()
         component_stddev = component_dist.stddev().numpy()
# 自定義觀測(cè)方程以及真實(shí)值y：
def plot_forecast_components(ci):
         component_forecasts = tfp.sts.decompose_forecast_by_component(ci.model, ci.posterior_dist, ci.model_samples)
       num_components = len(component_forecasts)
       mu, sig = ci.mu_sig if ci.mu_sig is not None else 0.0, 1.0
       for i, (component, component_dist) in enumerate(component_forecasts.items()):
          component_mean = component_dist.mean().numpy()
          component_stddev = component_dist.stddev().numpy()
# 生成模擬數(shù)據(jù)，包括一個(gè)實(shí)驗(yàn)組數(shù)據(jù)（有干預(yù)）以及兩條對(duì)照組數(shù)據(jù)（無(wú)干預(yù)）
observed_stddev, observed_initial = (tf.convert_to_tensor(value=1, dtype=tf.float32),tf.convert_to_tensor(value=0., dtype=tf.float32))
level_scale_prior = tfd.LogNormal(loc=tf.math.log(0.05 * observed_stddev), scale=1, name='level_scale_prior')  # 設(shè)置先驗(yàn)分布
initial_state_prior = tfd.MultivariateNormalDiag(loc=observed_initial[..., tf.newaxis], scale_diag=(tf.abs(observed_initial) + observed_stddev)[..., tf.newaxis], name='initial_level_prior')  # 設(shè)置先驗(yàn)分布 
ll_ssm = tfp.sts.LocalLevelStateSpaceModel(100, initial_state_prior=initial_state_prior,  level_scale=level_scale_prior.sample())  #訓(xùn)練時(shí)序模型
ll_ssm_sample = np.squeeze(ll_ssm.sample().numpy())
# 整合數(shù)據(jù)
x0 = 100 * np.random.rand(100)    # 對(duì)照組1
x1 = 90 * np.random.rand(100)     # 對(duì)照組2
y = 1.2 * x0 + 0.9 * x1 + ll_ssm_sample    #生成真實(shí)值y
y[70:] += 10     #設(shè)置干預(yù)點(diǎn)
data = pd.DataFrame({'x0': x0, 'x1': x1, 'y': y}, columns=['y', 'x0', 'x1'])

圖4-2：展示模擬數(shù)據(jù)。虛線表示干預(yù)發(fā)生的時(shí)間點(diǎn)，藍(lán)線表示受到干預(yù)的觀測(cè)數(shù)據(jù)；黃線與綠線表示沒(méi)有受到干預(yù)的兩組對(duì)照數(shù)據(jù)。

# 調(diào)用模型：
pre_period = [0, 69]    #設(shè)置干預(yù)前的時(shí)間窗口
post_period = [70, 99]  #干預(yù)后的窗口
ci = CausalImpact(data, pre_period, post_period)  #調(diào)用CausalImpact
# 對(duì)于causalImpact的使用我們核心需要填寫(xiě)三個(gè)參數(shù)：觀測(cè)數(shù)據(jù)data、干預(yù)前的時(shí)間窗口、干預(yù)后的時(shí)間窗口。
# 輸出結(jié)果：
ci.plot()
ci.summary()圖4-3：展示CausalImpact輸出的結(jié)果圖，圖1表示真實(shí)值與模型擬合值的曲線；圖2表示每個(gè)時(shí)刻真實(shí)值與預(yù)測(cè)值的差異；圖3表示真實(shí)值與預(yù)測(cè)值的累計(jì)差值。表3-1：展示CausalImpact輸出的結(jié)果表格，量化效應(yīng)值effect的估計(jì)及其置信區(qū)間，反映效應(yīng)值是否具有顯著性。107.71表示干預(yù)之后實(shí)際值的平均；96.25表示干預(yù)之后預(yù)測(cè)值的平均，3.28表示估計(jì)的標(biāo)準(zhǔn)差，[89.77,102.64]表示反事實(shí)估計(jì)的置信區(qū)間。11.46表示實(shí)際值與預(yù)測(cè)值的差距，[5.07,17.94]表示差值的置信區(qū)間，由于差距的置信區(qū)間在0的右側(cè)，表示干預(yù)有顯著的提升作用。

圖4-3：展示CausalImpact輸出的結(jié)果圖，圖1表示真實(shí)值與模型擬合值的曲線；圖2表示每個(gè)時(shí)刻真實(shí)值與預(yù)測(cè)值的差異，橙色陰影部分表示置信區(qū)間；圖3表示真實(shí)值與預(yù)測(cè)值的累計(jì)差值。

表3-1：展示CausalImpact輸出的結(jié)果表格，量化效應(yīng)值effect的估計(jì)及其置信區(qū)間，反映效應(yīng)值是否具有顯著性。107.71表示干預(yù)之后實(shí)際值的平均；96.25表示干預(yù)之后預(yù)測(cè)值的平均，3.28表示估計(jì)的標(biāo)準(zhǔn)差，[89.77,102.64]表示反事實(shí)估計(jì)的置信區(qū)間。11.46表示實(shí)際值與預(yù)測(cè)值的差距，[5.07,17.94]表示差值的置信區(qū)間，由于差距的置信區(qū)間在0的右側(cè)，表示干預(yù)有顯著的提升作用。

4.2 模型校驗(yàn)

對(duì)于模型擬合的結(jié)果，我們需要進(jìn)行類(lèi)似AB實(shí)驗(yàn)的“AA校驗(yàn)”。一般可以通過(guò)圖示的結(jié)果中的第二張圖，觀察干預(yù)之前真實(shí)值與預(yù)測(cè)值差值的置信區(qū)間是否包含0，如果包含0則說(shuō)明通過(guò)檢驗(yàn)，模型擬合效果不錯(cuò)。上圖中，置信區(qū)間均含0，說(shuō)明模型可用。

4.3 模型調(diào)整

• 過(guò)程參數(shù)：我們可以使用Tensorflow中的Decomposition來(lái)查看時(shí)序模型中各個(gè)結(jié)構(gòu)元素，包括周期性/季節(jié)性等等。

seasonal_decompose(data)

圖4-4展示了生成數(shù)據(jù)背后的狀態(tài)元素。第一張圖反映原數(shù)據(jù)的走勢(shì)；第二張圖反映局部趨勢(shì)因子；第三張圖反映季節(jié)性因子?？梢钥闯鰯?shù)據(jù)存在季節(jié)性結(jié)構(gòu)且呈單調(diào)上升趨勢(shì)。

• 自定義參數(shù)：我們可以自定義參數(shù)的先驗(yàn)分布；迭代次數(shù)；周期性的時(shí)間窗口長(zhǎng)度等等。往往參數(shù)調(diào)整會(huì)對(duì)結(jié)果輸出有影響，例如正確的選取先驗(yàn)分布會(huì)讓結(jié)果更準(zhǔn)確；迭代次數(shù)更多能保證MCMC收斂更穩(wěn)定（但也可能導(dǎo)致模型運(yùn)行時(shí)間較長(zhǎng)）等等。其中最重要的是對(duì)時(shí)間窗口長(zhǎng)度的設(shè)置，需要正確地反映觀測(cè)數(shù)據(jù)的周期性。如果是年維度數(shù)據(jù)以星期為周期則設(shè)置neasnotallow=52；如果是天維度數(shù)據(jù)以小時(shí)為周期則設(shè)置neasnotallow=24等等。

CausalImpact(...,  model.args = list(niter = 20000, nseasons = 24))

圖4-5展示CausalImpact包中各個(gè)參數(shù)含義及其默認(rèn)值。

• 自定義時(shí)序模型：causalImpact的包中默認(rèn)使用BSTS模型進(jìn)行訓(xùn)練，我們也可以改為其他的時(shí)序模型，但前提是需要對(duì)數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化。（如果使用默認(rèn)的BSTS則不一定需要標(biāo)準(zhǔn)化）

from causalimpact.misc import standardize
normed_data, _ = standardize(data.astype(np.float32)) #標(biāo)準(zhǔn)化數(shù)據(jù)
obs_data = normed_data.iloc[:70, 0]
# 使用tfp中的其他模型來(lái)訓(xùn)練時(shí)序數(shù)據(jù)
linear_level = tfp.sts.LocalLinearTrend(observed_time_series=obs_data)
linear_reg = tfp.sts.LinearRegression(design_matrix=normed_data.iloc[:, 1:].values.reshape(-1, normed_data.shape[1] -1))
model = tfp.sts.Sum([linear_level, linear_reg], observed_time_series=obs_data)
# 將自義定時(shí)序模型代入CausalImpact包中
ci = CausalImpact(data, pre_period, post_period, model=model)

五、業(yè)務(wù)場(chǎng)景實(shí)踐

用戶(hù)營(yíng)銷(xiāo)是促進(jìn)留存及轉(zhuǎn)化的重要方式，其中對(duì)用戶(hù)進(jìn)行消息觸達(dá)是一大核心手段，尤其是在節(jié)假日的購(gòu)票高峰期對(duì)用戶(hù)進(jìn)行推送，方式包括站內(nèi)push；微信生態(tài)中的小程序訂閱消息；公眾號(hào)或是企微環(huán)境等等，目的包括但不限于提醒用戶(hù)購(gòu)票；宣傳品牌功能；發(fā)放優(yōu)惠券吸引用戶(hù)轉(zhuǎn)化等等。在節(jié)假日之后，我們希望對(duì)這次的營(yíng)銷(xiāo)觸達(dá)進(jìn)行效果評(píng)估。

這是一個(gè)較為典型的不適合進(jìn)行AB實(shí)驗(yàn)的場(chǎng)景，首先因?yàn)楣?jié)假日是流量高峰時(shí)期，如果嚴(yán)格預(yù)留50%用戶(hù)不觸達(dá)，可能會(huì)損失一批潛在的轉(zhuǎn)化用戶(hù)；如果改將對(duì)照組預(yù)留很少的人數(shù)，例如對(duì)照組：實(shí)驗(yàn)組=1:9，那對(duì)于后續(xù)的轉(zhuǎn)化對(duì)比的科學(xué)性會(huì)產(chǎn)生影響。其次，節(jié)假日的推送策略往往非常精細(xì)化，總量達(dá)幾十條，我們較難保證對(duì)照組用戶(hù)的“純粹性”，用戶(hù)可能會(huì)被交叉觸達(dá)。

基于種種問(wèn)題，我們較難通過(guò)傳統(tǒng)AB 的手段來(lái)評(píng)估推送帶來(lái)的轉(zhuǎn)化效果。因此我們考慮使用因果推斷的方式來(lái)解決。常規(guī)可選的方法和潛在問(wèn)題如下：

• 如果使用PSM，需要在大盤(pán)中尋找與推送人群相似但是沒(méi)有被推送的用戶(hù)作為對(duì)照組。但一般節(jié)假日推送時(shí)都會(huì)有兜底策略，幾乎覆蓋了95%以上的平臺(tái)用戶(hù)，較難從中找到符合條件但未被推送的人群來(lái)進(jìn)行對(duì)照。
• 如果使用SCM，我們較難找到合適的對(duì)照組來(lái)合成。如評(píng)估度假BU的推送效果時(shí)，我們不太可能用火車(chē)、機(jī)票、酒店等各個(gè)產(chǎn)線合成一個(gè)“虛擬度假BU”，因?yàn)楸旧砀鱾€(gè)產(chǎn)線的用戶(hù)需求就不同，使用這樣合成的虛擬對(duì)照組來(lái)對(duì)比度假訂單的轉(zhuǎn)化率是不夠科學(xué)的。
• DID的方式也同理，我們很難找到一個(gè)滿足平行趨勢(shì)假設(shè)且業(yè)務(wù)場(chǎng)景相似的對(duì)照組來(lái)進(jìn)行推送前后的對(duì)比。

綜上所述，一些傳統(tǒng)的因果推斷方法縱使在技術(shù)上可行，在業(yè)務(wù)的解釋性上也有所欠缺。而且，以上三種方式都沒(méi)有考慮到用戶(hù)購(gòu)票行為的“時(shí)間周期性”。因此即使合成了對(duì)照組也不一定能夠匹配到實(shí)驗(yàn)組真正的結(jié)構(gòu)特點(diǎn)，進(jìn)而導(dǎo)致效應(yīng)值計(jì)算有偏。于是我們考慮首先驗(yàn)證用戶(hù)購(gòu)票的數(shù)據(jù)周期性；在定位到周期規(guī)律之后嘗試使用BSTS模型結(jié)合CausalImpact來(lái)進(jìn)行反事實(shí)值的預(yù)測(cè)。下文我們選擇2022年端午的火車(chē)票營(yíng)銷(xiāo)推送場(chǎng)景進(jìn)行實(shí)踐。

圖5-1展示端午期間對(duì)于用戶(hù)進(jìn)行不同策略的推送觸達(dá)。

5.1 數(shù)據(jù)選取

• 我們以小時(shí)為周期窗口，通過(guò)簡(jiǎn)單的圖像能夠看出大盤(pán)的火車(chē)票下單人數(shù)確實(shí)隨著時(shí)間推移呈現(xiàn)某種固定趨勢(shì)。

圖5-2展示選取端午周期內(nèi)（正端午前后10天）每小時(shí)的火車(chē)票大盤(pán)支付人數(shù)

• 考慮到端午作為節(jié)假日本身就有的自然流量增長(zhǎng)，支付人數(shù)的提升不能完全歸因于推送帶來(lái)的，因此訓(xùn)練時(shí)序模型的時(shí)候，選取了19年-21年所有的端午數(shù)據(jù)（正端午前后10天）輸入BSTS模型進(jìn)行訓(xùn)練，得到端午這個(gè)窗口內(nèi)的特有的結(jié)構(gòu)狀態(tài)，隨后用這個(gè)結(jié)構(gòu)化的模型來(lái)代入22年的端午數(shù)據(jù)，對(duì)2022年端午推送之后的轉(zhuǎn)化人數(shù)做出預(yù)測(cè)。
• 最后使用真實(shí)的轉(zhuǎn)化人數(shù)與預(yù)測(cè)人數(shù)作差體現(xiàn)本次營(yíng)銷(xiāo)推送的效果。

5.2 R-代碼實(shí)現(xiàn)

# 選取19-22年每年的端午窗口，按照小時(shí)劃分，共960個(gè)數(shù)據(jù)點(diǎn)
y_hour=c(x1,x2,x3,x4)
x_time_hour=c(1:960)
data_hour <- cbind(y_hour, x_time_hour)
pre.period <- c(1, 808)   # 2022年的推送發(fā)生在第808個(gè)時(shí)間點(diǎn)，故以此為干預(yù)節(jié)點(diǎn)。
post.period <- c(809, 960)
# nseasnotallow=24, 迭代次數(shù)2000，fit the model
impact_hour <- CausalImpact(data_hour, pre.period, post.period, model.args = list(niter = 20000, nseasons = 24))
summary(impact_hour)
plot(impact_hour)

圖5-3展示使用CausalImpact返回的結(jié)果圖。第一張圖表示真實(shí)支付人數(shù)與預(yù)測(cè)支付人數(shù)；第二張圖表示真實(shí)值與預(yù)測(cè)值的差值及置信區(qū)間；第三張圖是累計(jì)差值和置信區(qū)間。

圖像顯示模型能夠通過(guò) AA校驗(yàn)，模型有效。在干預(yù)點(diǎn)之后，實(shí)際值較預(yù)測(cè)值有所提升，但提升的置信區(qū)間含0，因此未達(dá)到顯著程度。體現(xiàn)出2022年端午營(yíng)銷(xiāo)策略對(duì)于轉(zhuǎn)化人數(shù)有所促進(jìn)作用，但是效果未達(dá)顯著。

六、方法優(yōu)缺點(diǎn)

相較于傳統(tǒng)因果推斷方法，BSTS模型有2個(gè)主要優(yōu)點(diǎn)：

• 能夠識(shí)別出數(shù)據(jù)背后的結(jié)構(gòu)化特征，更好的做出預(yù)測(cè)；
• 利用了貝葉斯估計(jì)的思想，得到參數(shù)的后驗(yàn)分布情況，計(jì)算效應(yīng)值時(shí)能夠給出置信區(qū)間。但第（2）點(diǎn)對(duì)于BSTS模型是一把“雙刃劍”，如果先驗(yàn)分布設(shè)置得不好，會(huì)影響MCMC的收斂速度和方向甚至最終的后驗(yàn)分布情況。因此對(duì)于先驗(yàn)分布的選取需謹(jǐn)慎。

七、方法拓展

本文介紹的結(jié)構(gòu)化時(shí)序模型將數(shù)據(jù)的周期特點(diǎn)拆分成了趨勢(shì)項(xiàng)、季節(jié)項(xiàng)、周期項(xiàng)等等，每種元素挨個(gè)探究。更進(jìn)一步，我們可以將時(shí)間序列按照周期性的長(zhǎng)短來(lái)進(jìn)行拆分，分為長(zhǎng)周期項(xiàng)（使用大滑動(dòng)窗口）、短周期項(xiàng)（使用小滑動(dòng)窗口）、季節(jié)項(xiàng)等等。這樣的好處是防止一些小窗口內(nèi)的周期情況被長(zhǎng)周期的信息平滑掉，能夠更好的體現(xiàn)出數(shù)據(jù)在不同程度上的周期特點(diǎn)。具體的方程可以拆分成如下形式：

其中表示不同時(shí)間點(diǎn)的狀態(tài)值；4個(gè)模塊分別代表長(zhǎng)周期項(xiàng)/短周期項(xiàng)/季節(jié)項(xiàng)/序列相關(guān)性項(xiàng)（帶有協(xié)變量X）；每個(gè)結(jié)構(gòu)模塊都有一個(gè)均值和一個(gè)標(biāo)準(zhǔn)差。

圖7-1展示了某時(shí)間序列背后的4個(gè)模塊：從上至下以此表示：原數(shù)據(jù)情況；季節(jié)性因子；短周期項(xiàng)；相關(guān)性項(xiàng)；長(zhǎng)周期項(xiàng)。短周期來(lái)看數(shù)據(jù)的波動(dòng)比較明顯；長(zhǎng)周期來(lái)看數(shù)據(jù)波動(dòng)較不明顯，因此這里需要考慮到短周期內(nèi)的數(shù)據(jù)結(jié)構(gòu)，避免被長(zhǎng)周期的數(shù)據(jù)平滑。

接下來(lái)對(duì)以上4種模塊分別進(jìn)行預(yù)測(cè)。針對(duì)長(zhǎng)周期和季節(jié)性，由于他們?cè)诙虝r(shí)間內(nèi)的變化不大，因此可以直接使用對(duì)應(yīng)方程中的μ和σ來(lái)進(jìn)行預(yù)測(cè)；針對(duì)短周期項(xiàng)和相關(guān)性項(xiàng)可以通過(guò)其他機(jī)器學(xué)習(xí)方式進(jìn)行預(yù)測(cè)。得到各個(gè)模塊的預(yù)測(cè)結(jié)果之后，結(jié)合各模塊特征進(jìn)行融合，得到整體的預(yù)測(cè)結(jié)果。參考文獻(xiàn)[4]中給出了更具體的預(yù)測(cè)方式和與傳統(tǒng)方式的對(duì)比結(jié)果。

圖7-2展示針對(duì)長(zhǎng)短周期不同的預(yù)測(cè)方式：長(zhǎng)周期項(xiàng)與季節(jié)項(xiàng)可以直接用μ表示預(yù)測(cè)；短周期及協(xié)變量相關(guān)項(xiàng)使用自定義的機(jī)器學(xué)習(xí)模型進(jìn)行預(yù)測(cè)。

依照上述方法得到的時(shí)序預(yù)測(cè)模型后，我們?cè)賹⑵浯隒ausalImpact的代碼中，調(diào)整參數(shù)model為自定義的時(shí)序模型即可。

八、總結(jié)

本文介紹了用因果推斷的方式評(píng)估某一政策作用于時(shí)間序列數(shù)據(jù)帶來(lái)的效應(yīng)，使用了BSTS的狀態(tài)空間模型來(lái)進(jìn)行反事實(shí)值的預(yù)測(cè)，并通過(guò)CausalImpact的代碼進(jìn)行實(shí)現(xiàn)。

不同于其他因果推斷方法的框架，本文中的方法對(duì)所有超參數(shù)進(jìn)行貝葉斯估計(jì)，再對(duì)后驗(yàn)分布進(jìn)行MC采樣得到反事實(shí)的預(yù)測(cè)值，主要優(yōu)勢(shì)是能夠根據(jù)最大程度考慮到所有對(duì)照組以及實(shí)驗(yàn)組自身的結(jié)構(gòu)特點(diǎn)，給出反事實(shí)預(yù)測(cè)值的估計(jì)以及置信區(qū)間，衡量效應(yīng)值的顯著性。

同時(shí)，本文介紹的方法主要聚焦于結(jié)構(gòu)化時(shí)序數(shù)據(jù)，利用BSTS模型識(shí)別觀測(cè)數(shù)據(jù)背后的狀態(tài)值以及各個(gè)狀態(tài)之間的轉(zhuǎn)化情況，進(jìn)而在進(jìn)行反事實(shí)預(yù)測(cè)時(shí)，盡可能消除由隱藏狀態(tài)帶來(lái)的影響。

需要注意的是，使用BSTS模型之前，需要驗(yàn)證數(shù)據(jù)是否真的有周期性特點(diǎn)以及結(jié)構(gòu)元素是怎樣的（是否是長(zhǎng)短周期等等），再挑選合適的時(shí)序模型來(lái)訓(xùn)練；同時(shí)對(duì)于參數(shù)的先驗(yàn)分布設(shè)置也需謹(jǐn)慎，盡可能使得最終的效應(yīng)估計(jì)值科學(xué)穩(wěn)定。