作者简介
Yiwen，携程数据分析师，专注用户增长、因果推断、数据科学等领域。
一、背景
如何科学地推断某个产品策略对观测指标产生的效应非常重要，这能够帮助产品和运营更精准地得到该策略的价值，从而进行后续方向的迭代及调整。
在因果推断框架下，效果评估的黄金准则一定是“AB实验”，因为实验的分流被认为是完全随机且均匀的，在此基础上对比实验组与对照组的指标差异就可以体现某个干预带来的增量值。但是很多场景下，我们较难进行严格的AB实验，例如对于酒店的定价；现金奖励的发放等等，不适宜向不同人群展现不同的内容。对于这些问题，我们会采取因果推断的方法来进行策略的效果评估。
本文主要介绍BSTS模型原理以及CausalImpact对模型的代码实现，旨在面对一些具有特定周期性特点的数据时，更精准科学地进行因果效应值的估计。下文将首先对模型原理进行简要阐释；随后利用模拟数据展示代码逻辑，最后在具体的业务场景中进行实践。
二、现有方法及潜在问题
大部分运营和产品在评估一些全量上线的策略效果时，最常用的方式就是看上线前后的效果差异。但这种方法最大的问题在于其假设前提：假设上线的功能是唯一影响效果的变量（即没有任何其他干预和混淆变量），但这个假设现实中往往很难实现。
于是我们尝试使用更多因果推断的方法，例如PSM（倾向分匹配法），在所有非实验组的用户群中，找到与实验组用户的特征非常相似的一群人，将他们的指标数据（例如下单率，订单收益等等）与实验组的用户进行对比，从而体现出干预带来的影响。但这个方法较为依赖选取的用户特征与最后的匹配效果。
再例如SCM（合成控制方法），利用一些未受干预的地区合成一个“类似的虚拟地区”来与“上线策略的地区”进行整体的对比。但这也需要一个关键假设：可以找到长期变化趋势高度同步的地区来进行对照，而这个条件往往也很难实现。
进而在传统SCM的基础上，我们企图通过类似集成学习的方法，将多个未干预的对照组作为输入值，再结合实验组自身长期的时间序列波动情况，拟合出一个未受干预的虚拟对照组，从而将“对照组与实验组高度同步”的强假设降为弱假设。本文介绍的BSTS模型就是用来刻画某种“长期的时间序列波动”的数据模型，CausalImpact是用来针对这样的数据进行因果效应值的估计。下文中我们将详细介绍这两个工具。
三、模型介绍
BSTS模型 （Bayesian Structured Time Series）称为“贝叶斯结构化时间序列”，正如其名，它的主要特点体现在：
- 适用于有结构特征的时间序列数据
- 利用贝叶斯的思想来进行参数估计
结构化的时间序列数据在日常生活中不少见，尤其像携程这样的OTA行业，平台的订单情况其实是有一定时间规律的，例如周末和节假日是订单高峰期；周中是订单平峰期等。另一方面，贝叶斯的思想是指在得到样本数据之前，即对要估计的参数有一些“先天认知”），随后基于这样的认知，结合样本数据再得到后验分布（如下方公式展示）
故BSTS模型主要即对结构化时序数据进行模型拟合及预测，在拟合的过程中使用到了贝叶斯的先验思想。其好处是能够给出预测值的置信区间，使得预测结果更科学可信。下文将对这几种思想逐一进行介绍。
3.1 状态空间模型
结构化的时间序列数据是指某一观测数据的背后其实隐藏着随时间变化而变化的不同状态，其中观测值与状态值之间有对应关系；不同时刻的状态之间也有转换关系。我们一般用以下状态空间模型来刻画这两种映射逻辑：
(1) 称为观测方程，反映观测值与其背后隐藏状态的关系；(2) 称为状态方程，反映随时间推移各个状态之间的转换。；
都是不同变量之间的“关系映射矩阵”；
是独立于其他变量且服从正态分布的噪声。所谓数据的“结构化”，主要包括：
- Linear Local Trend（局部趋势）：一定时间内的单调性（单调上升或下降）
- Seasonality（季节性因子）：固定长度的变化，类似于一年四季的温度变化
- Cyclical（周期性）：类似季节性但波动时间不固定，波动频率也不固定的变化
图3-1：观测数据及其结构化元素。第一张图体现原数据的波动情况；第二张体现季节性因子的情况；第三张图体现局部趋势的情况。
如果希望在映射关系中加入协变量X，可以将(1)拓展为：
其中
表示协变量X与观测数据之间的关系，如果协变量项表现很好（如有显著影响）的话，那对应的local trend就会相对较弱。上述三个方程中的参数将在后文中展示估计方式。
3.2 贝叶斯及MCMC（马尔可夫蒙特卡洛方法）
假设状态方程(2)中各个时刻的状态序列为
表示模型中所有的参数。我们现在希望对θ进行估计，核心步骤如下：
- 对θ设置先验分布
以及初始状态的分布 
- 构造马尔科夫链，用MCMC方法得到
- 通过贝叶斯公式计算得到参数的后验分布
下面对于各个步骤中用到的方法进行简要说明：
1）贝叶斯估计：BSTS模型的一大特点就是在参数估计上使用了贝叶斯估计，即在估计之前先给出参数设置先验分布，随后再结合样本数据给出参数的后验分布。不同类型的参数一般有一些常用的先验分布，例如均值一般使用正态分布，
，方差使用inverse-Gamma分布，
协方差矩阵可以使用IW分布等等。值得注意的是，先验分布的设置一定程度上会影响后续MCMC收敛的情况以及后验分布的准确性，因此并不能太过随意地设置先验分布，应尽可能多地根据实际数据推导出最合适的先验分布，或是比较各先验分布下后验分布和似然函数的值来进行选择。
2）MCMC方法：我们尝试构造一条马尔可夫链（一种特殊的序列，当前时刻的状态值仅与前一时刻的状态值有关，最终序列会收敛到某个稳定的分布），使得其最终收敛的稳态分布就是参数的后验分布。这一过程我们可以通过Gibbs采样实现：设置先验分布之后，从初始状态
出发，每次固定α采样θ；再固定θ采样α，逐渐一次次更新两组参数，最终形成一条服从马尔可夫性质的链路
，可以证明其稳态收敛的分布就是
，其中
代表所有的观测数据。
3）预测值估计：得到
之后，我们从该分布中对（α，θ）进行采样，再代入状态空间方程(1)中对y进行预测，得到
，其中
表示时间点n之后y的预测值。
图3-2：展示某结构化时序数据及其背后各个状态转换的过程。状态α 包含 Local trend:
（局部趋势); local level:
(局部趋势的均值) 以及协变量x,
表示所有的观测数据;
表示根据状态模型得到的预测数据。
分别表示
的标准差这些参数均通过MCMC的方式得到估计。
四、模型应用与代码实现
以上我们给出了BSTS模型及MCMC方法的简要理论推导及结果输出，核心目的就是对观测值y做出预测。接下来我们将介绍如何在因果推断场景中应用BSTS模型。
在对政策的效果评估上，我们核心想要的是观测对象“反事实值”，例如“如果没有这个广告投放，用户的浏览情况会怎样？”相较于传统的PSM或SCM方法，BSTS胜在其能够对于时间序列数据进行效果评估；同时利用贝叶斯估计输出反事实值y的预测，并给出预测值的置信区间，能一定程度上降低反事实值预测的波动性，提升效应评估的准确性与稳定性。
在实践应用上，可以通过谷歌开源的CausalImpact包来实现BSTS模型，在Python和R中均可调用，具体代码实现详见参考文献[7][8]。
图4-1：展示执行代码时的三个步骤：训练BSTS模型；反事实值预测；计算因果效应值，包括效应值的点估计及置信区间。
4.1 代码实现
下面通过模拟数据展示代码的具体命令
import tensorflow as tf
import tensorflow_probability as tfp
from causalimpact import CausalImpact
# 模型初始化 - 自定义时间序列数据：
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()
# 自定义观测方程以及真实值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()
# 生成模拟数据，包括一个实验组数据（有干预）以及两条对照组数据（无干预）
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') # 设置先验分布
initial_state_prior = tfd.MultivariateNormalDiag(loc=observed_initial[..., tf.newaxis], scale_diag=(tf.abs(observed_initial) + observed_stddev)[..., tf.newaxis], name='initial_level_prior') # 设置先验分布
ll_ssm = tfp.sts.LocalLevelStateSpaceModel(100, initial_state_prior=initial_state_prior, level_scale=level_scale_prior.sample()) #训练时序模型
ll_ssm_sample = np.squeeze(ll_ssm.sample().numpy())
# 整合数据
x0 = 100 * np.random.rand(100) # 对照组1
x1 = 90 * np.random.rand(100) # 对照组2
y = 1.2 * x0 + 0.9 * x1 + ll_ssm_sample #生成真实值y
y[70:] += 10 #设置干预点
data = pd.DataFrame({'x0': x0, 'x1': x1, 'y': y}, columns=['y', 'x0', 'x1'])
图4-2：展示模拟数据。虚线表示干预发生的时间点，蓝线表示受到干预的观测数据；黄线与绿线表示没有受到干预的两组对照数据。
# 调用模型：
pre_period = [0, 69] #设置干预前的时间窗口
post_period = [70, 99] #干预后的窗口
ci = CausalImpact(data, pre_period, post_period) #调用CausalImpact
# 对于causalImpact的使用我们核心需要填写三个参数：观测数据data、干预前的时间窗口、干预后的时间窗口。
# 输出结果：
ci.plot()
ci.summary()图4-3：展示CausalImpact输出的结果图，图1表示真实值与模型拟合值的曲线；图2表示每个时刻真实值与预测值的差异；图3表示真实值与预测值的累计差值。表3-1：展示CausalImpact输出的结果表格，量化效应值effect的估计及其置信区间，反映效应值是否具有显著性。107.71表示干预之后实际值的平均；96.25表示干预之后预测值的平均，3.28表示估计的标准差，[89.77,102.64]表示反事实估计的置信区间。11.46表示实际值与预测值的差距，[5.07,17.94]表示差值的置信区间，由于差距的置信区间在0的右侧，表示干预有显著的提升作用。
图4-3：展示CausalImpact输出的结果图，图1表示真实值与模型拟合值的曲线；图2表示每个时刻真实值与预测值的差异，橙色阴影部分表示置信区间；图3表示真实值与预测值的累计差值。
表3-1：展示CausalImpact输出的结果表格，量化效应值effect的估计及其置信区间，反映效应值是否具有显著性。107.71表示干预之后实际值的平均；96.25表示干预之后预测值的平均，3.28表示估计的标准差，[89.77,102.64]表示反事实估计的置信区间。11.46表示实际值与预测值的差距，[5.07,17.94]表示差值的置信区间，由于差距的置信区间在0的右侧，表示干预有显著的提升作用。
4.2 模型校验
对于模型拟合的结果，我们需要进行类似AB实验的“AA校验”。一般可以通过图示的结果中的第二张图，观察干预之前真实值与预测值差值的置信区间是否包含0，如果包含0则说明通过检验，模型拟合效果不错。上图中，置信区间均含0，说明模型可用。
4.3 模型调整
- 过程参数：我们可以使用Tensorflow中的Decomposition来查看时序模型中各个结构元素，包括周期性/季节性等等。
seasonal_decompose(data)
图4-4展示了生成数据背后的状态元素。第一张图反映原数据的走势；第二张图反映局部趋势因子；第三张图反映季节性因子。可以看出数据存在季节性结构且呈单调上升趋势。
- 自定义参数：我们可以自定义参数的先验分布；迭代次数；周期性的时间窗口长度等等。往往参数调整会对结果输出有影响，例如正确的选取先验分布会让结果更准确；迭代次数更多能保证MCMC收敛更稳定（但也可能导致模型运行时间较长）等等。其中最重要的是对时间窗口长度的设置，需要正确地反映观测数据的周期性。如果是年维度数据以星期为周期则设置neasnotallow=52；如果是天维度数据以小时为周期则设置neasnotallow=24等等。
CausalImpact(..., model.args = list(niter = 20000, nseasons = 24))
图4-5展示CausalImpact包中各个参数含义及其默认值。
- 自定义时序模型：causalImpact的包中默认使用BSTS模型进行训练，我们也可以改为其他的时序模型，但前提是需要对数据进行标准化。（如果使用默认的BSTS则不一定需要标准化）
from causalimpact.misc import standardize
normed_data, _ = standardize(data.astype(np.float32)) #标准化数据
obs_data = normed_data.iloc[:70, 0]
# 使用tfp中的其他模型来训练时序数据
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)
# 将自义定时序模型代入CausalImpact包中
ci = CausalImpact(data, pre_period, post_period, model=model)五、业务场景实践
用户营销是促进留存及转化的重要方式，其中对用户进行消息触达是一大核心手段，尤其是在节假日的购票高峰期对用户进行推送，方式包括站内push；微信生态中的小程序订阅消息；公众号或是企微环境等等，目的包括但不限于提醒用户购票；宣传品牌功能；发放优惠券吸引用户转化等等。在节假日之后，我们希望对这次的营销触达进行效果评估。
这是一个较为典型的不适合进行AB实验的场景，首先因为节假日是流量高峰时期，如果严格预留50%用户不触达，可能会损失一批潜在的转化用户；如果改将对照组预留很少的人数，例如对照组：实验组=1:9，那对于后续的转化对比的科学性会产生影响。其次，节假日的推送策略往往非常精细化，总量达几十条，我们较难保证对照组用户的“纯粹性”，用户可能会被交叉触达。
基于种种问题，我们较难通过传统AB 的手段来评估推送带来的转化效果。因此我们考虑使用因果推断的方式来解决。常规可选的方法和潜在问题如下：
- 如果使用PSM，需要在大盘中寻找与推送人群相似但是没有被推送的用户作为对照组。但一般节假日推送时都会有兜底策略，几乎覆盖了95%以上的平台用户，较难从中找到符合条件但未被推送的人群来进行对照。
- 如果使用SCM，我们较难找到合适的对照组来合成。如评估度假BU的推送效果时，我们不太可能用火车、机票、酒店等各个产线合成一个“虚拟度假BU”，因为本身各个产线的用户需求就不同，使用这样合成的虚拟对照组来对比度假订单的转化率是不够科学的。
- DID的方式也同理，我们很难找到一个满足平行趋势假设且业务场景相似的对照组来进行推送前后的对比。
综上所述，一些传统的因果推断方法纵使在技术上可行，在业务的解释性上也有所欠缺。而且，以上三种方式都没有考虑到用户购票行为的“时间周期性”。因此即使合成了对照组也不一定能够匹配到实验组真正的结构特点，进而导致效应值计算有偏。于是我们考虑首先验证用户购票的数据周期性；在定位到周期规律之后尝试使用BSTS模型结合CausalImpact来进行反事实值的预测。下文我们选择2022年端午的火车票营销推送场景进行实践。
图5-1展示端午期间对于用户进行不同策略的推送触达。
5.1 数据选取
- 我们以小时为周期窗口，通过简单的图像能够看出大盘的火车票下单人数确实随着时间推移呈现某种固定趋势。
图5-2展示选取端午周期内（正端午前后10天）每小时的火车票大盘支付人数
- 考虑到端午作为节假日本身就有的自然流量增长，支付人数的提升不能完全归因于推送带来的，因此训练时序模型的时候，选取了19年-21年所有的端午数据（正端午前后10天）输入BSTS模型进行训练，得到端午这个窗口内的特有的结构状态，随后用这个结构化的模型来代入22年的端午数据，对2022年端午推送之后的转化人数做出预测。
- 最后使用真实的转化人数与预测人数作差体现本次营销推送的效果。
5.2 R-代码实现
# 选取19-22年每年的端午窗口，按照小时划分，共960个数据点
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年的推送发生在第808个时间点，故以此为干预节点。
post.period <- c(809, 960)
# nseasnotallow=24, 迭代次数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返回的结果图。第一张图表示真实支付人数与预测支付人数；第二张图表示真实值与预测值的差值及置信区间；第三张图是累计差值和置信区间。
图像显示模型能够通过 AA校验，模型有效。在干预点之后，实际值较预测值有所提升，但提升的置信区间含0，因此未达到显著程度。体现出2022年端午营销策略对于转化人数有所促进作用，但是效果未达显著。
六、方法优缺点
相较于传统因果推断方法，BSTS模型有2个主要优点：
- 能够识别出数据背后的结构化特征，更好的做出预测；
- 利用了贝叶斯估计的思想，得到参数的后验分布情况，计算效应值时能够给出置信区间。但第（2）点对于BSTS模型是一把“双刃剑”，如果先验分布设置得不好，会影响MCMC的收敛速度和方向甚至最终的后验分布情况。因此对于先验分布的选取需谨慎。
七、方法拓展
本文介绍的结构化时序模型将数据的周期特点拆分成了趋势项、季节项、周期项等等，每种元素挨个探究。更进一步，我们可以将时间序列按照周期性的长短来进行拆分，分为长周期项（使用大滑动窗口）、短周期项（使用小滑动窗口）、季节项等等。这样的好处是防止一些小窗口内的周期情况被长周期的信息平滑掉，能够更好的体现出数据在不同程度上的周期特点。具体的方程可以拆分成如下形式：
其中
表示不同时间点的状态值；4个模块分别代表长周期项/短周期项/季节项/序列相关性项（带有协变量X）；每个结构模块都有一个均值和一个标准差。
图7-1展示了某时间序列背后的4个模块：从上至下以此表示：原数据情况；季节性因子；短周期项；相关性项；长周期项。短周期来看数据的波动比较明显；长周期来看数据波动较不明显，因此这里需要考虑到短周期内的数据结构，避免被长周期的数据平滑。
接下来对以上4种模块分别进行预测。针对长周期和季节性，由于他们在短时间内的变化不大，因此可以直接使用对应方程中的μ和σ来进行预测；针对短周期项和相关性项可以通过其他机器学习方式进行预测。得到各个模块的预测结果之后，结合各模块特征进行融合，得到整体的预测结果。参考文献[4]中给出了更具体的预测方式和与传统方式的对比结果。
图7-2展示针对长短周期不同的预测方式：长周期项与季节项可以直接用μ表示预测；短周期及协变量相关项使用自定义的机器学习模型进行预测。
依照上述方法得到的时序预测模型后，我们再将其代入CausalImpact的代码中，调整参数model为自定义的时序模型即可。
八、总结
本文介绍了用因果推断的方式评估某一政策作用于时间序列数据带来的效应，使用了BSTS的状态空间模型来进行反事实值的预测，并通过CausalImpact的代码进行实现。
不同于其他因果推断方法的框架，本文中的方法对所有超参数进行贝叶斯估计，再对后验分布进行MC采样得到反事实的预测值，主要优势是能够根据最大程度考虑到所有对照组以及实验组自身的结构特点，给出反事实预测值的估计以及置信区间，衡量效应值的显著性。
同时，本文介绍的方法主要聚焦于结构化时序数据，利用BSTS模型识别观测数据背后的状态值以及各个状态之间的转化情况，进而在进行反事实预测时，尽可能消除由隐藏状态带来的影响。
需要注意的是，使用BSTS模型之前，需要验证数据是否真的有周期性特点以及结构元素是怎样的（是否是长短周期等等），再挑选合适的时序模型来训练；同时对于参数的先验分布设置也需谨慎，尽可能使得最终的效应估计值科学稳定。
