多层回归模型和辛普森悖论

避免使用不当工具得出错误的结论

数据分析是数据科学家工作的一个重要组成部分，从描述性统计和简单的回归模型到复杂的机器学习方法都涉及其中。然而，这些方法必须谨慎处理，选择正确的方法远非易事。复杂数据通常包含隐藏的结构，如果不适当考虑，可能导致谬误，最终导致无效的结论。

在本文中，我想举一个辛普森悖论的例子，展示一个简单但目光短浅的分析如何导致看似有数据支持的错误结论，实际上只是一种误解。同时，我演示了使用多层次回归模型作为层次性（即嵌套）数据分析的适当方式。

问题

让我们开始吧！假设我们发布了一个智能手机应用程序，想要了解更多关于我们的用户以及他们对我们的应用的满意程度。因此，我们进行了一项小型调查，询问一些用户在1（非常不满意）到4（非常满意）的评分尺度上对我们的应用的满意程度。此外，我们还测量了他们在过去一周内在我们的应用程序中所花费的时间，并为了得到丰富的样本，我们询问了来自不同国家的用户。那么我们的数据可能是这样的（本文使用的是生成的数据）：

   满意程度  时间花费  国家0      2.140440    1.585295        01      2.053545    0.636235        02      1.589258    1.468033        13      1.853545    0.968651        24      1.449286    0.967104        2.      .            .              ..      .            .              ..      .            .              .

我们对在我们的应用程序中花费的时间和报告的满意程度之间的关系感兴趣。确切地说，我们想知道在我们的应用程序中花费更多时间是否与更高或更低的满意度有关，并且我们希望量化这种关联，即我们想要提出这样的观点：“在我们的应用程序中多花费一小时与满意度增加或减少x倍有关”。当我们查看数据时，我们可能已经有了第一直觉，即在应用程序中花费更多时间与满意度较低相关：

线性回归

让我们进行线性回归，看看我们是否正确。通过线性回归，我们尝试预测满意度，给定花费时间作为形式为满意度=截距+回归系数*时间花费的线性函数。我们可以使用statsmodels包中的OLS（最小二乘法）函数轻松实现这一点。

import statsmodels.api as smresult = sm.OLS(df["满意程度"], sm.add_constant(df["时间花费"])).fit()print(result.params)

add_constant方法只是我们使用的一个技术细节，用于告诉模型我们想要在方程中有一个截距（这是必需的，只要我们的数据没有标准化）。result.params给我们两个值，即截距（const）和变量时间花费的回归系数。

const         3.229412时间花费   -0.655470

也就是说，我们的模型告诉我们，满意度可以预测为3.229-0.655*时间花费。换句话说，应用程序中再多花费一小时，满意度会减少0.655个点（因为有负号）。然而，起点不是从零开始，而是一个人的第一印象的平均满意度（即时间花费=0）为3.229。我们还可以将其可视化为截距为3.229，斜率为-0.665的线：

当然，这个预测并不完美，但至少给了我们一个趋势。好了，案例很明确了，对吧？在应用程序中花费更多的时间会导致满意度下降，我们甚至可以量化这种下降。我们现在可以从中得出结论，思考如何改进应用程序（当然，我们希望用户在使用应用程序时更加满意），或者进行更详细的调查，找出用户不满意的原因。

嗯，不要这么快下结论！

按国家分组

请记住，我们收集了来自不同国家的用户数据。如果我们按国家分别查看数据会发生什么呢？在下面的图表中，我们看到与之前相同的数据点，但现在我们用不同的颜色突出显示每个国家。

从该图表中我们可以得出两个观察结果。首先，各个国家在满意度和应用程序使用时间上似乎存在差异。蓝色国家的用户在应用程序中花费的时间更长，但平均满意度较低。此外，当我们单独观察这三个国家时，我们可能会认为应用程序使用时间和满意度之间的关联确实是正向的。这是否与我们之前的分析相矛盾呢？

辛普森悖论

我们刚才看到的效果被称为辛普森悖论。当数据的相关性在不同群体和各群体内部时不同时，就会出现这种情况。尽管非常违反直觉，但这确实可能发生（正如我们刚才看到的），其原因是混杂变量。让我们用上面的例子来解释。当单独观察每个国家时，我们看到一个正向趋势：在应用程序中花费的时间越多，满意度越高。然而，正如我们已经看到的，各个国家在平均满意度和应用程序使用时间上存在差异。在蓝色国家中，平均满意度较低，但平均应用程序使用时间比橙色或绿色国家高；这与各个国家内部的趋势相反。然而，可能存在另一个导致这种情况的变量。例如，我们可以想象，在蓝色国家中，更多的人经常感到无聊，导致总体上的不满意度较低（因此对我们的应用程序心情较差），但却有更多的时间可以在应用程序上使用。当然，这只是一个可能的解释，还可能有其他许多解释。然而，目前正确的解释并不太重要。对我们来说重要的是要理解各个国家之间存在着系统性差异。

那么，为什么我们在之前的分析中没有发现这一点？我们在执行线性回归时犯了错误吗？是的，因为进行线性回归是错误的，因为线性回归的核心假设之一被违反了：线性回归假设所有数据点是独立采样且来自相同的分布。然而，在我们的例子中并非如此！显然，应用程序使用时间和满意度的分布在不同国家之间是不同的。现在，如果线性回归的假设被违反，那么线性回归就不是适合数据分析的正确工具。

层次模型

现在，我们可以如何以更合适的方式分析我们的数据呢？幸运的是，统计模型可以将线性回归的思想扩展到层次数据上。我们所说的层次数据是指我们采样的数据点嵌套在层次结构中，就像在我们的例子中，我们询问的人嵌套在国家中一样。这些统计模型被称为层次线性模型、多层模型或线性混合效应模型。这些模型通过引入所谓的固定效应和随机效应来考虑群体结构。在一个简单的例子中，我们想要根据一个变量来预测另一个变量（例如，我们想要根据在应用程序中花费的时间来预测满意度），固定效应由一个截距和一个斜率组成，适用于所有群体。到目前为止，这与线性回归是完全一样的。

现在随机效应可以为每个组单独引入一个偏差，从而偏离固定截距。例如，蓝色国家的截距可能会稍低，而绿色国家的截距可能会稍高于固定截距。这可以解释各个国家满意度平均水平的差异。

此外，随机效应可以为每个组引入斜率的偏差。例如，在橙色组中，斜率可能会高于固定斜率（即满意度与时间消耗之间的关联更强），而在绿色国家中，斜率可能会较低。

层次模型的实际应用

让我们通过实际操作来了解发生了什么。我们进行了一项新的分析，但现在我们使用statsmodels的mixedlm函数。我们明确指出我们要通过公式“Satisfaction ~ Time_spent”预测满意度（而不是反过来），并指示我们的数据框中的“Country”列确定了不同的组。此外，参数re_formula=“Time_spent”告诉模型我们要为每个组设置一个单独的斜率。如果没有这个参数，随机效应只会考虑组特定的截距，而不是组特定的斜率。

import statsmodels.formula.api as smfresult = smf.mixedlm("Satisfaction ~ Time_spent", data=df, groups=df["Country"], re_formula="Time_spent").fit()print(result.fe_params)print(result.random_effects)

如果我们打印fixed_effects（fe_params）和random_effects，我们会得到以下值：

Fixed effects  Intercept     2.286638  Time_spent    0.497657Random Effects  {0: Group -0.958805, Time_spent -0.018178,   1: Group 0.155233,  Time_spent 0.274222,   2: Group 0.803572,  Time_spent -0.256044}

那么这意味着什么呢？对于固定效应，我们有一个截距的值和一个变量time_spent的值。然而，对于随机效应，我们有两个值（每个国家0、1、2）：一个是截距（Group），一个是变量的斜率（Time_spent）。如上所述，随机效应描述了每个组相对于平均效应的偏差。对于我们的三个组，我们可以通过将随机效应添加到截距和斜率的固定效应上来构建三个不同的线性方程。

satisfaction_0 = (2.286638 - 0.958805) + (0.497657 - 0.018178) * time_spent = 1.327833 + 0.479479 * time_spentsatisfaction_1 = (2.286638 + 0.155233) + (0.497657 + 0.274222) * time_spent = 2.441871 + 0.771879 * time_spentsatisfaction_2 = (2.286638 + 0.803572) + (0.497657 - 0.256044) * time_spent = 3.090210 + 0.241613 * time_spent

我们可以看到，组0的随机截距为负数（-0.958），组2的随机截距为正数（0.803），因此组0在线性函数中的截距最低（1.327），组2最高（3.090）。换句话说，在国家0中，满意度的起始水平低于国家2。

我们还可以看到组之间的斜率不同。在组1中，斜率最高为0.771，而在组2中只有0.241。这意味着在国家1中，应用程序使用时间增加一小时会导致满意度增加0.771点（平均），而对于国家2，只会增加0.241点满意度。此外，所有的斜率都是正值，这与我们在开始时进行的线性回归的负斜率相矛盾。

现在我们可以为每个国家绘制一条回归线：

现在我们清楚地看到了每个国家的正向趋势和不同的截距（即线在time_spent=0时的位置）。

结论

在上面的例子中，我们看到了一个目光短浅的分析如何容易地导致错误的结论。忽视数据的嵌套结构，即来自不同国家的用户，我们可能很容易在线性回归之后停止并得出结论：在应用程序中花费更多时间与满意度较低相关。只有通过理解我们的数据不符合线性回归的核心假设，因为数据点不是从相同的分布中抽样得到的，我们才被激发进一步分析，揭示出事实上相反的情况：在应用程序中花费更多时间确实与更高的满意度相关。

因此，让我们从这个例子中得出一些要点：

• 在使用统计方法进行分析之前，应该用数据验证其假设。
• 对于嵌套数据，假设所有数据点都是从相同分布中抽样得到的可能并不总是成立。
• 数据整体上的趋势可能与构成该数据的单个组内的趋势不同。这被称为辛普森悖论。
• 多层线性模型是处理嵌套数据结构并避免辛普森悖论导致错误结论的一种方法。

进一步阅读

我们在statsmodels中使用了以下层次模型的实现：

• https://www.statsmodels.org/stable/mixed_linear.html

我使用了以下统计学教材（可惜只有德语版本）。

• Eid, M., Gollwitzer, M., & Schmitt, M. (2017). Statistik und Forschungsmethoden.

关于多层模型的背景信息还可以在这里找到：

• Snijders, T. A. B.; Bosker, R. J. (2011). Multilevel Analysis: an Introduction to Basic and Advanced Multilevel Modeling (2nd ed.). London: Sage. ISBN 9781446254332.

如果您想要重现结果，可以按照以下方式生成数据：

import numpy as npimport pandas as pdgroup_1_x = np.random.uniform(0.5, 1.8, 25)group_1_y = (1 + 0.3 * group_1_x) + np.random.rand(len(group_1_x))#start_2, end_2, step_2 = 0.3, 1.3, 0.04group_2_x = np.random.uniform(0.3, 1.3, 22)group_2_y = (2 + 0.7*group_2_x) + np.random.rand(len(group_2_x))#start_3, end_3, step_3 = 0, 1, 0.04group_3_x = np.random.uniform(0, 1, 32)group_3_y = (2.5 + 0.3*group_3_x) + np.random.rand(len(group_3_x))all_x = np.concatenate([group_1_x, group_2_x, group_3_x])all_y = np.concatenate([group_1_y, group_2_y, group_3_y])df = pd.DataFrame({"满意度": all_y, "花费时间":all_x, "国家":[0]*len(group_1_x) + [1]*len(group_2_x) + [2]*len(group_3_x)})

喜欢这篇文章吗？关注我以获取我的未来文章通知。