解释熊猫中带有复制警告的设置
如果你想知道是什么导致了SettingwithCopyWarning，这就是你要找的地方！
不管你和pandas共事多久，是一天还是一年，迟早你都有可能遇到臭名昭著的SettingWithCopyWarning。在本文中，我将解释是什么导致了这个问题，以及如何正确地解决这个问题。
警告不是错误
在我深入研究技术细节之前，我想强调一下SettingWithCopyWarning是——顾名思义——一个警告，而不是一个错误。所以我们正在执行的代码很可能不会中断并产生最终结果。然而，最终的结果可能不是我们真正想要的。
我想强调这种区别的原因是，当我们看到代码实际上成功返回结果时，我们可能会忽略警告。事实上，结果可能是正确的！最佳实践是格外小心，并真正理解基本原理。通过这种方式，我们通常可以节省大量时间来识别一个不明显的 bug，这是我们一开始就可以避免的。
视图与副本
与SettingWithCopyWarning相关的关键概念是视图和副本。pandas中的一些操作(和numpy一样)将返回原始数据的视图，而其他的副本。
简单地说，视图是链接到原始源的原始对象(DataFrame或Series)的子集，而副本是一个全新的对象。一般来说，当我们对副本进行的操作完成后，副本就会被丢弃。这种区别的结果是，当我们修改视图时，我们也修改了原始对象。副本不会发生这种情况，因为它们没有连接到原始对象。
描述了区别之后，SettingWithCopyWarning实际上是让我们知道我们写的代码可能做了一件事，而实际上我们想做的是另一件事。我将用一个真实的例子来说明这一点。想象一下有一个大的DataFrame。对于某些分析，您过滤(切分)了DataFrame，使其只包含完整数据的子集，例如，来自某个国家的用户。然后，您可能想要修改提取的DataFrame中的一些值，假设一个特性的最大值为 100。这是您可能遇到臭名昭著的警告的典型情况——您只想修改提取的帧，而最终却修改了源数据。你可以很容易地想象这不是你想做的事情，并可能导致以后的潜在问题。
注意:要了解一个框架是否是副本的视图，可以使用一个pandas.DataFrame的内部_is_view和_is_copy方法。第一个函数返回一个布尔值，而第二个函数要么是原来的DataFrame的weakref，要么是None。
警告的常见事件
在这一节中，我回顾了实践中发生SettingWithCopyWarning时最常见的情况。我将使用一个小的定制DataFrame来举例说明这些情况，因为这对于理解逻辑来说已经足够了。
为了准备数据，我运行以下代码:
运行代码打印出小的DataFrame:
为了清楚地理解正在发生的事情，对于下面的每一种情况，我们将从头开始—运行get_data函数的结果。
作为将来的参考，本文是使用pandas版本 1.0.3 编写的。
1.链式分配
为了解释链式赋值的概念，我们将依次讨论构造块。赋值操作(也称为设置操作)简单地设置一个对象的值。我们可以通过创建一个列表来说明这一点:
x_list = [1, 2, 3, 4, 5]
尽管第一个例子是基于列表的，但是同样的原则也适用于数组、Series和DataFrames(我们马上就会看到)。第二种类型的操作称为 get 操作，用于访问和返回对象的值。索引是一种 get 操作，我们可以通过运行返回的x_list[0:3]来索引列表
[1, 2, 3]
最后一个构建块叫做链接，本质上是指链接多个索引操作，比如x_list[0:3][1]，返回 2。
已经描述了所有的单个部分，通过链接赋值我们指的是链接和赋值的组合。是时候提及我们的玩具DataFrame了。首先，我们对DataFrame进行切片，以显示B特性值大于 12 的观察值。
X = get_data()
X[X['B'] > 12]
只有 2 行符合该标准。让我们用 999 替换C特性的值。
X[X['B'] > 12]['C'] = 999
X[X['B'] > 12]['C']
运行上面的行会导致臭名昭著的警告:
SettingWithCopyWarning:
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead
在生成的输出中，我们看到这些值没有被替换！
我们看到警告是因为我们链接了两个索引操作。这些连锁操作在后台独立执行。第一部分是 get 操作，它返回一个包含所有B值大于 12 的行的DataFrame。第二部分是设置的操作，并且是对新建的DataFrame进行由获取的操作。所以我们不是在修改原来的DataFrame！
当我们在一行中使用两个方括号时，这是非常明显的，然而，通过使用loc、iloc或访问列的点方法也会发生同样的情况。例如，运行X.loc[X[‘B’] > 12][‘C’] = 999会给出同样不正确的结果。
为了正确替换DataFrame中的值，我们需要以如下方式使用loc:
X.loc[X['B'] > 12, 'C'] = 999
X[X['B'] > 12]['C']
我们可以看到原来的DataFrame中的值被成功替换。
2.隐藏链接
隐藏链接可能是一个很难调试的问题，因为问题到底出在哪里通常不是很明显。我们将复习一个例子。首先，让我们加载数据，并使用前面案例中的知识，创建一个DataFrame，它是原始数据的子集。我们过滤掉所有具有大于 101 的C特性值的行。
经常发生的是，我们探索并进一步处理新的DataFrame。让我们想象运行几行代码来进一步检查temp对象，比如pandas.DataFrame的shape、describe、plot方法。
我们实际上没有在这里打印输出，因为这不是重要的部分。现在，再运行几行代码后，让我们用 999 替换第一行temp中C特性的值:
temp.loc[2, 'C'] = 999
在这样做的时候，我们遇到了我们的老朋友:
SettingWithCopyWarning:
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead
让我们检查原始的和提取的DataFrames中C的值:
print(f"New DataFrame: {temp.loc[2, 'C']}")
print(f"Original DataFrame: {X.loc[2, 'C']}")# New DataFrame: 999
# Original DataFrame: 102
所以我们真的发生了？警告的原因在于，链式索引可能出现在两行中，而不仅仅出现在一行中。当我们创建新的DataFrame时，我们使用了 get 操作的输出。这可能是原件DataFrame的副本，也可能不是。在我们检查之前无法知道。引用关于链式索引的pandas文档:
除了简单的情况之外，很难预测它[链式索引]将返回一个视图还是一个副本(这取决于数组的内存布局，对此 pandas 不能保证)，…
所以当我们索引temp来分配新值时，我们实际上使用了链式索引。因此，我们可能在修改temp的同时也修改了X。
棘手的部分是，在实际的代码库中，负责隐藏链式分配的两行代码可能被几十行代码分开，这使得识别潜在问题相当困难。
为了解决这个问题，我们可以通过使用copy方法直接指示pandas创建原始DataFrame的副本。下面，我们使用这种方法来避免隐藏链接:
我们看到运行代码生成了正确的输出。
3.一个假阴性的例子
最后，我们回顾一下[3]中提到的假阴性(pandas没有通知我们SettingWithCopyWarning，而它实际上应该通知我们)。当我们在对DataFrame的多个列进行切片时使用链式索引时，就会发生这种情况。我们来看两个简单的案例。
X = get_data()
X.loc[X['A'] > 2, ['A', 'B']]['A'] = 999
X
和
X = get_data()
X[['A', 'B', 'C']]['A'] = 999
X
两者都产生了下面的结果，但没有真正显示出SettingWithCopyWarning。
从上面的图片中我们可以看到，这些值并没有按照我们想要的那样被修改。由于在索引操作中包含多个列，因此没有显示警告。我们可以很容易地验证在单个列的情况下不会发生这种情况——运行X[[‘C’]][‘C’] = 999
会产生警告并且不会修改X。
关于警告来源的更深层次的背景
我们可以说pandas从numpy继承了视图和副本的概念。在引擎盖下，pandas使用numpy进行高效的数据存储和操作。在numpy中，视图和副本遵循一组特定的规则，并以可预测的方式返回(更多信息见【5】)。那么为什么pandas不是这样呢？问题在于numpy数组被限制为单一数据类型。正如我们所知，pandas的情况并非如此。
实际上，在多数据类型DataFrame上的索引( get 操作)将总是返回一个帧的副本。对单一类型框架的相同操作几乎总是返回基于单一numpy数组的视图，这是解决问题的最有效方式。然而，正如我们在文档的引用中已经看到的，返回视图依赖于对象的内存布局，不幸的是不能保证。
综上所述，pandas尽最大努力将其通用的索引方法(由于这种方法，它非常流行，并且基本上是用 Python 进行数据科学研究的先决条件)和使用底层numpy数组的效率结合起来。这导致了一些小麻烦，然而，这种权衡绝对是值得的，并且通过正确理解pandas在引擎盖下如何工作，这些问题可以被克服。
结论
在这篇文章中，我解释了 pandas 中副本和视图之间的区别，以及它们与臭名昭著的SettingWithCopyWarning之间的关系。主要的想法缩小到知道什么是链式索引以及如何成功地避免它。一般规则是:
- 如果要改变原来的
DataFrame，使用单赋值。 - 如果您想要制作一个
DataFrame的副本，请使用copy方法显式地这样做。
遵循这两条规则可以节省您调试一些奇怪情况的大量时间，尤其是在冗长的代码库中。
同样值得一提的是，SettingWithCopyWarning仅在我们使用 set 操作(赋值)时发生。然而，最好也避免对 get 操作进行链式索引。这是因为链式操作通常比较慢，如果您后来决定将 set 操作添加到代码中，可能会导致问题。
您可以在我的 GitHub 上找到本文使用的代码。一如既往，我们欢迎任何建设性的反馈。你可以在推特上或者评论里联系我。
参考
[1]https://pandas . pydata . org/pandas-docs/stable/user _ guide/indexing . html
https://github.com/pandas-dev/pandas/issues/9767
[4]https://www . practical data science . org/html/views _ and _ copies _ in _ pandas . html
[5]https://scipy-cookbook . readthedocs . io/items/viewsvscopies . html
有影响力的解释
活动讲座
谢尔登·费尔南德斯和迈克尔·圣·朱尔斯| TMLS2019c
关于扬声器:
DarwinAI 首席执行官 Sheldon Fernandez 是一位经验丰富的高管，也是技术和企业社区中备受尊敬的思想领袖。在他的职业生涯中，他将人工智能等新兴技术应用于企业客户的实际场景中。谢尔登也是一位有成就的作家和演说家。他在许多场合的许多会议上发表过演讲，包括湾区著名的智库 Singularity University，并撰写了许多主题的技术书籍和文章，包括人工智能和计算创造力。
— Michael St. Jules 是 DarwinAI 的高级研究开发人员，自 2018 年初以来一直在该公司工作。他于 2014 年和 2016 年分别获得卡尔顿大学和渥太华大学数学专业的 BMath 和 M.Sc 学位，主要研究数学分析、逻辑和计算机科学，硕士论文和出版物为量子密码学。然后他在滑铁卢大学的计算数学 MMath 中重点研究了机器学习和深度学习，于 2017 年毕业，并在加入 DarwinAI 之前担任 Alexander Wong 博士的研究助理。
关于演讲:
围绕人工智能的普遍进展引发了人们对可解释的人工智能(XAI)的浓厚兴趣，其目标是产生由机器学习算法做出的可解释的决策。鉴于深度神经网络的复杂性和“黑箱”性质，特别感兴趣的是对深度神经网络如何做出决策的解释。
鉴于该领域的起步阶段，对可解释方法的性能评估的探索有限，大多数评估集中在对当前方法的主观和视觉解释上。在本次演讲中，演讲者介绍了两种量化性能指标，用于通过一种新颖的决策影响分析来量化深度神经网络上可解释性方法的性能:
- 影响得分，评估具有强烈信心降低影响或决策改变影响的关键因素的百分比；和
- 影响范围，评估输入中受到不利影响的因素的百分比范围。我们进一步考虑使用这种方法与众多最新的可解释性方法进行综合分析。
订阅我们的月刊,直接在你的邮箱✉️中接收我们最好的文章、视频和播客
用 SHAP 和莱姆解释你的机器学习模型！
帮助您揭开一些人可能认为您的机器学习模型是“黑箱”的神秘面纱
大家好！欢迎再次回到另一个数据科学快速技巧。这篇特别的文章对我来说是最有趣的，不仅因为这是我们迄今为止处理的最复杂的主题，而且也是我刚刚花了几个小时自学的一个主题。当然，还有什么比想出如何把它教给大众更好的学习方法呢？
在开始之前，我已经把这篇文章中展示的所有作品上传到了一个奇异的 Jupyter 笔记本上。如果你想更深入地了解，你可以在我的个人 GitHub 找到它。
因此，尽管这是一个幕后的非常复杂的话题，但我还是打算尽可能地为尽可能多的观众淡化这个话题。尽管这最终是一篇为数据科学从业者设计的帖子，但我认为对于任何业务人员来说，理解他们为什么应该关注这个话题也同样重要。
在进入如何计算/可视化这些值之前，让我们先对为什么我们会关心这个话题建立一些直觉。
为什么你应该关心 ML 的可解释性
如果你对机器学习有所了解，你会知道你将一些数据放入一个预测模型，它会在另一端产生一个输出预测。但是让我们诚实地说:你知道在那个模型里面到底发生了什么吗？请记住，我们正在谈论的是在引擎盖下进行大量复杂数学运算的算法。实际上，这可能就是你的 ML 模型现在的样子:
你真的不知道那个黑盒里发生了什么。数学完全掩盖了这些特性在模型中的重要性。那么，如果特性 1 和 2 承担了所有的重量，而特性 3 和 4 没有给模型增加任何价值呢？您更希望看到的是类似这样的东西:
好的，在我们的例子中，我们显然使用了外行人的术语。但这一点很重要，因为如果你的模型依赖一两个特征来做所有的预测，这可能表明有问题。这也可能表明你正在使用某种可能被认为是不道德的功能。例如，今天许多公司会避免在他们的模型中使用性别或种族之类的东西，因为它们往往会在模型中产生许多不希望的偏见。
通过一个具体的例子来理解这些概念总是更有意义，所以我们将回到我们在最近几篇文章中使用的数据集来具体演示 SHAP 和莱姆:泰坦尼克号数据集。
用泰坦尼克号数据集建模
好吧，所以在我们使用泰坦尼克号数据集的最后几个帖子中，我们创建了非常差的模型，因为我们的重点更多的是学习新技能。为了有效地演示 SHAP 和莱姆，我们真的需要一个模型来显示某种程度的“尝试”(？)这样我们就可以看到 ML 的可解释性是如何发挥作用的。现在，我不打算在这篇文章中花时间向您展示我的数据清理和建模部分的所有代码(您肯定可以在我的 GitHub 中找到这里，但是我将快速总结产生最终模型的特性。
- PClass :记录这个人上船的船票等级。(1 档、2 档或 3 档。)
- 年龄箱:根据人的年龄，他们被分成一般人群。这些群体包括儿童、青少年、年轻人、成年人、老年人和未知人群。
- 性别(Gender) :注明该人是男是女。
- 船舱段:如果这个人有船舱，我就拉出他们可能在哪个舱。这些可以是 A 舱、B 舱、C 舱，也可以是无舱。
- 登船:记下此人最初出发的三个地点之一。
- SibSp(兄弟姐妹/配偶):记录此人在船上有多少兄弟姐妹和/或配偶。
- Parch (Parent / Child) :记录这个人在船上有多少父母/孩子。
对于实际的模型本身，我选择使用来自 Scikit-Learn 的一个简单的随机森林分类器。将数据分为训练集和验证集，我得到了验证集的以下指标:
这些指标并不是很好，但是对我们的学习来说很好。至少我们将在下一节看到这些特性有多重要。让我们继续前进！
用 SHAP 和石灰解释
在这一点上，我们目前还不知道我们的随机森林模型是如何使用这些特征进行预测的。我猜性别和年龄是影响因素，但我们还不太清楚。
这是聪明人基于博弈论开发复杂算法来生成这些东西的地方，以更好地解释我们当前的黑匣子。当然，我指的是 SHAP 和酸橙。现在，我要 100%诚实:我不太熟悉这些算法是如何在引擎盖下工作的，所以我甚至不打算在这篇文章中尝试这样做。相反，我会给你介绍这篇关于 SHAP 的文章和另一篇关于石灰的文章。
就像 Scikit-Learn 抽象出我们的随机森林分类器的底层算法一样，我们将使用一些简洁的 Python 库来抽象出 SHAP 和莱姆的内部工作原理。要使用它们，您所需要做的就是对两者进行简单的 pip 安装:
pip install shappip install lime
在高层次上，这两种工作方式都是将您的训练数据和模型交给“解释者”，然后您可以将任何观察结果传递给“解释者”，它会告诉您该模型的特性重要性。这听起来像是天书，但是我们会用泰坦尼克号的例子来说明这一点。
首先，我们需要从我们的验证集中选出两个人，我们知道他们分别没有幸存和幸存。这是因为我们将通过 SHAP 和莱姆的解释来了解这些人的哪些特征对我们的模型的生存能力有最大的影响。(还是那句话，你可以在我的笔记本里看到这个作品。)
通过 SHAP 解释
让我们从 SHAP 开始。这里的语法非常简单。我们将首先实例化 SHAP 解释器对象，使我们的随机森林分类器(rfc)适合该对象，并插入每个人以生成他们可解释的 SHAP 值。下面的代码向您展示了如何为第一个人做这件事。要对人员 2 进行同样的操作，只需交换适当的值。
# Importing the SHAP library
import shap# Instantiating the SHAP explainer using our trained RFC model
shap_explainer = shap.TreeExplainer(rfc)
shap.initjs()# Getting SHAP values for person 1
person_1_shap_values = shap_explainer.shap_values(person_1)# Visualizing plot of expected survivability of person 1
shap.force_plot(shap_explainer.expected_value[1], person_1_shap_values[1], person_1)# Visualizing impact of each feature for person 1
shap.summary_plot(person_1_shap_values, person_1)
这真的就是全部了！让我们在下面的截图中看看我们刚才所做的实际输出。首先，一号人物。
根据上图，我们的模型预测第一个人有 94%的机会幸存，这是正确的。阅读起来有点困难，但看看第二个视觉，看起来最有影响的因素包括这个人的性别，这个人是否有三等票，以及这个人是否是个孩子。在这种情况下，第一个人是一个女孩，所以这两个因素促使我们的模型说是的，这个人可能活了下来。
现在让我们看看第二个人。
我们的模型预测这个人可能没有活下来，不幸的是这是正确的。再看看这里的影响因素，看起来第二个人是一个没有小屋的成年男性。我不是历史爱好者，但我猜想妇女和儿童首先被救出船外，所以对我来说，第一个人(一个女孩)幸存而第二个人(一个成年男子)没有生还是很有意义的。
好吧，让我们继续我们的石灰解释！
通过石灰的可解释性
使用 LIME 的语法有点不同，但是从概念上讲，SHAP 和 LIME 做了很多相似的事情。
# Importing LIME
import lime.lime_tabular# Defining our LIME explainer
lime_explainer = lime.lime_tabular.LimeTabularExplainer(X_train.values, mode = 'classification', feature_names = X_train.columns, class_names = ['Did Not Survive', 'Survived'])# Defining a quick function that can be used to explain the instance passed
predict_rfc_prob = lambda x: rfc.predict_proba(x).astype(float)# Viewing LIME explainability for person 1
person_1_lime = lime_explainer.explain_instance(person_1.iloc[0].values, predict_rfc_prob, num_features = 10)
person_1_lime.show_in_notebook()
让我们来看看 LIME 为每个人展示了什么。
我在这里有点偏见，但我认为这些石灰值的用户界面更容易阅读 SHAP。我们看到每个人是否幸存的预测概率是一样的。在中间，我们看到了 10 大特征的影响力。与 SHAP 相比，莱姆在解释能力上有一点点不同，但它们大体上是相同的。我们再次看到，性别是一个巨大的影响因素，无论这个人是否是个孩子。在右边，它们很好地显示了这个特定观察的精确值。
那么 SHAP 和莱姆哪个更好呢？我会说都不是。虽然两者都可以很好地相互补充，但我认为向您的业务用户展示这两者是很重要的，这样可以让他们更好地了解哪些功能最容易解释。我认为没有人能肯定地说 SHAP 和莱姆哪个更准确，所以两者都要展示出来。无论如何，语法很容易理解！
这就结束了另一个帖子！希望你们都喜欢这个。我要说这是我迄今为止最喜欢学习和写的，我期待着在我自己的工作中应用这一点。在下一篇文章中与大家见面！
解释性建模和统计因果推理教程
通过对新冠肺炎死亡危险因素的案例研究
来源:马库斯·斯皮斯克
介绍
也许最近的疫情和相关的封锁的一个积极的副作用是，我们可以实现被搁置的项目。有一段时间，我想写一篇关于统计模型解释的文章，我需要的只是时间、动机和数据。这个 Kaggle 数据集由 2019 年 12 月至 2020 年 3 月期间在武汉市采样的 1085 例新冠肺炎病例组成，它是首批发布的包含患者级别记录的数据集之一。这篇 TDS 文章很好地概述了数据及其字段。拥有患者级别的记录意味着有足够的粒度来得出关于一个人的个人属性如何与暴露于新冠肺炎病毒的死亡风险相关的医学见解。
本解释性建模指南要求对以下主题有中级理解:
死亡率风险分析需要一个带有死亡事件记录的响应变量。因此，我们将尝试将这种反应表达为预测变量的函数，这些变量与新冠肺炎对人类健康的影响直接或间接相关。解释模型的目的是解释而不是预测死亡的结果，其中解释的目的是应用统计推断以便:
- 识别对死亡事件有统计显著影响的预测变量
- 估计重要预测因素的影响程度
- 量化我们估计的不确定性
换句话说，我们感兴趣的是一个用于推理而不是预测的模型，或者更准确地说，一个更倾向于因果推理而不是*预测推理的模型。*理解两个之间的差异至关重要，因为它决定了我们选择建模方法的基本原理。解释性分析需要一个不太灵活(高偏差)但更加可解释的模型，它使用概率推理来近似一个假设的数据生成过程。在执行了模型选择过程之后，得到的“最佳”模型可以包括对响应具有显著影响的预测器，尽管不一定提高预测准确度。
从表面上看，一个不错的模型选择是二项式 logit 模型，通常称为二元逻辑回归。考虑到模型选择后缺乏用于推理的工具以及我们验证因果假设的意图，我们将使用领域知识来精选假设对死亡事件有直接或间接影响的独立变量；然后，我们将通过将我们的数据拟合到我们的模型假设的数据生成过程中来量化它们的影响。本文的目的是实现一个解释性建模的管道，它展示了如何使用一个小而稀疏的数据集来获得关于新冠肺炎的医学见解。为此，我有选择地将理论概念与 R 编程语言中的代码片段结合起来。为了达到理想的结果，采取以下步骤:
- 广义线性模型的定义
- 数据插补与数据删除
- 渐近准则下的模型选择
- 基于 bootstrap 推理的模型验证
- 模型解释
本作品中使用的代码库可以在这里找到。
让我们从加载必要的库和数据开始。
*# Load required libraries:*
library(mice)
library(VIM)
library(jtools)
library(ggplot2)
library(gridExtra)
options(repr.plot.width = 14, repr.plot.height = 8)*# Load COVID-19 data:*
dat0 = read.csv('./data/COVID19_line_list_data.csv',
header = T, sep = ',')
print(names(dat0))
在预测建模场景中，我们可以在这一点上选择我们想要的任何数量的变量，然后让模型选择算法来决定我们应该考虑哪些变量。然而，在解释性建模场景中，算法模型选择之后的推断并不是一个可行的选项。一个激烈的变量选择过程预计会偏向系数 p 值，并剥夺我们使用模型的渐近性质。有各种方法来克服这种现象，但是，为了解释因果关系，我们必须依靠领域知识来隔离我们认为有影响的变量。例如，一个显而易见的变量选择是与患者的年龄、性别、原籍国、症状开始日期、住院日期以及他/她是否是武汉本地人相关的变量。诸如“症状”、“康复”、“医院 _ 就诊 _ 日期”和“暴露 _ 开始”等变量看似相关，但实在太稀疏，没有任何用处。
dat1 = dat0[c("location", "country", "gender", "age", "reporting.date", "symptom_onset", "hosp_visit_date", "visiting.Wuhan", "from.Wuhan", "death", "recovered", "symptom")]
print(head(dat1, 5))
隔离最相关的预测值后，我们的新数据集应该如下所示:
*# Binarize response:*
dat1$y = 0
dat1$y[dat1$death==1] = 1
dat1$y = as.factor(dat1$y)
我们有一个代表死亡事件的二元响应变量，我们准备更深入地研究建模过程。
模型假设
在不打算进入太多细节的情况下，我认为为了理解模型可解释性是如何获得的，经历一些基本步骤是很重要的。
回想一下，在二项式 logit 模型中，响应变量是遵循二项式或伯努利分布的随机向量:
其中P(Y =1)= P，对于所有观测值 i= 1 ，…，N.
这是与线性回归的一个重要区别，在线性回归中，随机因素来自协变量，而不是响应。我们的模型做出了其他几个假设，分析的第一步应该是确保我们的场景满足这些假设。
关于响应分布的第一个假设已经确定，但是可以扩展为假设其期望值E[Y =1】= NP，对于所有的 i 、可以与协变量的线性组合相关联:
这相当于:
对于观测值 i =1、…、 N ，解释变量 j =1、…、 M ，链接函数 g 。其次，假设解释变量 X 相互独立，不相关。第三，观测值也被假定为相互独立和不相关的。一些进一步的假设与最大似然估计对数据大小的限制有关。众所周知，当参数的数量增加超过一定限度时，最大似然估计往往会过度拟合，但幸运的是，我们的数据并没有接近那个危险区域。
一个更普遍的假设是，我们的数据是随机抽样的结果。这在预测的情况下尤其重要，因为响应变量(样本类别平衡)的分布对截距项的估计有直接影响，因此也对模型预测有直接影响。好的一面是，斜率系数不受影响。换句话说，逻辑回归中的非代表性类别平衡对预测有直接影响，但对解释没有影响。因为我们的目标是后者，所以我们可以不做平衡样本假设。
定义 GLM
我们现在可以将我们的模型定义为二项式 GLM。回想一下，任何 GLM 都由三个部分定义:
- 随机成分
- 系统的组成部分
- 链接功能
假设随机变量 Y 将新冠肺炎案例的二元结果映射为伯努利试验，使得:
其中 p 是观察到单个患者死亡结果 P ( Y =1)的概率， N 是具有 Y 观察值的患者数量， n 是每个患者的试验次数(在我们的例子中， n =1，对于所有 n )。这是我们的随机成分。
设 η 为 Y ，E[Y]=η的期望值，假设存在一个函数 g (。)使得 g ( η )是协变量的线性组合:
对于 X = x ，对于所有的 j ，上式的 RHS 就是我们的系统分量。
我们选择 g (。)= logit (。)作为反函数的环节，逻辑(。)，可以把我们模型的输出挤在[0，1]区间内，这也恰好是E[Y]=η= p，for n =1 的期望范围。因此，通过将η=NP=p插入我们的链接函数 g (。)，我们得到:
从中我们得出模型的最终结构:
为了简单起见，我们省略了索引:
随机分量现在可以改写为:
注意:最后一个表达式解释了为什么逻辑回归是字面上的回归和而不是分类。
数据准备
这一步包括检查数据集中的方差和稀疏性，以及生成新的预测值。定义了 GLM 并检查了其假设后，我们可以继续检查数据中的类别平衡:
table(dat1$y)
该数据集与我们对新冠肺炎死亡率的先验知识一致，这使我们的分析进入了罕见事件检测的领域。利用新冠肺炎的领域知识，我们可以立即考虑将’年龄’和’性别’作为相关的预测因子，但我们也可以尝试创建新的预测因子。从第一次症状出现到患者住院的时间可能会提供一些信息。我们可以很容易地通过计算住院日期和症状发作日期之间的天数来创建这个预测值。
*# Cast 'Date' types:*
dat1$hosp_visit_date = as.character(dat1$hosp_visit_date)
dat1$symptom_onset = as.character(dat1$symptom_onset)
dat1$hosp_visit_date = as.Date(dat1$hosp_visit_date, format="%m/%d/%y")
dat1$symptom_onset = as.Date(dat1$symptom_onset, format="%m/%d/%y")
dat1$from.Wuhan = as.factor(dat1$from.Wuhan)*# Create 'days before hospitalization' variable:*
dat1$days_before_hosp = as.numeric(as.Date(dat1$hosp_visit_date, format="%m/%d/%y") - as.Date(dat1$symptom_onset, format="%m/%d/%y"))
dat1$days_before_hosp[dat1$days_before_hosp < 0] = NA
因此，在稍后的模型选择过程中，可以将计算症状发作和住院之间天数的变量添加到我们的模型中。数据的当前状态似乎包括相当多的缺失值(显示为’ NA ')，因此可视化稀疏性前景将是有见地的。
# Plot missing values:
aggr(dat1, col = c('green','red'), numbers = TRUE, sortVars = TRUE,
labels = names(dat1), cex.axis = .5, gap = 2,
ylab = c("Proportion in variable","Proportion in dataset"))
每列缺失值的比例
上表显示了每个变量中缺失值的比例。下面的柱状图展示了同样的观点(左图)。
缺失的数据环境
右侧的网格显示了单个预测的缺失值占数据集中所有值的比例。例如，矩阵的最后一行，全绿色，告诉我们 39%的数据集案例在所有变量中都没有缺失值。前一行告诉我们，前三个变量(出现的顺序)中缺少 20%的事例，前一行显示前五个变量中缺少 15%的值，依此类推。
很明显，有很多缺失的数据，而且它们分散在各处，所以在剔除和插补之间的选择并不明显。在试图回答这个问题之前，我们应该尝试使用少量但大量的预测因子来拟合基线模型。我们希望在接近原始数据集的数据集上拟合一个模型，但我们将执行缺失值移除，而不是插补。具有稀疏性或低方差的变量(例如’国家’、‘医院 _ 访问 _ 日期’、‘症状 _ 发作’)应该被排除，并且为了最小化数据移除，包含在基线模型中的变量应该具有尽可能少的缺失值。
模型 0(基线模型)
模型选择过程将包括拟合几个候选模型，直到我们找到最接近“真实”模型的一个。在每次比较时，将根据特定的渐近标准(待解释)对候选模型进行评估。
我们将从一个非常简单的模型开始，该模型根据变量’年龄和’性别’回归死亡率。
将死亡率风险表示为年龄和性别的函数，这使得我们可以将等式[1]中的 GLM 改写为:
# Loading custom functions used in the analysis:
source('../input/glmlib/myglm-lib.r') #modify path as required
我们使用 R 的标准’ glm ‘函数拟合模型，并借助’ jtools ‘软件包和一些自定义函数检查拟合优度。我们使用包’ jtools ‘中的’ summ ‘函数来获得模型诊断的良好摘要。应该注意的是，R 的内置’ glm ‘函数默认情况下会删除所有带有’ NA’ 值的行，所以我们在摘要中被告知，在我们总共 1082 个观察值中，只有 825 个用于训练。我们还看到关于拟合优度统计和标准的报告，如 AIC、BIC 、偏差，它们可用于特定条件下的模型选择。
#---------------------------------------------
# Model 0: gender + age
#---------------------------------------------# Train model:
gfit0 = glm(y~gender+age, data = dat1, binomial)# Summarize and inspect:
summ(gfit0, confint = TRUE, pvals = TRUE, digits = 3)
让我们为这个输出提供一些背景。
回想一下，偏差是普通最小二乘回归中残差平方和的推广，因此它将当前模型与饱和 模型的距离表示为它们各自对数似然 l (。).因此，当前模型的(缩放)偏差由下式获得:
其遵循渐近分布:
自然，偏差总是正的，与拟合优度成反比，用 D = 0 表示完美拟合。零偏差 D0 是仅截距(零)模型的偏差——最差的可能模型:
—通过设置获得:
代入上面的公式。 Ds 是饱和模型的偏差——可能的最佳模型——对于 n 样本量具有 n 个参数的模型，其似然等于 1，偏差等于 0。
偏差可以用来评估任何两个模型是否相等。这可以通过测试较大模型的所有额外系数都等于零的零假设来实现。对于任何给定的嵌套模型对 m1 和 m2 ，可以通过统计进行比较:
换句话说，对于两个相等的嵌套模型，它们相应偏差的差应该遵循一个 χ 分布，其中 k 自由度等于它们相应自由度的差。因此，这意味着我们也可以测试任意模型在统计上是否等于零模型(非常差的拟合)或饱和模型(非常好的拟合)。这些比较的代码如下所示:
# Comparing current model to the null model (LTR):
pval_m0 = 1 - pchisq(gfit0$null.deviance - gfit0$deviance, gfit0$df.null - gfit0$df.residual)
round(pval_m0, 3)
上述比较检验了零假设:
即当前模型是否等于仅拦截模型。自然，我们希望这个空值被拒绝，以使当前模型不被认为是无用的，因此我们期望 p 值低于可接受的决策阈值(通常为 0.05)。
# Comparing current model to the saturated model (Deviance test):
pval_ms = 1 - pchisq(gfit0$deviance, df = gfit0$df.residual)
round(pval_ms, 3)
这里，我们测试空值:
即当前模型在统计上是否等于饱和模型。未能拒绝零是期望的结果，尽管它不能作为零接受的证据。
我们可以通过观察它们的 p 值来评估模型系数的统计显著性。这些天有很多关于 p 值可靠性的讨论，所以我认为应该解释一下这些 p 值是如何获得的，它们真正的含义是什么。回想一下，GLM 的推论是渐近的，拟合系数的真实分布是未知的。MLE 的渐近性质为我们提供了模型系数的高斯分布:
其中 I ( β )表示费希尔信息矩阵。求解 RHS 中的常态 (0，1)，我们导出 j =0，…， M 的 Wald 统计量:
我们用它来计算置信区间和应用假设检验(就像我们在线性回归中使用 t 统计一样)。因此，我们可以利用下式计算系数的 100(1α)%置信区间:
自然，我们也可以使用系数的渐近分布，通过检验零假设来进行关于其显著性的正式假设检验:
该方法使用 Wald 统计量，通常被称为 Wald 检验。
麦克法登的伪 R 是对解释变量的度量，定义为:
其中，当前模型的对数似然比 1 减去仅截距模型的对数似然比。该等式改编自线性回归 1SSE/SST中的 R 公式，但是，对伪 R 的解释是模型与最佳拟合的接近度，而不是模型解释的方差百分比。这里的一个经验法则是，0.2-0.4 之间的值通常表示“好”适合(范围从好到非常好)。
在经历了拟合优度的一些方面之后，我们意识到我们的基线模型客观上并没有那么差。我们得到了总体拟合的显著 p 值以及我们的一个系数(使用 Fisherian 95%置信水平)，并且我们得到了 0.20 的相当不错的伪R。此外，我们可以使用软件包’ jtools ‘中的函数’ plot_summs '来绘制模型系数的渐近分布，并通过注意β2(性别的系数)的尾部如何包含零值且不拒绝空值来可视化 Wald 测试的结果:
# Asymptotic distribution of model coefficients (model 0)
plot_summs(gfit0, scale = TRUE, plot.distributions = TRUE)
基线模型系数的渐近分布
上面的密度图可以帮助我们理解渐近推断是如何工作的。渐近置信区间和 p 值是通过估计该理论密度图上的积分来计算的，并且在适当的条件下，渐近推断的输出应该与任何客观非渐近方法的输出相匹配，例如(frequentist)蒙特卡罗、bootstrap 或平坦先验贝叶斯估计。
使用该模型的预测能力来可视化年龄对男性和女性死亡风险的影响将是有趣的。死亡率风险和年龄之间函数关系的重建如下所示:
# Mortality risk as a function of Age:
x1 = data.frame(gender="male", age=1:90)
plot(predict(gfit0, x1, type = "response"),
type = 'l', xlab = "Age", ylab = "Mortality risk")
x2 = data.frame(gender="female", age=1:90)
lines(predict(gfit0, x2, type = "response"), col="red")
作为年龄函数的预测死亡风险(黑色:男性，红色:女性)
这个模型可以作为基线，所以我们可以从这里开始构建。显而易见的下一步是在模型中加入更多的预测因子，看看拟合度是否会提高。如前所述，该数据集的瓶颈在于我们有太多的缺失值分散在我们的预测器中。如果我们继续添加变量并删除缺少值的行，我们的数据集注定会危险地减少。例如，我们在基线中使用的变量’年龄和’性别’，已经花费了我们 260 次观察。变量’ days_before_hosp '直观上是一个有趣的预测因子，但如果我们将它添加到基线模型中并应用逐行移除，我们数据集中的观察总数将减少到其原始大小的一半。从变量’也是如此。武汉’，表示被抽样的个体是否是武汉本地人(稍后将详细介绍该预测因子的含义)。因此，如果我们排除数据集中至少有一个缺失值的所有行，我们将得到一个更小但完整的数据集。另一方面，如果相信插补算法来“猜测”我们所有的缺失值，我们最终会发现我们的大部分数据是模拟观察的结果。鉴于目前的事态，如果我们要用更多的解释变量来丰富基线模型，我们面临两种选择:
- 数据删除:如果我们假设我们的数据是“随机缺失”(MAR)或“完全随机缺失”(MCAR)，那么应该考虑删除案例。然而，如果这个假设是错误的，那么我们就有可能在最终的精简数据集里引入偏差。
- 数据插补:假设数据“非随机缺失”(MNAR)，有多种插补方法可供选择。
这两种方法都不一定是错误的，我们也没有办法事先知道哪种方法最适合我们的情况，除非我们清楚地了解为什么数据会丢失。由于我们在这种情况下没有这种先验知识，明智的做法是使用两种方法，并生成简化数据集和估算数据集；然后将安装在前者上的模型与安装在后者上的模型进行比较。这种比较不是直截了当的，因为使用不同的数据集(就行输入而言)会导致模型无法通过基于可能性的标准进行比较，如 AIC/BIC、偏差和麦克法登的伪 R 。更好方法的选择将取决于“随机缺失”假设的验证。如果数据删除的结果与数据插补的结果相似，我们可以有把握地假设缺失的数据实际上是“随机缺失的”，在这种情况下，我们将倾向于插补的数据而不是减少的数据。
模型 1(估算数据集)
如前所述，可以添加到基线模式的两个有趣的变量是:
- 'days_before_hosp ':从首次出现新冠肺炎症状到患者随后住院(之前创建)之间的天数
- ‘从。武汉’:表示患者是否来自武汉的因素。同样，后一个变量对死亡率有显著影响有点违背直觉，但在这个过程中会有一个合理的解释。
应该注意的是，我们有意避免自动化这个选择过程，以尽量减少模型比较。这种做法的目的是避免陷入逐步选择的领域，逐步选择是一种广泛使用的技术，如今被明智地视为一种统计学罪过。
我们的下一个候选模型将具有以下嵌套形式:
此时，我们必须创建估算数据集。根据数据中缺失值的性质、数量和分布，有大量的插补方法可供选择。在温和的“随机缺失”假设下，我们决定我们场景的最佳方法是通过链式方程的多重插补，通常称为 MICE 算法。使用“mice”库，我们可以模拟所有缺失数据的值，如下所示:
# Create imputed dataset using 'mice':
dat1_obj <- mice(dat1[c("gender", "age", "days_before_hosp", "from.Wuhan", "country")], m=5, maxit = 50, method = 'pmm', seed = 400)
dat1_imp = complete(dat1_obj)
dat1_imp$y = dat1$y
现在我们在’ dat1_imp ‘中有了估算数据集。接下来，我们将在新估算的数据集上重新调整基线模型和两个新的候选模型，并总结它们的诊断结果。我们正在重新调整我们的基线模型，因为多重插补还填充了之前由’ glm '函数自动删除的关于性别和年龄的缺失数据。我们并不期望这会对基线模型的诊断造成巨大的变化，但为了使所有模型相互之间具有可比性，有必要将其重新调整到与新的候选模型完全相同的数据集。
# Train models on imputed data:
gfit1a = glm(y~gender+age, data = dat1_imp, binomial)
gfit1b = glm(y~gender+age+days_before_hosp, data = dat1_imp, binomial)
gfit1c = glm(y~gender+age+days_before_hosp+from.Wuhan, data = dat1_imp, binomial)# Compare fits:
export_summs(gfit1a, gfit1b, gfit1c, scale = F, error_format = "[{conf.low}, {conf.high}]", digits = 3, model.names = c("model 1a", "model 1b", "model 1c"))
plot_summs(gfit1a, gfit1b, gfit1c, scale = TRUE, plot.distributions = F, inner_ci_level = .95, model.names = c("model 1a", "model 1b", "model 1c"))
比较数据插补后模型系数的渐近分布(基线，住院前天数，来自武汉)
因此，现在我们有三个嵌套模型适用于相同的输入情况，我们可以使用基于可能性的标准直接进行比较。正如预期的那样，在所有方面都有明显的改善，因为我们看到 AIC 和 BIC 值都显著下降，而伪 R 在第三个模型中几乎翻了一番。可以有把握地得出结论，带有两个新变量的模型是我们应该关注的。
95%置信区间图可以被认为是系数渐近分布的“全景”视图，它允许我们评估估计的一致性。两个置信区间是否重叠让我们知道两个模型之间相同系数的估计值是否显著不同。这有助于我们检测多重共线性的存在，这可能是由于向模型中添加了新的变量而导致的。
在这一点上，有人可能想知道，使用基于可能性的度量来比较模型的最佳方式是什么。根据经验，Akaike 信息标准(AIC) 已被描述为留一交叉验证的渐近等价物，因此它是用于预测目的的更好指标。另一方面，贝叶斯信息标准(BIC)已经被在之前引用作为在一组候选模型中找到真模型的更好标准，这是我们为了解释的目的所需要的。请注意，“真实模型”是具有正确的函数形式和正确的回归变量集的模型，这些回归变量可以最好地解释响应，但它不一定是进行预测的最佳模型(解释响应所需的一些变量可能会影响估计的精度，而不会增加预测值)。考虑到这一点，我们将把 BIC 作为选择车型的标准。
我们继续将相同的候选模型拟合到已经经历了数据移除的数据集。
模型 2(精简数据集)
缩减的数据集将通过列表式删除来创建，即移除具有一个或多个缺失值的输入案例(观测值)。我们使用以下代码创建以下内容:
# Create reduced dataset for model 2:
# (filter out rows with NA values)
not_na_mask = !is.na(dat1$y) & !is.na(dat1$days_before_hosp) & !is.na(dat1$age) & !is.na(dat1$gender) & dat1$days_before_hosp >= 0
dat2 = dat1[not_na_mask, ]
table(dat2$y)
我们可以看到观察总数减少了 40%。现在让我们看看拟合模型的诊断结果会是什么样子。
#------------------------------------------------------------------
# Model 2: gender + age + time to hospitalization (on reduced
# dataset)
#------------------------------------------------------------------# Train model:
gfit2a = glm(y~gender+age, data = dat2, binomial)
gfit2b = glm(y~gender+age+days_before_hosp, data = dat2, binomial)
gfit2c = glm(y~gender+age+days_before_hosp+from.Wuhan, data = dat2, binomial)# Merge model summaries:
export_summs(gfit2a, gfit2b, gfit2c, scale = F, error_format = "[{conf.low}, {conf.high}]", digits = 3, model.names = c("model 2a", "model 2b", "model 2c"))# Compare asymptotic distributions of coefficients:
plot_summs(gfit2a, gfit2b, gfit2c, scale = TRUE, plot.distributions = F, inner_ci_level = .95, model.names = c("model 2a", "model 2b", "model 2c"))
比较数据删除后模型系数的渐近分布(基线，住院前天数，来自武汉)
考虑到我们如何将数据集缩减到原来的一半，这个模型的拟合度出奇的好。渐近分布的比较表明，在缩减数据集上拟合的三个模型的诊断与在估算数据集上拟合的诊断一样稳定。这是一个积极的结果，表明缺失数据实际上是如假设的那样随机缺失的，并且插补算法在填充缺失值方面做得很好。在下一节中，我们将进一步加强这一假设，比较两个数据集的最佳候选模型(就 BIC 而言)。
比较估算数据集和缩减数据集
# Merge summaries of two models (imputed vs. reduced)
export_summs(gfit1c, gfit2c, scale = TRUE, error_format = "[{conf.low}, {conf.high}]", digits = 3, model.names = c("model 1c", "model 2c"))# Compare asymptotic distributions of coefficients:
p1 = plot_summs(gfit1c, gfit2c, scale = TRUE, plot.distributions = T, model.names = c("model 1c", "model 2c"))
p2 = plot_summs(gfit1c, gfit2c, scale = TRUE, plot.distributions = F, inner_ci_level = .95, model.names = c("model 1c", "model 2c"))
grid.arrange(p1, p2, ncol=2)
比较模型系数的渐近分布(缩减数据集模型与估算数据集模型)
正如已经提到的，基于可能性的分数在不同数据集上拟合的模型之间是不可直接比较的，尽管如此，根据我们的经验法则，可以有把握地说 0.56 的麦克法登分数本身就表明非常好的拟合。将估算数据集模型系数的渐近分布绘制在简化模型系数的渐近分布旁边是一种比较标准误差并确认两种方法之间没有根本变化的好方法。所有预测因子的估计系数在两种方法中都被发现是显著的，所以我们很高兴地观察到他们的推断保持一致。性别系数的点估计值在拟合简化数据的模型中较高，但其置信区间与拟合估算模型的置信区间重叠，因此两个估计值之间的推断没有显著差异。
到目前为止，我们有足够的证据相信我们的缺失数据实际上是“随机缺失”(MAR)，即一个值缺失的概率只取决于观察值，而不取决于未观察值。如果 MAR 假设成立，那么 MICE 在产生正确表示“真实”缺失数据的模拟观察方面做得很好。因此，MICE 估算数据集上的模型与缩减(列表式删除)数据集上的模型一致，因此，可以安全地假设两个数据集都是人口的代表性随机子样本。理论上，我们可以使用两个数据集中的任何一个进行下一步的分析，但是出于信息增益的原因，我们将使用估算数据集。
解决了估算数据集后，我们可以进入模型选择的最后阶段，这需要创建和添加新的预测值。我们注意到，变量‘T2’的加入是从。武汉对模型的拟合通过减少 BIC 而显著提高，伪 R 增加约. 20。在这一点上，我们不清楚为什么一个表明一个人是否出生在武汉的变量会与死亡率有任何关系，更不用说有如此大的影响了。让我们通过测试武汉本地人与非本地人的人均死亡人数来进一步了解这一因素。
武汉本地人死亡率
武汉人和非武汉人的死亡率似乎有显著差异。如果我们相信我们的数据，这种反直觉的洞察力就不应该被忽视，但是这种差异背后的原因是什么呢？一个不幸但现实的假设是，来武汉出差或旅游的人平均收入高于本地人。我们可以通过比较武汉本地人和非本地人的死亡率来挑战这一假设。这可以通过简单的双样本𝑧-比例测试来实现，并可视化为柱状图。
#-------------------------------------------------------------------
#Inspecting proportion of mortality rate in Wuhan natives:
#-------------------------------------------------------------------
# Create dataframe:
tmp = data.frame(rbind(
table(dat1_imp$y[dat1_imp$from.Wuhan==1]),
table(dat1_imp$y[dat1_imp$from.Wuhan==0])
))
names(tmp) = c("total", "deaths")
tmp$death_rate = round(tmp$deaths/tmp$total, 3)
tmp$from_wuhan = as.factor(c(1,0))# Compare proportions (deaths per cases between groups):
se1 = sqrt(tmp$death_rate[1]*(1-tmp$death_rate[1])/tmp$total[1])
# Standard errors of proportions:
se2 = sqrt(tmp$death_rate[2]*(1-tmp$death_rate[2])/tmp$total[2])
tmp$prop_se = round(c(se1, se2), 4)
print(tmp)
print(prop.test(x = tmp$deaths, n = tmp$total, alternative = "greater"))# Barplot of proportions:
ggplot(tmp, aes( y=death_rate, x=from_wuhan)) +
geom_bar(position="dodge", stat="identity", width=0.4,
color="black", fill="cyan", alpha=.2) +
geom_errorbar(aes(ymin=death_rate - prop_se, ymax=death_rate +
prop_se), width=.1, position=position_dodge(.9))
显然，武汉本地人和非本地人之间的平均死亡率有显著差异。让我们探索不同协变量的相同测试，例如两组之间住院前的平均天数:
# Compare average number of days gone by before hospitalization
# between both groups:
d1 = dat1_imp$days_before_hosp[dat1_imp$from.Wuhan==1]
d2 = dat1_imp$days_before_hosp[dat1_imp$from.Wuhan==0]
sem1 = t.test(d1)$stderr
sem2 = t.test(d1)$stderr
tmp$avg_days = c(mean(d1), mean(d2))
tmp$mean_se = c(sem1, sem2)
print(tmp)
t.test(d1, d2, alternative = "greater")# Barplot:
b1 = ggplot(tmp, aes( y=avg_days, x=from_wuhan, fill=from_wuhan)) +
geom_bar(position="dodge", stat="identity", width = .4, alpha=1) +
geom_errorbar(aes(ymin=avg_days - mean_se, ymax=avg_days + mean_se), width=.1, position=position_dodge(.9)) # Boxplot:
df = data.frame(days = c(d1, d2), from_wuhan = as.factor(c(numeric(length(d1))+1, numeric(length(d2)) )) )
b2 = ggplot(df, aes( y=days, x=from_wuhan, fill=from_wuhan)) +
geom_boxplot(outlier.colour="black", outlier.shape=16,
outlier.size=2, notch=T) +
stat_summary(fun=mean, geom="point", shape=23, size=4)
grid.arrange(b1, b2, ncol=2)
武汉市本地人与非本地人平均住院天数的比较
正如所怀疑的，在这个测试中也有一个显著的不同。他们的平均值之间的距离可能很小(3.70 对 2.68)，但他们的差异在统计上是显著的，这进一步加强了我们关于因子’'与’之间潜在联系的假设。武汉’和个人的社会经济地位。至此我们可以推测，变量‘T4’来自何处。武汉作为另一个潜在变量的混杂因素。换句话说，可能有另一个携带所有社会经济信息的预测器。武汉百变。一个合乎逻辑的假设是，这个潜在的预测因素是“国家”，理由是某人的种族与其社会经济地位以及是否是武汉本地人直接相关。我们已经确定不能将’ country ‘添加到模型中，因为它的方差很低，但是我们可以’ from。假设武汉是其代理人。必须有另一个比’ from 更细粒度的潜在变量。‘武汉’与当时的‘国家’相关联。**
我们将引入一个国家层面的宏观经济变量，如 GDP (PPP)，它表示一个国家的人均国内生产总值(按购买力平价计算)。如果各组(来自武汉，而不是来自武汉)之间的每个国家的平均 GDP 值存在显著差异，那么我们关于患者的财富与死亡率相关联的最初假设将进一步加强。最终，这个变量也会提高我们模型的拟合度。
让我们为 GDP 创建一个数据向量，其中元素的顺序对应于国家名称的顺序(按’ levels(dat$country) '列出):
# Adding GDP (PPP):
dat = dat1_imp
dat$gdp = 0
gdp_per_country = c(2182, 16091, 54799, 55171, 51991, 50904, 5004,
52144, 20984, 29207, 14800, 49548, 48640, 55306,
66527, 9027, 17832, 40337, 41582, 46827, 67891,
15599, 34567, 3550, 10094, 30820, 105689, 46452,
43007, 14509, 55989, 67558, 57214, 21361, 70441,
48169, 67426, 8677)
for(i in 1:length(gdp_per_country)){
country = levels(dat$country)[i]
country_inds = dat$country == country
dat$gdp[country_inds] = gdp_per_country[i]
}
名为“ gdp 的新列现在是估算数据的一部分。让我们看看当我们把这个变量引入模型时会发生什么。
模型 3(包含 GDP 的模型)
因此，GDP 模型将具有以下形式:
#----------------------------------------------------------------
# Model 4: Adding GDP (PPP)
#----------------------------------------------------------------# Fit model with GDP:
gfit1d = glm(y~gender+age+days_before_hosp+from.Wuhan+gdp, data = dat, binomial)# Compare models with GDP with and without GDP:# Merge model summaries:
export_summs(gfit1c, gfit1d, scale = F, error_format = "[{conf.low}, {conf.high}]", model.names = c("without GDP", "with GDP"))# Compare asymtotic distributions:
f1 = plot_summs(gfit1c, gfit1d, scale = TRUE, plot.distributions = TRUE, model.names = c("without GDP", "with GDP"))
f2 = plot_summs(gfit1c, gfit2c, scale = TRUE, plot.distributions = F, inner_ci_level = .95, model.names = c("without GDP", "with GDP"))
grid.arrange(f1, f2, ncol=2)# Final model summary:
summ(gfit1d, scale = F, plot.distributions = TRUE, inner_ci_level = .9, digits = 3)
gfit = gfit1d #rename final model
比较模型系数的渐近分布(有国内生产总值与无国内生产总值)
正如我们推测的那样，将 GDP 加入到我们的模型中进一步提高了拟合度，如 BIC 的显著降低和伪 T2 R 的增加所示。它的系数可能很小(可能是由于‘T4’from 的存在)。武汉【我们模型中的 )但它的影响仍然是显著的。我们可以进一步比较每个群体的平均 GDP(武汉本地人与非本地人)，以评估他们的差异是否显著:
# Mean GDP per group (from Wuhan):
d3 = dat$gdp[dat$from.Wuhan==1]
d4 = dat$gdp[dat$from.Wuhan==0]
t.test(d3, d4)
sem3 = t.test(d3)$stderr
sem4 = t.test(d4)$stderr
tmp$avg_gdp = c(mean(d3), mean(d4))
tmp$mean_se_gdp = c(sem3, sem4)# Barplot:
f3 = ggplot(tmp, aes( y=avg_gdp, x=from_wuhan, fill=from_wuhan)) + geom_bar(position="dodge", stat="identity", width = .5) + geom_errorbar(aes(ymin=avg_gdp-mean_se_gdp, ymax=avg_gdp+mean_se_gdp), width=.1, position=position_dodge(.9)) # Boxplot:
df = data.frame(days = c(d1, d2), gdp = c(d3, d4), from_wuhan = as.factor(c(numeric(length(d1))+1, numeric(length(d2)) )) )
f4 = ggplot(df, aes( y=gdp, x=from_wuhan, fill=from_wuhan)) +
geom_boxplot(outlier.colour="black", outlier.shape=16,
outlier.size=2, notch=FALSE) +
stat_summary(fun=mean, geom="point", shape=23, size=4)
grid.arrange(f3, f4, ncol=2)
武汉本地人与非本地人人均 GDP /国家比较
我们的假设检验证实，来自武汉的死亡患者的平均 GDP (PPP)与来自非武汉的死亡患者的平均 GDP (PPP)显著不同。这进一步证实了我们模型的系数 p 值。
到目前为止，我们已经确认估算数据集是“最佳”数据集，并且具有预测值’年龄’、性别’、天 _ 之前 _ 住院日’、的模型来自。武汉’，GDP’，是最接近“真实”模型的一个。作为建模过程的最后一步，我们将使用计算推断来评估系数稳定性和样本的同质性。
自助系数稳定性
到目前为止，我们已经经历了一个保守的模型选择过程，这导致了一个我们认为接近真实模型的模型。验证过程的最后一步是通过自举程序测试模型系数的稳定性。如果我们的系数被发现是“稳定的”,那么我们的样本可以被认为是同质的，我们的估计是可靠的。这里使用两个标准来评估系数的稳定性:
- 计算/袋装系数估计的统计显著性(即当经验分布不包含零时评估)
- 计算/袋装估计值与渐近估计值的接近程度(即当渐近和计算平均值之间的差异分布包含零时进行评估)
请注意，在这种情况下，引导聚合(打包)和重采样的目的与随机森林的目的并不完全相同。我们将使用带替换的随机子采样来训练多个模型，但我们的目标是为模型系数而不是模型预测创建袋装估计和自举置信区间。Bagging 将用于产生模型参数(袋装系数)的点估计，而经验自助分布和那些袋装估计的相关自助置信区间将允许我们评估估计的模型稳定性。在其最简单的形式中，这种评估包括检验袋装估计等于渐近估计的假设，这归结为比较同一参数的渐近分布和 bootstrap 分布。实现这一点的一种方法是测试袋装估计和渐近估计之差等于零的零值。如果我们可以证明渐近估计可能与计算(bootstrap)估计没有不同，我们可以有把握地假设我们已经收敛到“真实”系数。
我们使用自定义函数’ bagged_glm ‘来获得模型参数的 bootstrap 分布。’ alpha 参数用于模型汇总统计，’ class_balancing 标志以 50%对欠采样类进行重新采样(未使用)，’ bag_pct ‘参数确定训练包的大小，标志’ ocv '激活出坏交叉验证(相当于随机森林的出包错误)。
*#-------------------------------------------------*
*# Coefficient stability (final model):*
*#-------------------------------------------------*
bagged_results = bagged_glm(dat = dat, gfit = gfit, alpha = 0.05, class_balancing = F, bag_pct = 0.80, ocv = T)
coef_bags_df = bagged_results$coef_bags_df
让我们来看看计算推理的结果:
*# Bootstrap means of transformed coefficients:*
cat(
"\n Bootstrap means of back-transformed coefficients:\n",
"Intercept: ", mean(coef_bags_df$`(Intercept)`),
"Gender: ", mean(coef_bags_df$gendermale),
"Age: ", mean(coef_bags_df$age),
"Days before hospitalization: ",
mean(coef_bags_df$days_before_hosp),
"From Wuhan: ", mean(coef_bags_df$from.Wuhan1),
"\n"
)
*# Compare with asymptotic values: (merge with table above)*
coefficients(gfit)
confint(gfit)
*# Test null distribution directly: H0: m2-m1 = 0*
par(mfcol=c(3,2))
for(i in 1:length(bagged_results$coef_bags_df)){
d = bagged_results$coef_bags_df[[i]] - coefficients(gfit)[i]
qtt = quantile(d, c(.975, .025))
att = t.test(d)
hist(d, col="lightblue", main = names(coefficients(gfit)[i]),
breaks = 30)
abline(v = 0, col="blue", lwd=3, lty=2)
abline(v = qtt[1], col="red", lwd=3, lty=2)
abline(v = qtt[2], col="red", lwd=3, lty=2)
}
dev.off()
检验模型系数的渐近分布和计算分布之间的差异
正如我们在上面的直方图中看到的，袋装估计值 μ 2 - μ 1 之间的差值 d 的分布在自举置信区间(红线)内包含零(蓝线)。因此，零假设:
在 95%的显著性水平上未被拒绝。尽管未能拒绝零假设并不意味着接受零假设，但这一发现与我们之前的所有发现相结合，应该足以假设我们的估计是可靠的，并允许我们继续进行模型解释。
模型解释
作为分析中期待已久的步骤，我们终于可以继续解释我们模型的系数了。回想一下，任何二项式 logit 模型(等式[1])都可以表示为:
求解 β0 产量:
求解 βj 产量:
最后:
表达式[2]中的 β0 的值可以被解释为当所有连续回归变量为 0 且分类回归变量处于其基础水平时，响应的概率为 1。在我们的场景中，我们可以说“住院前天数、患者年龄和 GDP 为零时的平均死亡概率，并且患者不是武汉本地人”。显然，患者年龄和 GDP 等于 0 是一种无意义的解释，但是，如前所述，我们对该分析中的 β0 的直接解释不感兴趣(否则我们可以使用截距校正技术)。我们对表达式[3]中的 βj 的解释很感兴趣，如果 x 连续，它可以被一般地解释为回归变量 j 每增加一个单位，或者如果 x 是绝对的，它可以被解释为从基础水平到当前水平的转变；其他条件不变。为了使这种解释有意义，我们必须使它适应我们的场景，并简单地把 βj 看作:
死亡率概率的百分比变化，
对于与回归量 j 相对应的任何变化，其他都是相同的。我们可以将这种解释应用于最终模型的系数:
因此，让我们总结一下原始形式和转换形式的模型系数:
glm.transform(gfit, alpha = .10, logodds_str = "exp(beta)%", ci = T, stdout = F)
glm.plot(gfit, alpha = 0.10, logodds_str = "exp(beta)%", stdout = F)
对死亡风险的可变影响的大小
我们可以将转换后的系数作为一个因子，并将 Y 的几率解释为比 X 的几率高几倍，或者我们可以使用百分比变化形式，并说 Y 的概率按照所解释的 *X、*的每一次变化而变化一定的比例。就个人而言，我认为后一种解释更直观，下面我们可以看到它在 90%置信水平下对我们模型的每个系数的应用:
- []—exp(β1):
年龄 每增加一岁，新冠肺炎患者的死亡风险大约增加 11%。(CI: 8.36%-14.97%)
表达这一观点的一种更直观的方式是将转换后的系数乘以因子 10，并得出年龄增加 10 岁的几率增加 114%。因此，我们可以说:
年龄每增加 10 岁，新冠肺炎的死亡风险几乎翻倍。 - []—exp(β2):
新冠肺炎死亡风险男性比女性高 120%。(置信区间:8.67%-367%)** - [days before hosp]—exp(β3):
新冠肺炎患者在他/她出现第一个症状 的第一天之后，患者不住院的每一天，死亡风险都会增加 12%。(置信区间:2.7%-20.9%) - [ 来自武汉]—exp(β4):
武汉本地人的新冠肺炎死亡风险大约高 6 倍 。(置信区间:231%-1577%) - [GDP]—exp(β5):
新冠肺炎的死亡风险间接受到患者原籍国 GDP (PPP)的影响。虽然这种影响太弱，无法量化，但结合武汉的’'，它 很可能暗示了患者的生存机会与其个人财富 之间的依赖关系。
利用模型训练和预测之间的差异。
减少内存占用，提高部署模型的速度和可移植性。
这篇文章的代码可以在这里找到。
最近几个月，我们已经帮助许多公司在各种环境中部署他们的 AI / ML 模型。我们为医疗保健领域的模型部署做出了贡献，在过去的几个月里，我们已经帮助几家公司将经过训练的模型迁移到不同类型的物联网设备上。特别是在后一种情况下，要求通常很严格:计算周期的数量和可用内存通常都是有限的。
在这篇文章中，我澄清了我们如何确保使用标准 ML 库训练的模型，如 PyTorch 、 Scikit-learn 和 Tensorflow 可以在各种边缘设备上有效地部署*。为了使事情变得具体，我们将检查一个简单的逻辑回归模型的训练和部署。然而，我们在这里讨论的大部分直接转移到更复杂的模型。*
模特培训
为了说明模型训练和部署之间的区别，让我们从模拟一些数据开始。下面的代码根据下面的简单模型生成 1000 个观察值:
(一个简单的逻辑回归模型用作本例的 DGP)
*import numpy as np
np.random.seed(66) # Set seed for replication# Simulate Data Generating Process
n = 1000 # 1000 observations
x1 = np.random.uniform(-2,2,n) # x_1 & x_2 between -2 and 2
x2 = np.random.uniform(-2,2,n)
p = 1 / (1 + np.exp( -1*(.75 + 1.5*x1 - .5*x2) )) # Implement DGPy = np.random.binomial(1, p, n) # Draw outcomes# Create dataset and print first few lines:
data = np.column_stack((x1,x2,y))
print(data[:10])*
生成数据后，我们可以专注于拟合模型。我们简单地使用sklearn的LogisticRegression()函数来实现:
*from sklearn.linear_model import LogisticRegression
mod = LogisticRegression().fit(data[:,[0,1]], np.ravel(data[:,[2]]))*
近距离观察
在这一点上，有必要暂停一下，简要地考虑一下引擎盖下正在发生的事情。逻辑回归模型，就像许多其他有趣的 ML 模型一样，被反复训练。为了训练这个模型，sklearn(或者任何其他提供类似功能的包)必须实现几个功能:**
- 某种分数函数表示模型的适合度。这可能是一个误差函数，或最大似然函数。
- 从一次迭代到下一次迭代更新拟合模型参数的函数。
训练过程将有效地重复使用这两个函数:最初，模型的参数是随机实例化的。接下来，检查模型的分数。如果分数被认为是不够的(通常是因为它与前一次迭代相比有所改进)，则更新模型参数并重复该过程。
即使对于这个简单的模型，sklearn也需要多次遍历数据集。下面的代码给出了迭代次数(对于这个种子选择是 7 次):
**# Print the number of iterations
print(f'The number of iterations is: {mod.n_iter_}.')**
因此，为了训练一个模型，我们需要访问数据，几个效用函数，并且我们需要多次迭代/遍历数据集。一般来说，这种训练过程在计算上要求很高，这解释了为什么对于复杂的模型，我们求助于并行计算和 GPU 或 NPU 加速来在合理的时间内完成它。然而，幸运的是，当训练模型时，我们使用的各种 ML 库抽象掉了这样做所需的相当复杂的逻辑。
生成预测
将这与从已经拟合的模型中生成预测进行比较(通常称为推论，但我发现后一个术语令人困惑，因为它在统计学中的用法不同，所以我坚持使用预测)。在拟合模型时，在这种情况下，实际上我们需要生成预测的只是逻辑回归函数(与我们在上面的示例中用于生成数据的数学函数相同)和拟合模型的三个参数。这些很容易检索:
**b = np.concatenate((mod.intercept_, mod.coef_.flatten()))
print(b)**
并且参数最终与我们用于数据生成的值相对接近:[0.84576563 1.39541631 -0.47393112]。
此外，在大多数部署情况下，我们通常只使用一个输入来评估模型:在本例中，是一个长度为 2 的数字向量。所以真的，如果我们想要部署一个模型，我们不需要拟合函数，我们不需要数据，我们不需要迭代。为了生成预测，我们只需要简单有效地实现相关的数学函数。
利用训练和预测之间的差异进行(边缘)部署
“那又怎么样？”你可能会问。当现代模型训练工具抽象掉所有这些细节时，为什么还要关心训练和预测中涉及的本质细节呢？嗯，因为当你希望你的模型被有效地部署的时候，例如当你需要它们在小型设备上快速运行的时候，你最好利用上面描述的差异。
为了便于讨论，请对比以下两种模型部署方法(即，将训练好的模型投入生产，以便您可以使用其预测):
- ****sklearn 作为 REST 服务的 docker 容器部署:这种方法简单且常用:我们启动一个 Docker 映像，其中包含 python 堆栈和用于训练的工具:在上面的示例逻辑回归模型中，
sklearn。接下来，我们创建一个 REST 端点，它使用拟合模型的mod.predict()函数来生成结果。 - ****可伸缩 WebAssembly 部署:最近，但一点也不困难的是将 fitted 模型转换为 WebAssembly(使用类似于可伸缩提供的服务)，并部署。WASM 二进制只包含在最小 WebAssembly 运行时预测所需的逻辑。(自动生成的)二进制文件将只包含必要的逻辑函数(在这种情况下)和估计的参数。二进制文件可能部署在服务器上(因此类似地通过 REST 调用使用)，但是，使用各种可用的运行时，它也可以运行在几乎任何边缘设备上。
显然，第一个部署过程接近数据科学家的“我们所知道的”。直接使用我们熟悉的工具工作是很好的。在许多方面，它是可行的:我们可以通过调用 REST 端点来生成预测。第二种解决方案与我们的标准实践相去甚远，并且对于模型训练毫无用处(也就是说，没有“用于训练模型的 WebAssembly 包……”，如果那句话有任何意义的话)。但是，我们仍然认为应该优先选择它:第二种设置利用了训练和预测之间的差异，在几个方面使模型部署更好:
- 内存占用:**上面两个选项中的第一个将需要至少 75Mb 的容器(让容器变得那么小需要大量的工程设计，更常见的是容器的大小接近 1Gb)。在这种情况下，存储的模型本身很小(~2Kb)，因此容器包含了部署的最大内存块(注意，对于大型神经网络来说，这可能不是真的)。相反，WebAssembly 运行时可以减少到小于 64 Kb 。无可否认，WebAssembly 二进制文件本身比存储的
sklearn模型大(~50kb)，但是它现在包含了生成预测所需的全部内容。因此，第一个部署选项至少需要 75Mb，而第二个部署选项需要不到. 1Mb。** - ****速度:与高效的 WebAssembly 部署相比，消耗运行在 Docker 容器中的 REST 端点(它启动了训练所需的所有东西)在执行时间方面并不占优势。这里是各种模型的一些速度比较，但是，不用说，利用训练和预测之间的差异，仅仅将预测的基本需求投入生产，会将这些预测生成的速度提高一个数量级。
因此，内存占用更少，执行速度更快。这很好，有几个原因:一个原因是，我们可能希望为巴黎协定做出贡献，从而有效地部署模型，而不会在每次生成预测时浪费能量。但是，小的占用空间和快速的执行也很有吸引力，因为这正是我们在将模型投入生产时所需要的:祝您在 ESP32 MCU 板上部署 Docker 容器好运。有了 WebAssembly，这简直是小菜一碟。
放弃
值得注意的是我自己的参与:我是 Jheronimus 数据科学院 的数据科学教授，也是Scailable的联合创始人之一。因此，毫无疑问，我对 Scailable 有既得利益；我有兴趣让它成长，这样我们就可以最终将人工智能投入生产并兑现它的承诺。这里表达的观点是我自己的。
开放式电子健康档案的探索
路易斯·梅伦德斯在 Unsplash 上的照片
放弃
这篇文章的内容是面向学术读者的，我不作任何保证，也不承担任何责任。作为读者，您在探索或以其他方式使用此处包含的材料时，风险完全由您自己承担。您有责任遵守您所在地区、州、省、国家/地区和/或司法管辖区的所有适用隐私法，尤其是关于本文中提及或提及的任何及所有材料。你被警告了。
介绍
在这篇文章中，我将讨论一个公开的电子健康记录(EHR)数据集。作为一个乐观主义者，我希望有一天去身份化的 EHR 能被广泛使用。掌握这些信息的研究人员应该能够更好地测试他们的医疗保健应用。潜在的应用包括训练非线性编程模型，以帮助预测给定某些治疗和程序的患者结果。其他包括改善医疗机构内部和之间的互操作性和信息交换。第三个应用是各种后端的基准测试。最后但同样重要的是，大数据的使用应该有助于为运营和管理改进提供信息，包括更好的护理和服务交付模式。结果是:越来越多的人得到了最佳的医疗保健！
直截了当地说，在写这篇文章的时候，公开的电子病历还很少。我是说，真的很少。许多网站初看起来很有帮助，提供了大量的报告和各种其他研究工具。不幸的是，用于生成这些产品的大部分基础数据是高度模糊的。堪比那么多冰山的尖端！除非您是合作机构或政府机构的成员，否则获取源数据几乎是不可能的。
尽管最初遇到了一些挫折，但我并没有放弃寻找电子病历公共数据集的旅程。我的收集标准包括以下内容:
- 隐私。必须对数据进行适当的去识别，以保护患者及其家属的隐私。除了什么都不做这一显而易见的选择之外，模拟数据是帮助缓解隐私担忧的一个不错的选择。参见此链接了解模拟 EHR。不幸的是，根据研究人员的预期应用，模拟数据可能无法提供足够的真实世界保真度、相关性、噪声或其他标准的分数。第三种选择包括去除身份的患者记录，医疗机构或团体同意在他们通过隐私“嗅探测试”后发布记录。在这些情况下，对于研究人员来说，有这样的决定伴随他们的工作是谨慎的。这可以包括经批准的讨论记录、参考号或其他形式的可审计证据，从隐私角度证明数据的去标识性和可发布性。见以上免责声明。你最好小心行事，因为在你的国家，侵犯病人健康记录的隐私会受到严厉的惩罚。
- 容易接近，如没有任何附加条件的接近。任何提供“合作伙伴计划”和“免费试用”的网站都被排除在我的搜索之外。我的排除列表中的其他网站包括要求人们提供信息请求的网站，由底层组织决定请求是否有效。这些网站通常会在提供数据之前宣传几周到几个月的延迟，如果有的话。
- 这项研究工作的第三个标准涉及互操作性和通过一种或多种已知且受支持的格式集成到数据科学应用程序的容易程度。用更通俗的话来说:如果我能快速、轻松地将其融入 Python 生态系统，那就万事大吉了。这部分是我喜欢称之为“mungieness”的部分。
- 最后，我对包含代表性样本的数据集感兴趣。如果可能的话，它应该基于跨越数年的大量人口。它还应该具有与企业级应用程序生成的文件相当的大小，例如，中型到大型医疗记录系统(MRS)的大小。
我重申一下:这篇文章的目的是什么？以获得包含去识别的 EHR 的代表性数据集。有了这样一个数据集，我打算通过制定和测试几个以医疗保健研究主题为基础的假设来解决问题。我觉得有趣的一些主题包括:
- 一个机构能在多大程度上为患者提供及时有效的服务？
- 管理人员和其他员工是否可以使用任何指标来梳理效率，从而帮助优化护理的提供？
剧透:我不幸发现的数据集并没有帮助我回答这些问题。根据一篇关于 ORBDA:一个用于电子健康记录服务器性能评估的 openEHR 基准数据集的论文的作者，他们研究中的数据旨在用于 OpenEHR 系统的基准测试:“我们不鼓励在基准评估以外的其他情况下使用 ORBDA。”[1]尽管如此，庞大的信息量(以 GB 为单位)意味着我至少可以使用我认为非常现实的医疗保健信息技术副产品来测试我的一些分析工具集。请查看此链接以访问数据集。
程序
这些数据可以在几个文件中下载，并有附带的 SQL 脚本来帮助研究人员将其加载到 PostgreSQL 数据库中。我的测试是在运行 Ubuntu 20.04 的笔记本电脑上进行的。
首先，您需要安装 PostgreSQL:
sudo apt-get install postgresql
我建议打开两个终端窗口，一个运行 shell 命令，另一个在登录到数据库后从 postgresql 提示符运行命令。
从这个网站下载 ORBDA 数据文件和脚本，并将它们放在一个合适的文件夹中。
接下来，解压缩脚本和数据文件。
为了让 PostgreSQL 在您的本地机器上启动并运行，可以在这里找到一个有用的教程。使用默认用户名登录。这一步需要您的 Ubuntu 帐户的管理员密码:
sudo su postgres
在第二个终端窗口中执行上述步骤。登录后，运行:
psql
提示符应该变成类似如下的内容:
postgres=#
使用这个新提示，您现在可以输入 PostgreSQL 命令来与您的数据库进行交互。
在创建数据库之前，打开名为“orbdaCreateDB.sql”的脚本，并根据您所在的地区更改排序规则，否则会出现错误。对于我所在的地区，我必须将其从:
LC_COLLATE = ‘en_US.UTF-8’
LC_CTYPE = ‘en_US.UTF-8’
收件人:
LC_COLLATE = ‘en_CA.UTF-8’
LC_CTYPE = ‘en_CA.UTF-8’
更改后，保存脚本。切换回第一个终端窗口，运行以下脚本。
psql -U postgres -f orbdaCreateDB.sql
如果它成功完成，您现在应该有一个名为“orbda”的数据库。在 PostgreSQL 提示符的第二个终端窗口中，键入:
\c orbda
您应该会收到一条消息，说明:
You are now connected to database “orbda” as user “postgres”.
orbda=#
您在该提示符下运行的后续命令应该在 orbda 数据库上执行。
切换回第一个终端并运行表创建脚本:
psql -U postgres -d orbda -f orbdaCreateTables.sql
如果一切顺利，您可以切换到第二个终端来查看新创建的表的列表。在 postgresql 提示符下键入’ \dt '以查看表列表:
orbda=# \dt
List of relations
Schema | Name | Type | Owner
— — — — + — — — — — — — — -+ — — — -+ — — — — —
public | bariatrics | table | postgres
public | chemotherapy | table | postgres
public | hospitalisation | table | postgres
public | medication | table | postgres
public | miscellaneous | table | postgres
public | nephrology | table | postgres
public | radiotherapy | table | postgres
(7 rows)
如果你已经达到了这篇文章的这一点，你就可以开始玩一些有趣的游戏了，叫做等待游戏。在此之前，请打开“importingFiles.sql”脚本，将目录更改为您在上述步骤中放置文件的位置。打开“importingFiles.sql ”,并根据需要将路径修改为文件的绝对路径。如果不这样做，在尝试导入数据时会出现错误。
更新并保存脚本后，导航到第一个终端并运行:
psql -U postgres -d orbda -f importingFiles.sql
导入脚本运行后，您应该会看到连续写入数据库的条目数量。将鼠标移到第二个终端上，在导入数据时查看表的详细信息。它们的大小应该增加。在短暂停顿之间，在 PostgreSQL 提示符下执行以下命令:
\dt+
一旦导入开始，您可能需要等待几个小时才能完成。您的里程可能会有所不同。
当我发现其中一个文件夹位置命名错误时，一个令人高兴的小后果发生了，导致我在纠正错误后不得不重新运行“导入文件”脚本。我注意到，不仅我试图追加的表变大了，许多其他表也变大了。显然，这些表的一些记录在上次导入时丢失了。第三次运行该脚本导致数据库中没有新条目。
如果一切顺利，您应该有一个 PostgreSQL 数据库，其中有几十亿字节的去标识 EHR。在 postgresql 提示符下使用以下命令检查其中一个表:
SELECT * FROM nephrology LIMIT 1;
结论
这是一个快速入门教程，可以通过链接获得一组公开的电子健康记录。我将重申我上面的友好提醒，这些 EHR 中包含的数据不应用于研究或其他目的，因为研究人员的预期范围是关于 OpenEHR 系统的基准测试。关于将信息从 PostgreSQL 的关系数据库格式转换成 OpenEHR 的更多细节包含在他们的网站上。如果您希望继续转换，建议您遵循附加步骤。
本文到此结束，祝您在医疗保健信息技术领域的数据科学应用一切顺利！
参考
[1] ORBDA: An open EHR 电子健康记录服务器性能评估基准数据集，Douglas Teodoro，Erik Sundvall，Mario joo Junior，Patrick Ruch，Sergio Miranda Freire，2018。
意大利新冠肺炎数据的探索性数据分析
在 Unsplash 上 Caleb Stokes 拍摄的照片
不幸的是，意大利是受疫情病影响最严重的国家之一，所以我试图在现有的数据中找到一些见解。
我使用的数据来自每天更新的官方 git 库“Protezione Civile”，可在这里获得。
用于执行该分析的代码可以在这里找到。
如果你想直接运行笔记本只需点击这个链接:https://my binder . org/v2/GL/acalax % 2f covid-19-EDA-Italy/master？filepath=eda_italy.ipynb
笔记本会下载数据，所以你可以随时更新图表！
免责声明:我没有试图预测或模拟任何事情，只是看着数据，并相应地建立一些图表。
分析覆盖了整个国家，但也可以很容易地扩展到地区和省份，只需修改代码。
我们开始吧。
数据
关于数据的所有信息都可以在存储库的 README.md 中找到，但是可以帮助澄清一些特性是相关的
有些是其他值的和，几乎所有的值都是日复一日累积的。
在代码中，我使用了意大利语的特性名称，而没有翻译它们
我们来详细看看:
被发现呈阳性的人可能处于以下情况之一:
- ricoverati _ con _ sintomi**(有症状住院患者)**
- terapia_intensiva (重症监护)
他们的总数是住院病人总数
- isolamento_domiciliare (居家隔离)
总和为 totale_positivi(当前阳性病例总数)
结果可能是:
- 地美西 _ 瓜里提(已恢复)
- deceduti (死亡)
所有阳性和结果的总和是 total_casi(阳性病例总数)
测试的计数器
- 棉塞(已进行测试)
昨天增加了一个新的字段“nuovi_positivi ”,每天都有新的案例，但我直接通过移动和减去数据得出所有的差异
让我们从显示所有案例总体情况的两个图表开始。
这些图像与 4 月 1 日更新的数据有关。
总阳性病例的总累积计数
并且结果尊重正面的案例
阳性、死亡和恢复的总累积计数
这就是疫情的样子，阳性病例和相关结果大幅增加，幸运的是死亡人数低于痊愈人数。
现在，让我们来看一下阳性病例的累计计数
阳性条件的累积计数明细
3 月 20 日以后，家庭分娩超过住院分娩，而重症监护增加。
让我们看看测试和阳性是如何联系的
测试和阳性的总累积计数
测试和阳性反应越来越多，希望阳性反应会“稳定下来”
现在让我们推导出几个新的特征，死亡与总病例数之间的比率(死亡/总阳性)以及恢复与总阳性之间的比率(恢复/总阳性)
死亡/总阳性比率和康复/总阳性比率
这是最“奇怪”的趋势之一，因为死亡和康复非常接近，而在其他国家，差距更大。
现在，我们来分析每天的增量，即一天与前一天的差值。这让我们更好地了解发生了什么和当前的趋势。
每天新增阳性病例
但愿下降趋势正在发生
案件总数的细目
增量阳性分解
家庭分娩在 3 月 20 日达到高峰，住院分娩在同一天达到最低点
这是 delta outcomes 的细分
Delta 结果分解
回收的价值很高，趋势似乎不错。
让我们看看测试和住院的趋势
重症监护是很小的一部分，让我们来看更多细节
趋势似乎不错，希望是因为没那么必要了。
最后，对所有相关增量进行统计分析。
三角洲的箱线图分析
这是对所发生事情的最好描述。
添加测试缩小了其他数据，但给出了所执行的努力的概念，在一天内执行了超过 35K 个测试，作为最大数量(26/03)
三角洲的箱线图分析
结论
这些只是可能的报道的简单例子，希望可以帮助更容易理解正在发生的事情的规模。
就我个人而言，我没有一个清晰的见解，即使是因为数据的收集方式可能会影响结果，也因为在这个数据集中，目前有一些缺失的功能可能会有所帮助，如关于阳性、康复和死亡的年龄间隔的统计数据。
但是我也在想，使用可用的东西总比什么都没有好，尤其是当我们只是使用“原样”而不猜测任何东西的时候。
最后一个想法:数字可能是不可思议的无菌，但我们不能忘记它们涉及到生命，医院里遭受痛苦的人，哀悼亲人的家庭，处于严格封锁中等待知道自己是否生病的人。
然后还有其他人，包括我自己，没有被任何数据捕获，但仍然经历着这一切。我们真的是同舟共济，齐心协力才能克服这种局面。
注意安全！
编者按: 走向数据科学 是一份以数据科学和机器学习研究为主的中型刊物。我们不是健康专家或流行病学家，本文的观点不应被解释为专业建议。想了解更多关于疫情冠状病毒的信息，可以点击 这里 。
成分数据的探索性分析:农业土壤的案例研究
组成数据是向量，其中所有的组成部分都是正数，并且受到闭包形式的约束，这意味着它们的总和是一个常数。在本文中，我使用土壤矿物部分作为一个案例研究的探索性分析这样的数据。
图片由 Zbysiu Rodak 提供
介绍
背景
农业行业的数据驱动决策在满足不断增长的食品、饲料和原材料需求方面发挥着关键作用，同时确保自然资源和环境的可持续利用。例如，精准农业涵盖了广泛的监测技术，包括数据科学、地质统计学和预测分析等跨学科领域，为农民、农学家、土壤科学家和所有其他利益相关者提供了大量信息，这些信息是我们从未见过的，可用于准确监测土壤和作物状况并分析处理方案。
图片由若昂·马塞洛·马克斯拍摄
在这方面，在目前市场上可用的许多监测技术中，使用台式或手持设备的激光诱导击穿光谱(LIBS)越来越受欢迎，主要是因为其速度快、样品制备最少且成本相对较低。像在农业和食品工业中广泛使用的近红外光谱一样，LIBS 提供了对土壤化学成分的估计，可用于评估其肥力——即土壤中可用的矿物养分，这是为大田作物推荐肥料的基础。
然而，无论是在实验室还是在野外，测量准确度和精度对被分析的土壤类型非常敏感:细质地、中等质地和粗质地土壤。为了克服这些缺点，可以将土壤筛分成细粒或中等粒级，这样可以消除大颗粒对分析的任何影响。另一个选择是把大颗粒打碎。第三种选择是根据土壤质地调整设备参数并进行单独分析。后者是美国宇航局的 LIBS 仪器所采用的，名为 ChemCam ，它有一个机载相机来通知纹理，并配备了一个直径为 300-550 微米的聚焦激光束，可以对大颗粒和非常细颗粒大小的火星岩石和土壤进行元素分析。
美国宇航局火星科学实验室任务的好奇号火星车，使用位于 LIBS 的 ChemCam 仪器分析火星岩石表面的化学成分。图片来源:美国宇航局/JPL 加州理工学院。
在地球上，仍然悬而未决的问题是这些选择将如何影响分析——例如，植物养分的可用性受到包括土壤颗粒大小在内的许多因素的影响——或者它们将如何影响 LIBS 吸引人的特征，如快速响应和低入侵性。在很大程度上，本文将搁置这些问题的细节。相反，主要主题是提供探索性分析的概述，显示土壤颗粒大小如何影响 LIBS 测量或任何基于激光的光谱测量，以便确定未来的挑战，并帮助确定最有效的解决方案。
土壤粗密度
对于任何从事农业的人来说，无论规模大小，了解土壤的质地都是非常重要的，因为它会极大地影响土壤养分和持水能力。土壤质地通常由沉降法确定——使用移液管或比重计——提供土壤矿物颗粒的相对比例，称为土壤颗粒大小分离:沙子(2.0-0.05 毫米)、淤泥(0.05-0.002 毫米)和粘土(<0.002 mm). Sandy soils are often referred to as coarse-textured soils whilst clayey soils are referred to as fine-textured soils. The term “coarse” is used because sandy soils are mostly made up of larger or coarse particles. Clayey soils are mainly composed of fine flaky-shape particles. The term “medium” is used for silty and loamy soils. The latter represents soil that does not have a dominant particle size fraction.
Image by 保罗·莫康
成分数据
众所周知，土壤颗粒大小部分构成了所谓的成分数据，这对其统计分析有很大的影响。事实上，自 20 世纪 80 年代初，在约翰·艾奇逊的工作—**—成分数据的统计分析之后，成分数据被地质学家和统计学家以及其他人熟知为向量，其中分量是作为整体的一部分的严格正数，这意味着每个分量仅携带相对部分的信息。这些数据通常被称为封闭数据，因为它们的总和是一个常数，例如 1%或 100%。从数学上讲，这意味着一个 D 部分的成分数据占据了一个有限的空间，在几何上用( D -1)单纯形表示，例如，成分只能从 0%变化到 100%，就像土壤矿物部分的情况一样。
因此， D 部分组成数据集𝐗表示为:
其中组件xj∀j∈{ 1，2，…， D }受到以下约束:
例如，对于 D = 2 或 1-单形，数据可以表示为一条线段，对于 D = 3 或 2-单形，数据可以表示为一个三角形(所谓的三元图，用于将土壤结构等级表示为沙子、淤泥和粘土百分比的函数)。
探索性分析
在接下来的几节中，我将采取几个步骤来更好地了解数据集。这些步骤包括:(1)收集有关分析样本性质的信息，(2)计算土壤颗粒大小部分的描述性统计，(3)分析土壤质地对 LIBS 信号的影响，以及(4)可视化 PCA 双标图以发现模式。
其余软件包:
*library(data.table)
library(tidyverse)
library(soiltexture)
library(ggtern)
library(compositions)
library(robCompositions)
library(FactoMineR)
library(factoextra)*
土壤样本
我们的数据集包含 LIBS 分析的 172 个土壤样本，其中包含土壤质地信息，即沙子、淤泥和粘土的百分比。
*> soil.data <- fread(file="data.csv") %>%
unique() %>%
as_tibble() %>%
str()
Classes ‘tbl_df’,‘tbl’ and 'data.frame': 172 obs. of 7156 variables:
$ SampleID : Factor w/172 levels "S18-0001",..: 1 2 3 4 5 ...
$ Clay : num 41 35.1 81 69.1 39.1 ...
$ Silt : num 30 41.9 12 23.9 41.9 ...
$ Sand : num 29 23 7 7 19 19 26.9 ...
$ 199.3771616: num 659 659 664 663 661 ...
$ 199.4644141: num 664 668 660 657 662 ...
$ 199.5516666: num 674 669 660 661 659 ...
$ 199.6389192: num 663 671 658 667 661 ...
[list output truncated]*
我们的第一步是收集有关分析样本的信息，以了解数据集的性质和分布。下图显示了一个三元图(或土壤结构三角形),它告诉我们分配给每个样品的结构类别以及它们在不同类别中的分布情况。
三元图用于报告土壤颗粒大小分离的相对比例。
然后，我们可以根据土壤样品的名称来确定它们的不同结构类别。根据美国农业部的分类系统，共有 12 种土壤质地类别，土壤质地类别符号见下表。
*> TT.classes.tbl(class.sys="USDA-NCSS.TT") abbr name points
[1,] "C" "clay" "24, 1, 5, 6, 2"
[2,] "SIC" "silty clay" "2, 6, 7"
[3,] "SC" "sandy clay" "1, 3, 4, 5"
[4,] "CL" "clay loam" "5, 4, 10, 11, 12, 6"
[5,] "SICL" "silty clay loam" "6, 12, 13, 7"
[6,] "SCL" "sandy clay loam" "3, 8, 9, 10, 4"
[7,] "L" "loam" "10, 9, 16, 17, 11"
[8,] "SIL" "silty loam" "11, 17, 22, 23, 18, 19, 13, 12"
[9,] "SL" "sandy loam" "8, 14, 21, 22, 17, 16, 9"
[10,] "SI" "silt" "18, 23, 26, 19"
[11,] "LS" "loamy sand" "14, 15, 20, 21"
[12,] "S" "sand" "15, 25, 20"*
我们还可以检查分配给土壤样品的不同结构类别的数量。结果是我们有 11 种土壤质地——我们没有砂质粘土。如下图所示，最丰富的是粘壤土、壤土和砂壤土，其次是粉质粘壤土、粉质壤土和粘土。
*soil.texture %>%
ggplot() +
geom_bar(aes(x=class, fill=class, weight=num.sample)) +
geom_text(aes(x=class, y=num.sample, label=num.sample), nudge_x=-0.1, nudge_y=1) +
scale_fill_hue() +
labs(x=" ", y="Total number of soil samples") +
theme_grey(base_size=14) +
theme(legend.position="none")*
分配给每个质地类别的土壤样本总数。
每一类可以根据它们的纹理进行分组，例如粗糙、中等或精细，这对于主成分分析是有用的。
分为粗、中和细质地的土壤质地类别的数量。
此外，我们还可以从下面的三元图中看到，我们的大多数土壤样品具有大致相等比例的沙子、淤泥和粘土大小的颗粒，这解释了不同壤土类型结构在结构类别中的优势。
*soil.data %>%
ggtern(aes(x=Sand, y=Clay, z=Silt)) +
stat_density_tern(aes(fill = ..level.., alpha = ..level..), geom = "polygon") +
scale_fill_gradient2(low="red", high="blue") +
theme_bw(base_size=14) +
guides(color="none", fill="none", alpha="none")*
密度等值线的可视化。
土壤颗粒级分
由于土壤颗粒的成分性质，土壤颗粒数据存在固有的依赖性。应用标准统计方法——例如算术平均值、标准偏差等——无法解释数据的这种受限属性。这就是为什么像探索性数据分析(EDA)一样，开发了一套称为成分数据分析(CoDA)的程序来分析这种数据。在 CoDA 中，对数比变换应用于成分数据，因此数据不受限制，可以在真实空间中取任何值，并且可以使用任何标准统计方法进行分析。不同类型的对数比变换包括:
加法对数比变换(alr)
居中对数比变换(clr)
等距对数比变换(ilr)
对于 D = 3 部分组成数据，加法、居中和等距对数比转换由下式给出
因此，下文通过箱线图矩阵总结了土壤颗粒大小部分的描述性统计数据，其中每个箱线图都是成对的对数比。
成对对数比的箱线图矩阵。
或者，我们可以计算一些描述性的统计数据。
*> comp.data %>% mean() # Compositional (geometric) mean
Clay Silt Sand
0.2413 0.3579 0.4008> comp.data %>% variation() # Variation matrix
Clay Silt Sand
Clay 0.0000 0.4174 1.2623
Silt 0.4174 0.0000 0.9059
Sand 1.2623 0.9059 0.0000> comp.data %>% var() # Variance matrix of the clr-transform
Clay Silt Sand
Clay 0.3074 0.0927 -0.4001
Silt 0.0927 0.2371 -0.3298
Sand -0.4001 -0.3298 0.7298*
土壤质地对 LIBS 光谱的影响
在查看了我们的土壤样本及其结构数据后，我现在将考虑 LIBS 光谱数据。这里，要考虑的两个最重要的参数是光谱的整体强度及其基于重复测量的相对标准偏差——我们对每个分析样品进行了 24 次重复测量。
下图说明了从不同土壤质地类别获得的平均强度，并根据土壤质地类型进行分组。
*spec.class %>%
arrange(norm.int) %>%
mutate(class=factor(class, levels=class)) %>%
ggplot(aes(x=class, y=norm.int, group=texture, color=texture)) +
geom_point() +
geom_line() +
labs(x=" ", y="Normalized averaged intensity") +
theme_bw(base_size=14)*
作为土壤质地函数的 LIBS 光谱的标准化平均强度。
类似地，比较信号的可变性。
作为土壤质地函数的 LIBS 信号的相对标准偏差。
主成分分析
对光谱数据进行 PCA，其清楚地显示了样品按结构类型聚类，尽管按结构类别聚类不太明显。
*set.seed(1234)
pca.model <- soil.data[ ,-1] %>%
select(Clay, Silt, Sand, colnames(spec.diff[selec.wave])) %>%
PCA(scale.unit=TRUE, graph=FALSE)fviz_pca_biplot(pca.model,
axes = c(1,3),
geom.ind = "point",
geom.var = c("arrow", "text"),
pointshape = 21,
pointsize = 4,
label = "var",
select.var = list(name=c("Clay", "Silt", "Sand")),
fill.ind = pca.data$texture,
col.ind = "black",
col = "darkred",
palette = "jco",
legend.title= "Texture",
repel = TRUE) +
labs(x="t1 (69.3%)", y="t3 (5.8%)") +
theme_bw(base_size=14)*
按土壤质地类型分组的光谱数据的 PCA 双标图。
按土壤结构类别分组的光谱数据的 PCA 双标图。
摘要
我提供了一个快速概览，介绍了使用 LIBS 仪器对土壤样品进行元素成分分析的准确性如何受到不同类型土壤质地的影响。我们以前认为，特别注意土壤颗粒大小的不均匀性会产生更好的测量精度。LIBS 受欢迎的原因之一是其速度和多功能性，因此将样本接地是一个最佳的选择，但可能更具理论性。
**更新:R 代码可从以下 GitHub 链接访问，https://GitHub . com/ChristianGoueguel/EDA-of-composition-Variables
对 1959 年至 2016 年美国宇航局宇航员的探索性分析
当我们还是孩子的时候，这个世界似乎充满了无限、无尽和奇妙的可能性，我们被问及长大后想做什么。虽然有许多不同的反应，如医生，消防队员，警察，或美国总统，其他人想成为宇航员，为美国航天局工作。虽然成为一名宇航员不是一件容易的事，但这篇文章将探索那些已经成为宇航员的人，并希望鼓励那些正在努力实现他们太空雄心的人。
我们会看母校的，军衔，军兵种，军种分布，性别分布，最致命的任务，最多的太空飞行，在太空呆的时间最多。
宇航员的性别分布是怎样的？
作者图片
大多数是男性，但女性不要气馁！
哪些学校培养了最多的宇航员？
作者图片
我们可以看到美国海军学院&美国空军学院培养了最多的宇航员。普渡大学和麻省理工学院紧随其后。
宇航员的本科学位是什么？
作者图片
许多宇航员拥有不止一个本科学位。
他们的研究生学位呢？
作者图片
许多宇航员也有不止一个研究生学位。
宇航员中有百分之多少是军人？
作者图片
大多数是军人，但也有许多不是军人。
哪些军种？
作者图片
似乎美国空军和美国海军生产的宇航员最多。
军衔？
作者图片
他们中的许多人过去/现在都是军队中的高级成员。
五次或更多次太空飞行的宇航员。
作者图片
富兰克林·r·张-迪亚兹 & 杰里·l·罗斯并列第一。
太空飞行时间最多的前 20 名宇航员。
作者图片
杰弗里·n·威廉姆斯赢得最多太空飞行小时。
最后，哪个任务是最致命的？
作者图片
1986 年 1 月 28 日，STS-51L 上的挑战者号及其七名宇航员失踪。他们的名字分别是格雷戈里·b·贾维斯，s·克丽斯塔·麦考利夫，罗纳德·e·麦克奈尔，埃里森·s·奥尼祖卡，朱迪思·a·雷斯尼克，弗朗西斯·r·斯科比，迈克尔·j·史密斯。愿所有在任务中丧生的宇航员安息。
总之，
根据我的发现，如果你正在考虑成为一名宇航员，这可能会增加你的机会，加入军队并真正将自己融入高等教育是明智的。虽然工程、数学和科学主要是在宇航员中发现的，但我相信 NASA 可以找到许多对他们的任务有用的不同领域。
这个数据集来源于Kaggle.com的美国宇航局宇航员年鉴。
探索性数据分析
视觉化！视觉化！视觉化！哦，不要浪费墨水！
图片由皮克斯拜的 Gerd Altmann 提供
数据质量
我们以前写过关于数据质量的问题。当然，许多问题可以在一开始就迅速而清晰地被发现，例如那些缺失的河流水位值(即使您不能立即确定如何处理它们)。但是，许多问题只有在你稍微深入挖掘数据，并开始考虑变量内部和变量之间的分布和关系时才会浮出水面。例如，一个人的身高可能非常合理地介于 50 厘米(学步儿童)和 1.80 米之间。他们的体重也可能非常合理地介于 10 千克和 100 千克之间。但是只有把两者放在一起考虑，你才能意识到有一个 100 公斤重的“蹒跚学步的孩子”或者一个 15 公斤重的高个子成年人。你手上突然有了一些离群值！因此，探索性数据分析(EDA)是数据质量问题仍处于探索阶段的一个步骤。
复杂的关系
我们刚刚举了一个例子，在这个例子中，你需要一致地比较两个变量，以便发现异常值。言下之意，这很容易；你只需在 x 轴上画出一个变量，在 y 轴上画出另一个变量，然后看看“突出”的是什么。如果你也有年龄呢？有一些绘图工具可以让你定义一个 z 轴，甚至可以动态地旋转 z 轴来显示不同方向的数据。如果你也有薪水呢？现在变得棘手了，不是吗？你可以在多张图上画出变量的不同组合，比如 x 和 y。但即使这样也只是考虑了成对的关系。
如果你的数据有 10 个变量呢？20?1000?幸运的是，有许多技术可以减少数据的维数并可视化模式。一种常见的技术是主成分分析(PCA)，尽管还有许多其他技术。
形象化
谈到 EDA，一项我们怎么强调都不为过的技能是数据可视化。这并不总是一个数字！如果它是合理的少量汇总统计数据，那么表格通常是最有效的交流方式。现在我们不太倾向于打印这么多，但是想象一下，打印一个中等大小的条形图，仅仅传达三个数据值，要用多少墨水。也许这将是一个很好的候选，在一个简单的表中提出。
现在软件包中有越来越多的绘图类型，提供了许多不同的方式来表示您的数据。但是仅仅因为它们的存在，并不意味着你需要使用它们。不要像 datavis 那样，使用他能找到的每一个动画和声音效果进行 powerpoint 演示，因为它们就在那里！仔细考虑你想在你的形象中引出什么信息，并据此精心制作。让您入门的一些基本类型:
- 散点图—非常适合显示两个变量之间的关系。可以有效地利用点的颜色、大小甚至形状。考虑将哪些功能映射到哪种美学。除非你确信你的读者有非常敏锐的眼光，否则不要把一个有 20 个不同值的特征映射到形状上！你真的想用 20 种不同的值来表示吗？或许只给前五个类别涂上颜色，其余的归入“其他”类别。不同的颜色可以提供离散(分类)特征的信息。当心有太多！色调(单一颜色，例如从黑色到越来越浅的灰色)可以为你的身材增加一个额外的连续变量。许多图形库会根据数据类型选择颜色或色调。注意你的绘图函数认为一个特征是什么数据类型！1/0 代表真/假是一个经典的例子，它可以表现为一种色调，而当它作为不同的颜色会更清晰。最重要的是，探索不同的选择，有效地展现你的信息。如果你有大量密集的数据点，就试试 alpha 值，问问自己是否真的需要将每个数据点都放在页面上——你能通过在每一千个数据点中绘制一个来显示关于数据分布的相同信息吗？
- 条形图——显示一组离散值的常用工具，例如类别计数。如果 x 轴也是一个类别，不要因为它产生了一道美丽的彩虹就分别给每个类别上色！例如，假设 x 轴是月，y 值是蝴蝶数量。不要按月份给酒吧涂颜色。但是也许应该考虑按季节上色。这可能会增加一些东西。或者，根据一个月是特别潮湿还是特别干燥，它可能与颜色有关。这将是对色彩的有效利用。如果你有太多的类别(在 x 轴上)以至于标签难以辨认，问问你自己是否真的需要显示所有的类别。也许仅仅一个显示前 20 个类别的数字就能表达你的观点？在这里，我们来到了“浪费墨水”点。如果你有两个类别，比如男性和女性，那么也许只需要把它们放在一个表格里。
- 直方图——非常适合汇总分布。一定要探索不同的容器宽度，以便在数据中提取出适量的有意义的结构。不，我不能告诉你那是什么。
- 折线图-非常适合绘制连续的数值序列。如果中间的 x 值没有意义，例如在处理类别时，请避免使用它们。眼睛自然会认为墨水连接点的轨迹是有效的插值。如果你没有很多 x 值，重新考虑它们的用途。包含三个顶点的线图信息不多。一个个人的 bugbear 是由许多值(点)组成的线，但它们形成了非常非常接近直线的东西。该线图是否仅包含两点，两端各一点，在这种情况下，这是一个很差的图，给人一种高相关性的错觉，还是有 1000 个点，这是(非常可疑的)相关性的更强有力的证据。散点图在这里会更好吗？或者至少是有点有线的剧情！用不同类别着色的多条线组成的单个图可以有效地突出不同类别之间的不同趋势。
这是一个令人难以置信的情节类型的节略选择，但它们是很好的中坚力量，通过一些仔细的思考，你可以一次又一次地有效利用。最重要的是，这一点怎么强调都不为过，想想你在你的情节中试图传达的是什么，并试图有效地传达它，只有它，用最少的笔墨和无关的大惊小怪和麻烦。
关于这篇文章
这是一个链接系列的第四篇文章，旨在简单介绍如何开始数据科学过程。这里可以找到介绍，这里可以找到上一篇，这里可以找到系列的下一篇。
探索性数据分析:DataPrep.eda 与 Pandas-Profiling
使用正确的工具进行探索性数据分析(EDA)
探索性数据分析(EDA)是每个数据科学项目的重要组成部分。EDA 的目的是实现对数据的理解，并获得关于数据所代表的现象的见解。
Pandas-profiling (2016)被誉为进行 EDA 的典范工具[1，2，3]。然而，pandas-profiling 的一个显著缺点是它给出了数据集的轮廓！EDA 是一个 迭代 过程，数据科学家将在其中质疑、理解、处理、转换数据，并重复[4、5、6]。pandas-profiling 的僵化结构与当前的 EDA 最佳实践相悖。
DataPrep.eda (2020)是由 SFU 的数据科学研究小组制作的用于做 eda 的 Python 库。DataPrep.eda 支持迭代和以任务为中心的分析——这是 eda 应该做的。(关于 DataPrep.eda 的全面介绍，请参见本文
在这篇文章中，我们将分析 DataPrep.eda 比 pandas-profiling 更适合做 eda 的四个原因:
“‘探索性数据分析’是一种态度，一种灵活的状态，一种寻找我们认为不存在的东西以及我们认为存在的东西的意愿。”
—《探索性数据分析》的作者约翰·图基
1.更好的 API 设计
我们将使用 DataPrep.eda 中的plot()函数。为了理解如何使用该函数有效地执行 eda，下面给出了数据科学家意图的函数调用的语法:
plot(df):“我想要数据集的概览”plot(df, “col_1”):“我想了解栏目col_1”plot(df, “col_1”, “col_2”):“我想了解col_1和col_2列的关系。”
为了了解这一点，我们将使用一个包含南韩新冠肺炎患者记录的数据集。让我们从数据集的概述开始:
from dataprep.eda import plot
import pandas as pddf = pd.read_csv("PatientInfo.csv")
df["confirmed_date"] = pd.to_datetime(df["confirmed_date"])
plot(df)
每列的频率分布
请注意列birth_year似乎具有双峰分布，让我们了解一下该列的更多信息:
# transform birth_year to age in 2020 for simplicity
df["age"] = 2020 - df["birth_year"]
plot(df, "age", bins=26)
以多种方式绘制的一列分布图
使用工具提示，我们可以检查模式的边界，并使用各种图来全面了解这种分布。
接下来，让我们调查一下新冠肺炎男性和女性的年龄分布。为此，我们只需将列sex添加到前面的函数调用中:
plot(df, "age", "sex", bins=26)
按性别比较年龄分布的图表
我们看到在感染新冠肺炎病毒的男性和女性中，年龄分布有很大的差异——所有的差异都来自一行简单的代码！
这说明了 EDA 应该如何进行——提问、想象、理解、重复。用不到 0.5 秒的时间完成上述每个命令，DataPrep.eda 是一个高效的 eda 工具。
pandas-profiling 能支持这个简单的 EDA 任务吗？
我们运行以下代码，其中 pandas-profiling 花费了 50 秒来生成报告:
from pandas_profiling import ProfileReport
ProfileReport(df).to_widgets()
为了分析单变量分布，pandas-profiling 具有以下输出:
熊猫 Jupyter 笔记本中的剖析 API
用户需要切换每一列来查看它的信息。虽然我们可以发现birth_year的双峰分布，但是 pandas-profilin g 不支持对这一见解的进一步研究，也不支持分析age和sex之间的关系。
对于这个简单的 EDA 场景，DataPrep.eda 能够发现重要的见解，但是 pandas-profiling 是不够的。
交互式功能—工具提示
用户与数据的交互对于有效理解数据是必不可少的[7]。DataPrep.eda 使用交互式可视化库 Bokeh 创建所有图，并包含工具提示，可准确读取可视化的每个组件。然而，熊猫烧香不支持工具提示。
2.速度快 100 倍
DataPrep.eda 比 pandas 快 100 倍——针对 eda 的分析
回想一下，在第 1 部分的中，完成每项任务需要不到 0.5 秒的时间，而 pandas-profiling 需要 50 多秒才能生成一份报告。
DataPrep.eda 比 pandas 更快，原因有二:
- DataPrep.eda 使用 Dask ，一个并行计算库，用于所有的数据处理。然而，熊猫简介使用了熊猫。
- DataPrep.eda 通过创建与当前 eda 任务相关的可视化来避免不必要的计算，而 pandas-profiling 只分析整个数据集。
DataPrep.eda 比 pandas 快 5 倍以上——即使对于数据分析也是如此
DataPrep.eda 可用于生成 pandas-profiling 报告的组件，从而实现直接的性能比较。下图显示了针对 DataPrep.eda 的三个组件运行 pandas-profiling 的ProfileReport的结果:单变量分析(plot(df))、相关矩阵(plot_correlation(df))和缺失值图(plot_missing(df))。
DataPrep.eda 与 pandas-profiling 的等待时间(使用汽车数据集，通过复制放大)
使用 Dask 代替 Pandas 是 DataPrep.eda 比 pandas-profiling 更快的主要原因。影响性能的更具体的因素是
- DataPrep.eda 并行处理单变量分析，而 pandas-profiling 顺序计算单变量统计。
- 使用 Dask 的 DataPrep.eda 支持逐块计算，而 Pandas-profiling 在整个数据集上执行计算(对于大型数据集很重要)。
3.智能可视化
DataPrep.eda 的一些智能特性包括
- 选择正确的图来可视化每个 EDA 任务的数据；
- 列类型推断(数字、分类和日期时间)；
- 为每个图找到合适的时间单位(用户也可以指定)；
- 输出具有最高计数的分类值，以便视觉清晰(用户也可以指定)。
为了看到这些功能的作用，让我们理解人们是如何随着时间的推移收缩新冠肺炎的，即列confirmed_date和infection_case之间的关系。我们跑
plot(df, "confirmed_date", "infection_case")
描绘日期时间和分类列之间关系的图
我们可以很容易地看到哪些收缩方法是重要的，在哪些时期！
我们可以用熊猫轮廓来完成这项任务吗？
不要！Pandas-profiling 只支持相关矩阵形式的交互(也受 DataPrep.eda 支持)和用于两个连续变量的双变量分析的热图。在数据集分析框架中，有效的双变量分析代价太高，因为必须对每一对列进行计算，即使用户可能只对一小部分关系感兴趣。
4.处理大量数据
使用 Dask 的 DataPrep.eda 适用于大于内存的数据集。Dask 支持核外并行处理，因此可以高效地评估大型数据集上的计算。
使用 Pandas 的 Pandas-profiling 只有用于内存分析的数据结构；在大型数据集上，pandas-profiling 的性能会显著下降。
结论
探索性数据分析是一个迭代循环，其步骤包括
- 质疑数据
- 通过处理和可视化数据来回答问题
- 在获得新的理解后提炼之前的问题，或者提出新的问题
没有适合全面 EDA 的通用数据配置文件。
与 pandas-profiling 相比，DataPrep.eda 是更好的 eda 工具，原因有四:
- 更好的 API 设计
DataPrep.eda 的 API 是为 eda 设计的，而不是为数据剖析设计的 - 快 100 倍
DataPrep.eda 并行执行计算 - 智能可视化
DataPrep.eda 将自动选择合适的图来可视化数据 - 处理大数据
DataPrep.eda 支持核外处理
是时候从生成数据概要文件开始，以 DataPrep.eda 的方式执行 EDA 了。
一个笔记本的代码从这篇文章中可以找到 这里 。要安装 DataPrep.eda，并了解有关参与该项目的信息，请访问 此处 。一个 DataPrep.eda 教程视频可以在 这里 。别忘了给 GitHub 上的 项目 上星★。
参考
[1] M. Deep，快速探索性数据分析:熊猫简介 (2020)，中等
[2] L. Frei，使用 Pandas-Profiling (2019)加速您的探索性数据分析，走向数据科学
[3] R. Rei， EDA 使用 Panda 的剖析 (2020)，走向数据科学
[4] D. Bourke，探索性数据分析的温和介绍 (2019)，走向数据科学
[5] J. Wei，探索性数据分析:结构化数据实用指南和模板 (2019)，走向数据科学
[6] G. Grolemund 和 H. Wickham 著，R for Data Science(2016 年 12 月)，在线图书
[7] W. Koehrsen，Python 中数据可视化的下一个级别 (2019)，走向数据科学
使用 EarthPy 对卫星图像进行探索性数据分析
中间导轨
在本文中，我们将使用 EarthPy 来处理卫星图像，并有效地进行探索性数据分析(EDA)
杰西·艾伦和罗伯特·西蒙利用美国地质调查局地球探测器的数据拍摄的陆地卫星图像。
目录
- 卫星图像介绍
- 安装
- 如何下载卫星图片
- 卫星图像上的探索性数据分析(EDA)
- 最终想法
- 参考文献
让我们开始吧，✨
卫星图像介绍
卫星图像有着广泛的应用，它融入了人类生活的方方面面。特别是遥感经过多年的发展，已经解决了不同领域的很多问题。在遥感中，高光谱遥感器以其高光谱分辨率被广泛用于监测地球表面。
高光谱图像(HSI) 数据通常包含同一空间区域内的数百个光谱带，这些光谱带为识别各种材料提供了有价值的信息。在 HSI 中，每个像素可以被视为一个高维向量，其条目对应于从可见光到红外的光谱反射率。
遥感的一些最佳应用是矿物勘探、国防研究、生物化学组成、森林健康状况、农业等。使用下面的研究论文可以更直观地了解高光谱遥感的应用。
植被指数叶面积指数高光谱数据悬浮泥沙含量土地类这些…
link.springer.com](https://link.springer.com/chapter/10.1007/978-3-662-47456-3_9)
使用下面的文章获得使用 python 进行高光谱影像分析的实践经验。
使用 Python 进行高光谱图像分析的演练。
towardsdatascience.com](/hyperspectral-image-analysis-getting-started-74758c12f2e9)
装置
让我们看看 EarthPy ，它是一个开源 python 包，使用开源工具可以更容易地绘制和处理空间栅格和矢量数据。Earthpy 依赖于 geopandas，后者侧重于矢量数据，而 rasterio 则有助于栅格数据文件的输入和输出。它还需要 matplotlib 进行绘图操作。使用以下命令安装 EarthPy。
如何下载卫星图像
EarthpPy 包中有 14 个可用的数据集，让我们看看可供下载的数据集。
数据集-按作者分类的图像
让我们看看如何下载可用的数据集。“get_data”方法用于使用数据集的名称下载数据。
输出将是:
图片 bt 作者
卫星图像上的 EDA
在本文中，我们使用了’ vignette Landsat’ 数据集。该数据集包含 2017 年 2 月 21 日的 Landsat 8 数据，用于科罗拉多州尼德兰附近冷泉火边界周围的区域。它还包含 GeoMAC 提供的冷泉火边界。Landsat 波段已被裁剪为覆盖冷泉火边界，并且不包括原始图像东南边缘的数据值以及云层。
让我们来看看如何读取数据集:
使用 EarthPy 的空间*模块中的“堆叠”方法选择并堆叠条带。*上面的代码读取数据并打印元数据。输出如下所示。
Landsat8 数据集的元数据-作者提供的图像
数据集的形状为(2158，1941)，有 7 个波段，包含 41，88，678 个像素。
绘图区带
正如我们所讨论的，Landsat8 数据集有 7 个波段。让我们使用 earthpy 包中的内置方法’ plot_bands '来绘制波段。plot_bands 方法采用带和图的堆栈以及自定义标题，这可以通过使用title=参数将每个图像的唯一标题作为标题列表传递来完成。
让我们看看绘制 Landsat8 数据集波段的代码。
结果图如下所示:
landsat8 数据集中所有波段的图-作者提供的图像
RGB 合成图像
这些超光谱图像有多个波段，包含从可见光到红外的数据。所以对人类来说很难将数据可视化。因此，创建 RGB 合成图像可以更容易有效地理解数据。为了绘制 RGB 合成影像，您将绘制红色、绿色和蓝色波段，分别是我们根据 Landsat8 数据创建的影像堆栈中的波段 4、3 和 2。Python 使用从零开始的索引系统，所以你需要从每个索引中减去值 1。因此，红色波段的指数为 3，绿色为 2，蓝色为 1。让我们看看绘制 RGB 合成图像的代码。
生成的 RGB 合成图像如下所示。
作者图片
拉伸合成图像
如果像素亮度值偏向零值，我们创建的合成图像有时会很暗。这种类型的问题可以通过拉伸图像中的像素亮度值来解决，使用自变量stretch=True将这些值扩展到潜在值的整个 0-255 范围，以增加图像的视觉对比度。此外，str_clip参数允许您指定想要剪掉多少数据尾部。数字越大，数据将被拉伸或变亮的越多。
让我们看看如何起诉 earthpy。
应用拉伸后的 RGB 合成图像如下所示:
应用 Strech 后的 Landsat8 数据集 RGB 合成图像-作者提供的图像
绘制光谱
让我们看看如何绘制光谱，这有助于我们理解像素的性质。以下代码用于绘制 Landsast8 数据集第一行的 7 个光谱。
输出如下所示:
光谱图—图片由作者提供
直方图
可视化高光谱图像数据集的波段有助于我们理解波段值的分布。“eathpy.plot”中的hist方法通过为我们之前创建的数据集/堆栈的波段绘制直方图来完成工作。我们还可以修改单个直方图的列大小，标题，颜色。让我们看看绘制直方图的代码。
产生的图是
波段直方图-按作者分类的图像
归一化差异植被指数(NDVI)
为了确定一块土地上的绿色密度，研究人员必须观察植物反射的可见光(VIS)和近红外(NIR)阳光的不同颜色(波长)。归一化差异植被指数(NDVI)通过测量植被强烈反射的近红外和植被吸收的红光之间的差异来量化植被。NDVI 的范围总是从-1 到+1。
NDVI =(近红外-可见光)/(近红外+可见光)
例如，当你有负值时，很可能是水。另一方面，如果 NDVI 值接近+1，则很有可能是浓密的绿叶。但是当 NDVI 接近零度时，就没有绿叶，甚至可能成为城市化地区。
以上代码用于计算归一化差异植被指数(NDVI)并显示生成的图像。
NDVI 陆地卫星 8 号—作者提供的图像
NDVI 的分类
基于高光谱图像数据，归一化差异植被指数(NDVI)结果被分类为有用的类别。0 以下的值将一起归类为无植被。将为裸露区域和低、中等和高植被区域创建附加类别。让我们看看代码:
上述代码对 NDVI 进行分类，并绘制结果数据。分类后的结果图像如下所示。
归一化差异植被指数(NDVI)类别-图片由作者提供
不是但不是最不重要的，我在博客中写的代码可以从下面的 GitHub 链接或使用 google collaboratory 笔记本访问。
[## syamkakarla 98/超光谱 _ 图像 _ 分析 _ 简化
permalink dissolve GitHub 是超过 5000 万开发人员的家园，他们一起工作来托管和审查代码，管理…
github.com](https://github.com/syamkakarla98/Hyperspectral_Image_Analysis_Simplified/blob/master/Articles/Exploratory_Data_Analysis%28EDA%29_on_Satellite_Imagery_Using%C2%A0EarthPy.ipynb) [## 使用 EarthPy 对卫星图像进行探索性数据分析
colab.research.google.com](https://colab.research.google.com/drive/1kpG7Vp_gg5uXUEP0BO3GwN1hMDpPIJKR?usp=sharing)
最后的想法
本文涵盖了使用 EarthPy 分析卫星/高光谱图像的不同方法，但还有更多方法，如降维(DR) 和分类等。使用下面的文章进行详细解释并使用 python 进行实际操作。
高光谱图像的降维技术。
towardsdatascience.com](/dimensionality-reduction-in-hyperspectral-images-using-python-611b40b6accc) [## 超光谱图像分析分类
使用 python 对高光谱图像(HSI)进行分类的演练。
towardsdatascience.com](/hyperspectral-image-analysis-classification-c41f69ac447f)
参考
[## 测量植被(NDVI 和 EVI)
为了监测植被的主要波动并了解它们如何影响环境，科学家们使用…
earthobservatory.nasa.gov](https://earthobservatory.nasa.gov/features/MeasuringVegetation/measuring_vegetation_2.php) [## 地球引擎|地球引擎数据目录中的 Landsat 8 数据集
来自 Landsat 8 OLI/TIRS 传感器数据集的大气校正地表反射率可用性:2013 年 4 月…
developers.google.com](https://developers.google.com/earth-engine/datasets/catalog/landsat-8)
探索性数据分析(EDA): Python
学习使用 Python 和 Numpy、Matplotlib 和 Pandas 进行探索性数据分析的基础知识。
什么是探索性数据分析(EDA)？
如果我们想用简单的术语解释 EDA，这意味着试图更好地理解给定的数据，以便我们可以从中获得一些意义。
我们可以在 维基 中找到更正式的定义。
在统计学中，探索性数据分析是一种分析数据集以总结其主要特征的方法，通常采用可视化方法。可以使用或不使用统计模型，但 EDA 主要是为了查看数据可以告诉我们什么，而不仅仅是正式的建模或假设测试任务。
Python 中的 EDA 使用数据可视化来绘制有意义的模式和见解。它还包括通过消除数据中的不规则性来准备用于分析的数据集。
基于 EDA 的结果，公司也做出商业决策，这些决策会在以后产生影响。
- 如果 EDA 没有正确完成，那么它会妨碍机器学习模型构建过程中的进一步步骤。
- 如果做得好，可能会提高我们接下来做的一切事情的功效。
在本文中，我们将了解以下主题:
- 数据来源
- 数据清理
- 单变量分析
- 双变量分析
- 多变量分析
1.数据来源
数据源是查找数据并将其加载到我们的系统中的过程。概括地说，我们有两种方法可以找到数据。
- 私人数据
- 公共数据
私人数据
顾名思义，私人数据是由私人机构给出的。它有一些安全和隐私问题。这种类型的数据主要用于组织内部分析。
公开数据
每个人都可以获得这种类型的数据。我们可以在政府网站和公共组织等找到这一点。任何人都可以访问这些数据，我们不需要任何特殊的许可或批准。
我们可以从以下网站获取公共数据。
- https://data.gov
- https://data.gov.uk
- https://data.gov.in
- 【https://www.kaggle.com/
- https://archive.ics.uci.edu/ml/index.php
- https://github.com/awesomedata/awesome-public-datasets
EDA 的第一步是数据源，我们已经了解了如何访问数据并加载到系统中。现在，下一步是如何清理数据。
2.数据清理
完成数据源后，EDA 过程的下一步是数据清理。在将数据输入我们的系统后，去除不规则性并清理数据是非常重要的。
不规则性是不同类型的数据。
- 缺少值
- 格式不正确
- 不正确的标题
- 异常/异常值
为了执行数据清理，我们使用一个样本数据集，它可以在 这里 找到。
我们正在使用 Jupyter 笔记本进行分析。
首先，让我们导入必要的库并将数据存储在我们的系统中以供分析。
现在，数据集看起来像这样，
营销分析数据集
如果我们观察上面的数据集，前 2 行的列标题有一些差异。正确的数据来自索引号 1。因此，我们必须修复前两行。
这被称为**固定行和列。**让我们忽略前两行，再次加载数据。
现在，数据集看起来像这样，它更有意义。
修复行和列后的数据集
以下是在固定行和列时要采取的步骤:
- 删除数据集中的汇总行和汇总列。
- 删除每页上的页眉和页脚行。
- 删除多余的行，如空行、页码等。
- 如果有助于更好地理解数据，我们可以合并不同的列
- 同样，我们也可以根据我们的需求或理解，将一个列拆分成多个列。
- 添加列名，将列名添加到数据集中非常重要。
现在，如果我们观察上面的数据集，customerid列对我们的分析不重要，而且jobedu列同时包含了job和education的信息。
因此，我们要做的是，我们将删除customerid列，并将jobedu列拆分为另外两列job和education，之后，我们也将删除jobedu列。
现在，数据集看起来像这样，
删除Customerid 和 jobedu 列，添加 job 和 edu 列
缺失值
如果在进行任何统计分析之前数据集中有缺失值，我们需要处理这些缺失值。
缺失值主要有三种类型。
- MCAR(完全随机缺失):这些值不依赖于任何其他特征。
- MAR(随机缺失):这些值可能依赖于其他一些特征。
- MNAR(非随机缺失):这些缺失值有一些缺失的原因。
让我们看看数据集中哪些列缺少值。
# Checking the missing values
data.isnull().sum()
输出将是，
数据集中的空值
正如我们所看到的，有三列包含缺失值。让我们看看如何处理丢失的值。我们可以通过删除丢失的记录或输入值来处理丢失的值。
删除丢失的值
让我们处理age列中缺失的值。
现在让我们检查数据集中缺失的值。
处理年龄列后缺少值
让我们对 month 列的缺失值进行估算。
因为 month 列是一个对象类型，所以让我们计算该列的模式，并将这些值归入缺失值。
现在输出是，
# Mode of month is
**'may, 2017'**# Null values in month column after imputing with mode
**0**
处理响应列中的缺失值。因为我们的目标列是 Response 列，如果我们将值归入这个列，就会影响我们的分析。因此，最好从响应列中删除缺失的值。
#drop the records with response missing in data.
**data = data[~data.response.isnull()].copy()**# Calculate the missing values in each column of data frame
**data.isnull().sum()**
让我们检查数据集中缺失的值是否已被处理，
所有丢失的值都已被处理
我们还可以将缺失值填充为**‘NaN’**，这样在进行任何统计分析时，都不会影响结果。
处理异常值
我们已经看到了如何修复缺失值，现在让我们看看如何处理数据集中的异常值。
离群值是指远离下一个最近的数据点的值。
有两种异常值:
- **单变量异常值:**单变量异常值是基于一个变量，其值位于期望值范围之外的数据点。
- **多元异常值:**在绘制数据时，一个变量的某些值可能不会超出预期范围，但当您绘制其他变量的数据时，这些值可能会远离预期值。
因此，在理解了这些异常值的原因之后，如果更有意义的话，我们可以通过丢弃那些记录或者用值进行估算或者让它们保持原样来处理它们。
标准化值
要对一组值执行数据分析，我们必须确保同一列中的值应该在相同的范围内。例如，如果数据包含不同公司汽车的最高速度值，那么整列应该以米/秒或英里/秒为单位。
现在，我们已经清楚了如何获取和清理数据，让我们看看如何分析数据。
3.单变量分析
如果我们对数据集中的单个变量/列进行分析，这就是所谓的单变量分析。
分类无序单变量分析:
无序变量是没有定义顺序的分类变量。以我们的数据为例，数据集中的职位列被分成许多子类别，如技术人员、蓝领、服务、管理等。“作业列中的任何值都没有权重或度量。
现在，让我们通过使用图来分析工作类别。由于 Job 是一个类别，我们将绘制条形图。
输出看起来像这样，
通过上面的柱状图，我们可以推断，与其他类别相比，数据集包含更多的蓝领工人。
分类有序单变量分析:
有序变量是那些具有自然排序的变量。我们数据集中分类有序变量的一些示例如下:
- 月份:一月，二月，三月……
- 教育:小学、中学……
现在，让我们从数据集中分析教育变量。既然我们已经看到了条形图，让我们看看饼图是什么样子的。
输出将是，
通过以上分析，我们可以推断出，该数据集有大量属于中等教育之后的那一级和下一级小学。此外，其中很小一部分不为人知。
这就是我们分析单变量分类分析的方法。如果列或变量是数值型的，那么我们将通过计算它的平均值、中值、标准差等来进行分析。我们可以通过使用 describe 函数获得这些值。
data.salary.describe()
输出将是，
4.双变量分析
如果我们从数据集中考虑两个变量/列来分析数据，这就是所谓的双变量分析。
a)数值-数值分析:
分析数据集中的两个数值变量称为数值-数值分析。我们可以从三个不同的角度来分析。
- 散点图
- 配对图
- 相关矩阵
散点图
让我们从我们的数据集中选取三列“余额”、“年龄”和“薪水”，看看我们可以通过绘制salary balance和age balance之间的散点图来推断什么
现在，散点图看起来像，
散点图
配对图
现在，让我们为我们在绘制散点图时使用的三列绘制成对图。我们将使用 seaborn 库来绘制 Pair 图。
结对图看起来像这样，
年龄、平衡、工资的配对图
相关矩阵
因为在散点图和配对图中，我们不能使用两个以上的变量作为 x 轴和 y 轴，所以很难在单个图形中看到三个数值变量之间的关系。在这些情况下，我们将使用相关矩阵。
首先，我们使用年龄、工资和余额创建了一个矩阵。之后，我们使用矩阵的 seaborn 库绘制热图。
b)数值分类分析
分析数据集中的一个数值变量和一个分类变量称为数值分类分析。我们主要使用均值图、中值图和箱线图来分析它们。
让我们从数据集中取出salary和response列。
首先使用groupby检查平均值
#groupby the response to find the mean of the salary with response no & yes separately.**data.groupby('response')['salary'].mean()**
输出将是，
使用平均值的响应和薪金
根据薪水的不同，回答“是”和“否”没有太大区别。
让我们计算一下中位数，
#groupby the response to find the median of the salary with response no & yes separately.**data.groupby('response')['salary'].median()**
输出将是，
无论是平均值还是中间值，我们都可以说，不管一个人的工资多少，回答“是”和“否”的答案都是一样的。但是，它真的是那样表现的吗，让我们画出它们的方框图，并检查它的表现。
#plot the box plot of salary for yes & no responses.**sns.boxplot(data.response, data.salary)
plt.show()**
箱形图看起来像这样，
正如我们所看到的，当我们绘制箱线图时，它描绘了一个与平均值和中值非常不同的画面。给予肯定回答的顾客的 IQR 偏高。
这就是我们分析数字分类变量的方法，我们使用均值、中值和箱线图来得出某种结论。
c)分类—分类分析
由于我们的目标变量/列是回复率，我们将看到不同的类别如教育、婚姻状况等。，与响应列相关联。因此，我们将把它们转换成“1”和“0”，而不是“是”和“否”，这样我们就可以得到“回应率”。
输出看起来像这样，
让我们来看看婚姻状况中不同类别的回复率是如何变化的。
图表看起来像这样，
通过上图，我们可以推断出，对于数据集中的单一状态成员，正面响应更多。类似地，我们可以绘制贷款对响应率、住房贷款对响应率等图表。
5.多变量分析
如果我们从一个数据集中考虑两个以上的变量/列来分析数据，这就是所谓的多元分析。
让我们看看“教育程度”、“婚姻状况”和“回应率”之间的差异。
首先，我们将创建一个包含三列的数据透视表，之后，我们将创建一个热图。
数据透视表和热图如下所示，
根据热图，我们可以推断，受过小学教育的已婚人士不太可能对调查做出积极回应，而受过高等教育的单身人士最有可能对调查做出积极回应。
类似地，我们可以绘制工作对婚姻对反应、教育对收入对反应等图表。
结论
这就是我们如何进行探索性数据分析。探索性数据分析(EDA)帮助我们透过数据看问题。我们对数据探索得越多，从中得出的见解就越多。作为一名数据分析师，我们几乎 80%的时间都会花在理解数据和通过 EDA 解决各种业务问题上。
感谢您阅读和快乐编码！！！
在这里查看我以前关于 Python 的文章
- 使用 Python 对熊猫数据帧进行索引
- Seaborn: Python
- 熊猫:蟒蛇
- Matplotlib:Python
- NumPy: Python
- 数据可视化及其重要性:Python
- Python 中的时间复杂度及其重要性
- Python 中的递归或递归函数
参考
- 探索性数据分析:https://en.wikipedia.org/wiki/Exploratory_data_analysis
- Python 探索性数据分析:https://www . Data camp . com/community/tutorials/explorative-Data-Analysis-Python
- 利用 Python 进行探索性数据分析:https://www . active state . com/blog/explorative-Data-Analysis-using-Python/
- 单变量和多变量异常值:https://www . statistics solutions . com/Univariate-and-multi variable-Outliers/
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
