提交 ef228d90 编写于 作者: W wizardforcel

2022-07-05 10:13:20

上级 7352b98a
此差异已折叠。
此差异已折叠。
此差异已折叠。
此差异已折叠。
# 12 分类关系建模
到目前为止,我们已经讨论了统计建模和假设检验的一般概念,并将它们应用到一些简单的分析中。在本章中,我们将重点介绍 _ 分类 _ 关系的建模,通过建模,我们可以表示在名义(有时是序数)尺度上测量的变量之间的关系。这些数据通常用计数来表示;也就是说,对于变量的每个值(或多个变量的值的组合),有多少观测值取这个值?例如,当我们计算每个专业有多少人在我们班上时,我们将根据数据拟合一个分类模型。
\ No newline at end of file
到目前为止,我们已经讨论了统计建模和假设检验的一般概念,并将它们应用到一些简单的分析中。在本章中,我们将重点介绍 _ 分类 _ 关系的建模,通过建模,我们可以表示在名义(有时是序数)尺度上测量的变量之间的关系。这些数据通常用计数来表示;也就是说,对于变量的每个值(或多个变量的值的组合),有多少观测值取这个值?例如,当我们计算每个专业有多少人在我们班上时,我们将根据数据拟合一个分类模型。
## 12.1 示例:糖果颜色
假设我买了一袋 100 颗糖果,上面写着有 1/3 巧克力、1/3 甘草和 1/3 胶球。当我数包里的糖果时,我们得到以下数字:
```r
candyDf <-
tibble(
candyType = c("chocolate", "licorice", "gumball"),
count = c(30, 33, 37)
)
pander(candyDf)
```
<colgroup><col style="width: 16%"> <col style="width: 9%"></colgroup>
| 糖果型 | 计数 |
| --- | --- |
| 巧克力 | 30 个 |
| 甘草 | 33 |
| 胶球 | 37 岁 |
因为我更喜欢巧克力,而不是甘草或树胶球,我觉得有点被撕了。我想知道的是:如果每种糖果的真实概率是每种糖果平均 1/3 的比例,那么这种计数的可能性有多大?
## 12.2 皮尔逊卡方检验
皮尔逊卡方检验为我们提供了一种方法来检验观察到的计数数据是否与定义零假设的某些特定预期值不同:
![](img/a39dfdf63e525b1d15f99c58b2ae7a2e.jpg)
在我们的糖果例子中,无效假设是每种糖果的比例是相等的。我们可以计算我们观察到的糖果计数的卡方统计,如下所示:
```r
# compute chi-squared statistic
nullExpectation <- c(1 / 3, 1 / 3, 1 / 3) * sum(candyDf$count)
chisqVal <-
sum(
((candyDf$count - nullExpectation)**2) / nullExpectation
)
```
这个分析的卡方统计结果是 0.74,这本身是不可解释的,因为它取决于加在一起的不同值的数量。然而,我们可以利用这样一个事实:卡方统计量是根据零假设下的特定分布分布分布的,这就是所谓的 _ 卡方 _ 分布。这种分布被定义为一组标准正态随机变量的平方和;它有若干自由度,等于被加在一起的变量数。分布的形状取决于自由度的数量。图[12.1](#fig:chisqDist)显示了几种不同自由度的分布示例。
![Examples of the chi-squared distribution for various degrees of freedom.](img/file74.png)
图 12.1 不同自由度的卡方分布示例。
让我们通过模拟来验证卡方分布是否准确地描述了一组标准正态随机变量的平方和。
```r
# simulate 50,000 sums of 8 standard normal random variables and compare
# to theoretical chi-squared distribution
# create a matrix with 50k columns of 8 rows of squared normal random variables
d <- replicate(50000, rnorm(n = 8, mean = 0, sd = 1)**2)
# sum each column of 8 variables
dMean <- apply(d, 2, sum)
# create a data frame of the theoretical chi-square distribution
# with 8 degrees of freedom
csDf <-
data.frame(x = seq(0.01, 30, 0.01)) %>%
mutate(chisq = dchisq(x, 8))
```
[12.2](#fig:chisqSim)显示,理论分布与重复将一组随机正态变量的平方相加的模拟结果非常吻合。
![Simulation of sum of squared random normal variables. The histogram is based on the sum of squares of 50,000 sets of 8 random normal variables; the blue line shows the values of the theoretical chi-squared distribution with 8 degrees of freedom.](img/file75.png)
图 12.2 平方随机正态变量和的模拟。柱状图是基于 5 万组 8 个随机正态变量的平方和;蓝线显示了 8 个自由度下理论卡方分布的值。
对于糖果的例子,我们可以计算在所有糖果的相同频率的零假设下,我们观察到的卡方值为 0.74 的可能性。我们使用自由度等于 k-1(其中 k=类别数)的卡方分布,因为我们在计算平均值以生成预期值时失去了一个自由度。
```r
pval <- pchisq(chisqVal, df = 2, lower.tail = FALSE) #df = degrees of freedom
sprintf("p-value = %0.3f", pval)
```
```r
## [1] "p-value = 0.691"
```
这表明,观察到的糖果数量并不是特别令人惊讶的,基于印刷在糖果袋上的比例,我们不会拒绝等比的无效假设。
## 12.3 应急表及双向试验
我们经常使用卡方检验的另一种方法是询问两个分类变量是否相互关联。作为一个更现实的例子,让我们来考虑一个问题,当一个黑人司机被警察拦下时,他们是否比一个白人司机更有可能被搜查,斯坦福公开警务项目([https://open policing.stanford.edu/](https://openpolicing.stanford.edu/))研究了这个问题,并提供了我们可以用来分析问题的数据。我们将使用来自康涅狄格州的数据,因为它们相当小。首先清理这些数据,以删除所有不必要的数据(参见 code/process_ct_data.py)。
```r
# load police stop data
stopData <-
read_csv("data/CT_data_cleaned.csv") %>%
rename(searched = search_conducted)
```
表示分类分析数据的标准方法是通过 _ 列联表 _,列联表显示了属于每个变量值的每个可能组合的观测值的数量或比例。
让我们计算一下警察搜索数据的应急表:
```r
# compute and print two-way contingency table
summaryDf2way <-
stopData %>%
count(searched, driver_race) %>%
arrange(driver_race, searched)
summaryContingencyTable <-
summaryDf2way %>%
spread(driver_race, n)
pander(summaryContingencyTable)
```
<colgroup><col style="width: 15%"> <col style="width: 11%"> <col style="width: 11%"></colgroup>
| 已搜索 | 黑色 | 白色 |
| --- | --- | --- |
| 错误的 | 36244 个 | 239241 个 |
| 真的 | 1219 年 | 3108 个 |
使用比例而不是原始数字查看应急表也很有用,因为它们更容易在视觉上进行比较。
```r
# Compute and print contingency table using proportions
# rather than raw frequencies
summaryContingencyTableProportion <-
summaryContingencyTable %>%
mutate(
Black = Black / nrow(stopData), #count of Black individuals searched / total searched
White = White / nrow(stopData)
)
pander(summaryContingencyTableProportion, round = 4)
```
<colgroup><col style="width: 15%"> <col style="width: 12%"> <col style="width: 12%"></colgroup>
| searched | Black | White |
| --- | --- | --- |
| FALSE | 0.1295 年 | 0.855 个 |
| TRUE | 0.0044 美元 | 0.0111 个 |
Pearson 卡方检验允许我们检验观察到的频率是否与预期频率不同,因此我们需要确定如果搜索和种族不相关,我们期望在每个细胞中出现的频率,我们可以定义为 _ 独立。_ 请记住,如果 x 和 y 是独立的,那么:
![](img/ae00b41c7d126a5646eb9e06bf53e695.jpg)
也就是说,零独立假设下的联合概率仅仅是每个变量的 _ 边际 _ 概率的乘积。边际概率只是每一个事件发生的概率,与其他事件无关。我们可以计算这些边际概率,然后将它们相乘,得到独立状态下的预期比例。
| | 黑色 | 白色 | |
| --- | --- | --- | --- |
| 未搜索 | P(ns)*P(b) | P(ns)*P(w) | P(纳秒) |
| 已搜索 | P(S)*P(B) | P(S)*P(W) | P(S) |
| | P(B) | P(宽) | |
我们可以使用称为“外积”的线性代数技巧(通过`outer()`函数)来轻松计算。
```r
# first, compute the marginal probabilities
# probability of being each race
summaryDfRace <-
stopData %>%
count(driver_race) %>% #count the number of drivers of each race
mutate(
prop = n / sum(n) #compute the proportion of each race out of all drivers
)
# probability of being searched
summaryDfStop <-
stopData %>%
count(searched) %>% #count the number of searched vs. not searched
mutate(
prop = n / sum(n) # compute proportion of each outcome out all traffic stops
)
```
```r
# second, multiply outer product by n (all stops) to compute expected frequencies
expected <- outer(summaryDfRace$prop, summaryDfStop$prop) * nrow(stopData)
# create a data frame of expected frequencies for each race
expectedDf <-
data.frame(expected, driverRace = c("Black", "White")) %>%
rename(
NotSearched = X1,
Searched = X2
)
# tidy the data frame
expectedDfTidy <-
gather(expectedDf, searched, n, -driverRace) %>%
arrange(driverRace, searched)
```
```r
# third, add expected frequencies to the original summary table
# and fourth, compute the standardized squared difference between
# the observed and expected frequences
summaryDf2way <-
summaryDf2way %>%
mutate(expected = expectedDfTidy$n)
summaryDf2way <-
summaryDf2way %>%
mutate(stdSqDiff = (n - expected)**2 / expected)
pander(summaryDf2way)
```
<colgroup><col style="width: 15%"> <col style="width: 19%"> <col style="width: 12%"> <col style="width: 15%"> <col style="width: 15%"></colgroup>
| searched | 车手比赛 | N 号 | 预期 | 标准平方差 |
| --- | --- | --- | --- | --- |
| FALSE | 黑色 | 36244 | 36883.67 个 | 2009 年 11 月 |
| TRUE | Black | 1219 | 579.33 条 | 第 706.31 条 |
| FALSE | 白色 | 239241 | 238601.3 条 | 1.71 条 |
| TRUE | White | 3108 | 3747.67 美元 | 109.18 条 |
```r
# finally, compute chi-squared statistic by
# summing the standardized squared differences
chisq <- sum(summaryDf2way$stdSqDiff)
sprintf("Chi-squared value = %0.2f", chisq)
```
```r
## [1] "Chi-squared value = 828.30"
```
在计算了卡方统计之后,我们现在需要将其与卡方分布进行比较,以确定它与我们在无效假设下的期望相比有多极端。这种分布的自由度是![](img/206275e6a74b69f71d35c279b3e2d1c0.jpg)——因此,对于类似于这里的 2x2 表,![](img/ba8975aad6841ab43375c0463cf596e7.jpg)。这里的直觉是计算预期频率需要我们使用三个值:观察总数和两个变量的边际概率。因此,一旦计算出这些值,就只有一个数字可以自由变化,因此有一个自由度。鉴于此,我们可以计算卡方统计的 p 值:
```r
pval <- pchisq(chisq, df = 1, lower.tail = FALSE)
sprintf("p-value = %e", pval)
```
```r
## [1] "p-value = 3.795669e-182"
```
![](img/0a218a4dadec40dff9f8d2264a3c9d39.jpg)的 p 值非常小,表明如果种族和警察搜查之间真的没有关系,观察到的数据就不太可能,因此我们应该拒绝独立性的无效假设。
我们还可以使用 r 中的`chisq.test()`函数轻松执行此测试:
```r
# first need to rearrange the data into a 2x2 table
summaryDf2wayTable <-
summaryDf2way %>%
dplyr::select(-expected, -stdSqDiff) %>%
spread(searched, n) %>%
dplyr::select(-driver_race)
chisqTestResult <- chisq.test(summaryDf2wayTable, 1, correct = FALSE)
chisqTestResult
```
```r
##
## Pearson's Chi-squared test
##
## data: summaryDf2wayTable
## X-squared = 800, df = 1, p-value <2e-16
```
## 12.4 标准化残差
当我们发现卡方检验的显著效果时,这告诉我们,在无效假设下,数据是不可能的,但它并没有告诉我们 _ 数据有什么不同。为了更深入地了解数据与我们在零假设下预期的差异,我们可以检查模型的残差,该残差反映了数据(即观察到的频率)与每个单元中模型的偏差(即预期频率)。而不是查看原始残差(仅根据数据中观察的数量而变化),更常见的是查看其他 _ 标准化残差 _,其计算如下:_
![](img/fc6dfc834654c8a6d8b2323c1bd6b480.jpg)
其中![](img/31df9c730e19ca29b59dce64b99d98c1.jpg)和![](img/d8fdd0e28cfb03738fc5227885ee035a.jpg)分别是行和列的索引。我们可以为警察局的数据计算这些数据:
```r
# compute standardized residuals
summaryDf2way <-
summaryDf2way %>%
mutate(stdRes = (n - expected)/sqrt(expected))
pander(summaryDf2way)
```
<colgroup><col style="width: 15%"> <col style="width: 19%"> <col style="width: 12%"> <col style="width: 15%"> <col style="width: 16%"> <col style="width: 11%"></colgroup>
| 已搜索 | 车手比赛 | N 号 | 预期 | 标准平方差 | 标准普尔 |
| --- | --- | --- | --- | --- | --- |
| 错误的 | 黑色 | 36244 个 | 36883.67 个 | 2009 年 11 月 | -3.33 条 |
| 真的 | Black | 1219 年 | 579.33 条 | 第 706.31 条 | 26.58 美元 |
| FALSE | 白色 | 239241 个 | 238601.3 条 | 1.71 条 | 1.31 条 |
| TRUE | White | 3108 个 | 3747.67 美元 | 109.18 条 | -10.45 美元 |
这些标准化的残差可以解释为 z 分数——在这种情况下,我们看到,基于独立性,对黑人个体的搜索次数大大高于预期,而对白人个体的搜索次数大大低于预期。这为我们提供了解释显著的卡方结果所需的上下文。
## 12.5 优势比
我们还可以使用前面介绍的优势比在应急表中表示不同结果的相对可能性,以便更好地了解影响的大小。首先,我们表示每一场比赛被阻止的几率:
![](img/e6acd79788a8dbe24193b20151db6a12.jpg)
![](img/9aa2622d4176179ddfde24206732f146.jpg)![](img/be906c2bfa49b7ba6d831a99b9f47c71.jpg)
根据这一数据集,优势比显示,黑人和白人驾驶员被搜索的几率要高出 2.59 倍。
## 12.6 贝叶斯系数
![](img/ed955fccd67d36cab5ec92cfcd8c1742.jpg)
我们在关于贝叶斯统计的前一章中讨论了贝叶斯因子——你可能记得它代表了两种假设下数据的可能性比:贝叶斯因子在某种程度上类似于 p 值和影响大小,也就是说它们的解释是某种意义上的。t 主观。他们的解释有各种各样的指导方针——这是来自 Kass 和 Raftery(1995)的一个:
| 高炉 | 解释 |
| --- | --- |
| 1 到 3 | 不值得一提 |
| 3 至 20 | 积极的 |
| 20 至 150 | 坚强的 |
| 150 及以上 | 非常强壮 |
我们可以使用 BayesFactor 包中的`contingencyTableBF()`函数计算警察搜索数据的 Bayes 因子:
```r
# compute Bayes factor
# using independent multinomial sampling plan in which row totals (driver race)
# are fixed
bf <-
contingencyTableBF(as.matrix(summaryDf2wayTable),
sampleType = "indepMulti",
fixedMargin = "cols"
)
bf
```
```r
## Bayes factor analysis
## --------------
## [1] Non-indep. (a=1) : 1.8e+142 ±0%
##
## Against denominator:
## Null, independence, a = 1
## ---
## Bayes factor type: BFcontingencyTable, independent multinomial
```
这表明,在这个数据集中,支持驾驶员种族和警察搜索之间关系的证据非常有力。
## 12.7 超出 2 x 2 表的分类分析
分类分析也可以应用于应急表,其中每个变量有两个以上的类别。
例如,让我们看一下 nhanes 的数据,比较变量 _depressed_,它表示“参与者感到沮丧、沮丧或绝望的自我报告天数”。此变量编码为`None``Several``Most`。让我们测试这个变量是否与 _sleeptrouble_ 变量相关,这个变量报告个人是否向医生报告了睡眠问题。
```r
# summarize depression as a function of sleep trouble
depressedSleepTrouble <-
NHANES_adult %>%
drop_na(SleepTrouble, Depressed) %>%
count(SleepTrouble, Depressed) %>%
arrange(SleepTrouble, Depressed)
depressedSleepTroubleTable <-
depressedSleepTrouble %>%
spread(SleepTrouble, n) %>%
rename(
NoSleepTrouble = No,
YesSleepTrouble = Yes
)
pander(depressedSleepTroubleTable)
```
<colgroup><col style="width: 16%"> <col style="width: 23%"> <col style="width: 23%"></colgroup>
| 沮丧的 | 无阻力 | 是的,可重复 |
| --- | --- | --- |
| 无 | 2614 个 | 676 个 |
| 几个 | 418 个 | 249 个 |
| 大多数 | 138 个 | 145 个 |
简单地看一下这些数据,我们就可以知道这两个变量之间可能存在某种关系;特别是,尽管睡眠问题患者的总数比没有睡眠问题的患者要少很多,但对于大多数时间都处于抑郁状态的患者来说,睡眠问题患者的数量更大。比没有的要水。我们可以直接使用卡方检验对其进行量化;如果我们的数据框只包含两个变量,那么我们可以简单地将数据框作为参数提供给`chisq.test()`函数:
```r
# need to remove the column with the label names
depressedSleepTroubleTable <-
depressedSleepTroubleTable %>%
dplyr::select(-Depressed)
depressedSleepChisq <- chisq.test(depressedSleepTroubleTable)
depressedSleepChisq
```
```r
##
## Pearson's Chi-squared test
##
## data: depressedSleepTroubleTable
## X-squared = 200, df = 2, p-value <2e-16
```
这项测试表明,抑郁和睡眠问题之间有很强的关系。我们还可以计算贝叶斯因子来量化有利于替代假设的证据的强度:
```r
# compute bayes factor, using a joint multinomial sampling plan
bf <-
contingencyTableBF(
as.matrix(depressedSleepTroubleTable),
sampleType = "jointMulti"
)
bf
```
```r
## Bayes factor analysis
## --------------
## [1] Non-indep. (a=1) : 1.8e+35 ±0%
##
## Against denominator:
## Null, independence, a = 1
## ---
## Bayes factor type: BFcontingencyTable, joint multinomial
```
在这里,我们看到贝叶斯系数非常大,这表明支持抑郁和睡眠问题之间关系的证据非常有力。
## 12.8 注意辛普森悖论
上述应急表是对大量观察结果的总结,但有时会产生误导。让我们以棒球为例。下表显示了 1995-1997 年间德里克·杰特和大卫·贾斯汀的击球数据(击数/击数和平均击球数):
| 玩家 | 1995 年 | | 1996 年 | | 1997 年 | | 合并 | |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| 基特 | 12 月 48 日 | 0.250 | 183/582 年 | .314 条 | 190/654 年 | .291 条 | 385/1284 年 | **.300** |
| 大卫·正义 | 104/411 号 | **.253** | 45/140 分 | **.321** | 163/495 年 | **.329** | 312/1046 年 | .298 条 |
如果你仔细观察,你会发现有些奇怪的事情正在发生:在每一年,正义比杰特有一个更高的击球平均值,但当我们结合所有三年的数据,杰特的平均值实际上高于正义!这是一个被称为 _ 辛普森悖论 _ 的现象的例子,在这种现象中,组合数据集中的模式可能不存在于数据的任何子集中。当有另一个变量可能在不同的子集之间发生变化时,就会发生这种情况——在这种情况下,AT 蝙蝠的数量随着时间的推移而变化,1995 年司法部的击球次数更多(击球平均数较低时)。我们把它称为一个潜伏变量(htg2),每当我们检查分类数据时,注意这些变量总是很重要的。
\ No newline at end of file
# 12 分类关系建模
到目前为止,我们已经讨论了统计建模和假设检验的一般概念,并将它们应用到一些简单的分析中。在本章中,我们将重点介绍 _ 分类 _ 关系的建模,通过建模,我们可以表示在名义(有时是序数)尺度上测量的变量之间的关系。这些数据通常用计数来表示;也就是说,对于变量的每个值(或多个变量的值的组合),有多少观测值取这个值?例如,当我们计算每个专业有多少人在我们班上时,我们将根据数据拟合一个分类模型。
## 12.1 示例:糖果颜色
假设我买了一袋 100 颗糖果,上面写着有 1/3 巧克力、1/3 甘草和 1/3 胶球。当我数包里的糖果时,我们得到以下数字:
```r
candyDf <-
tibble(
candyType = c("chocolate", "licorice", "gumball"),
count = c(30, 33, 37)
)
pander(candyDf)
```
<colgroup><col style="width: 16%"> <col style="width: 9%"></colgroup>
| 糖果型 | 计数 |
| --- | --- |
| 巧克力 | 30 个 |
| 甘草 | 33 |
| 胶球 | 37 岁 |
因为我更喜欢巧克力,而不是甘草或树胶球,我觉得有点被撕了。我想知道的是:如果每种糖果的真实概率是每种糖果平均 1/3 的比例,那么这种计数的可能性有多大?
## 12.2 皮尔逊卡方检验
皮尔逊卡方检验为我们提供了一种方法来检验观察到的计数数据是否与定义零假设的某些特定预期值不同:
![](img/a39dfdf63e525b1d15f99c58b2ae7a2e.jpg)
在我们的糖果例子中,无效假设是每种糖果的比例是相等的。我们可以计算我们观察到的糖果计数的卡方统计,如下所示:
```r
# compute chi-squared statistic
nullExpectation <- c(1 / 3, 1 / 3, 1 / 3) * sum(candyDf$count)
chisqVal <-
sum(
((candyDf$count - nullExpectation)**2) / nullExpectation
)
```
这个分析的卡方统计结果是 0.74,这本身是不可解释的,因为它取决于加在一起的不同值的数量。然而,我们可以利用这样一个事实:卡方统计量是根据零假设下的特定分布分布分布的,这就是所谓的 _ 卡方 _ 分布。这种分布被定义为一组标准正态随机变量的平方和;它有若干自由度,等于被加在一起的变量数。分布的形状取决于自由度的数量。图[12.1](#fig:chisqDist)显示了几种不同自由度的分布示例。
![Examples of the chi-squared distribution for various degrees of freedom.](img/file74.png)
图 12.1 不同自由度的卡方分布示例。
让我们通过模拟来验证卡方分布是否准确地描述了一组标准正态随机变量的平方和。
```r
# simulate 50,000 sums of 8 standard normal random variables and compare
# to theoretical chi-squared distribution
# create a matrix with 50k columns of 8 rows of squared normal random variables
d <- replicate(50000, rnorm(n = 8, mean = 0, sd = 1)**2)
# sum each column of 8 variables
dMean <- apply(d, 2, sum)
# create a data frame of the theoretical chi-square distribution
# with 8 degrees of freedom
csDf <-
data.frame(x = seq(0.01, 30, 0.01)) %>%
mutate(chisq = dchisq(x, 8))
```
[12.2](#fig:chisqSim)显示,理论分布与重复将一组随机正态变量的平方相加的模拟结果非常吻合。
![Simulation of sum of squared random normal variables. The histogram is based on the sum of squares of 50,000 sets of 8 random normal variables; the blue line shows the values of the theoretical chi-squared distribution with 8 degrees of freedom.](img/file75.png)
图 12.2 平方随机正态变量和的模拟。柱状图是基于 5 万组 8 个随机正态变量的平方和;蓝线显示了 8 个自由度下理论卡方分布的值。
对于糖果的例子,我们可以计算在所有糖果的相同频率的零假设下,我们观察到的卡方值为 0.74 的可能性。我们使用自由度等于 k-1(其中 k=类别数)的卡方分布,因为我们在计算平均值以生成预期值时失去了一个自由度。
```r
pval <- pchisq(chisqVal, df = 2, lower.tail = FALSE) #df = degrees of freedom
sprintf("p-value = %0.3f", pval)
```
```r
## [1] "p-value = 0.691"
```
这表明,观察到的糖果数量并不是特别令人惊讶的,基于印刷在糖果袋上的比例,我们不会拒绝等比的无效假设。
## 12.3 应急表及双向试验
我们经常使用卡方检验的另一种方法是询问两个分类变量是否相互关联。作为一个更现实的例子,让我们来考虑一个问题,当一个黑人司机被警察拦下时,他们是否比一个白人司机更有可能被搜查,斯坦福公开警务项目([https://open policing.stanford.edu/](https://openpolicing.stanford.edu/))研究了这个问题,并提供了我们可以用来分析问题的数据。我们将使用来自康涅狄格州的数据,因为它们相当小。首先清理这些数据,以删除所有不必要的数据(参见 code/process_ct_data.py)。
```r
# load police stop data
stopData <-
read_csv("data/CT_data_cleaned.csv") %>%
rename(searched = search_conducted)
```
表示分类分析数据的标准方法是通过 _ 列联表 _,列联表显示了属于每个变量值的每个可能组合的观测值的数量或比例。
让我们计算一下警察搜索数据的应急表:
```r
# compute and print two-way contingency table
summaryDf2way <-
stopData %>%
count(searched, driver_race) %>%
arrange(driver_race, searched)
summaryContingencyTable <-
summaryDf2way %>%
spread(driver_race, n)
pander(summaryContingencyTable)
```
<colgroup><col style="width: 15%"> <col style="width: 11%"> <col style="width: 11%"></colgroup>
| 已搜索 | 黑色 | 白色 |
| --- | --- | --- |
| 错误的 | 36244 个 | 239241 个 |
| 真的 | 1219 年 | 3108 个 |
使用比例而不是原始数字查看应急表也很有用,因为它们更容易在视觉上进行比较。
```r
# Compute and print contingency table using proportions
# rather than raw frequencies
summaryContingencyTableProportion <-
summaryContingencyTable %>%
mutate(
Black = Black / nrow(stopData), #count of Black individuals searched / total searched
White = White / nrow(stopData)
)
pander(summaryContingencyTableProportion, round = 4)
```
<colgroup><col style="width: 15%"> <col style="width: 12%"> <col style="width: 12%"></colgroup>
| searched | Black | White |
| --- | --- | --- |
| FALSE | 0.1295 年 | 0.855 个 |
| TRUE | 0.0044 美元 | 0.0111 个 |
Pearson 卡方检验允许我们检验观察到的频率是否与预期频率不同,因此我们需要确定如果搜索和种族不相关,我们期望在每个细胞中出现的频率,我们可以定义为 _ 独立。_ 请记住,如果 x 和 y 是独立的,那么:
![](img/ae00b41c7d126a5646eb9e06bf53e695.jpg)
也就是说,零独立假设下的联合概率仅仅是每个变量的 _ 边际 _ 概率的乘积。边际概率只是每一个事件发生的概率,与其他事件无关。我们可以计算这些边际概率,然后将它们相乘,得到独立状态下的预期比例。
| | 黑色 | 白色 | |
| --- | --- | --- | --- |
| 未搜索 | P(ns)*P(b) | P(ns)*P(w) | P(纳秒) |
| 已搜索 | P(S)*P(B) | P(S)*P(W) | P(S) |
| | P(B) | P(宽) | |
我们可以使用称为“外积”的线性代数技巧(通过`outer()`函数)来轻松计算。
```r
# first, compute the marginal probabilities
# probability of being each race
summaryDfRace <-
stopData %>%
count(driver_race) %>% #count the number of drivers of each race
mutate(
prop = n / sum(n) #compute the proportion of each race out of all drivers
)
# probability of being searched
summaryDfStop <-
stopData %>%
count(searched) %>% #count the number of searched vs. not searched
mutate(
prop = n / sum(n) # compute proportion of each outcome out all traffic stops
)
```
```r
# second, multiply outer product by n (all stops) to compute expected frequencies
expected <- outer(summaryDfRace$prop, summaryDfStop$prop) * nrow(stopData)
# create a data frame of expected frequencies for each race
expectedDf <-
data.frame(expected, driverRace = c("Black", "White")) %>%
rename(
NotSearched = X1,
Searched = X2
)
# tidy the data frame
expectedDfTidy <-
gather(expectedDf, searched, n, -driverRace) %>%
arrange(driverRace, searched)
```
```r
# third, add expected frequencies to the original summary table
# and fourth, compute the standardized squared difference between
# the observed and expected frequences
summaryDf2way <-
summaryDf2way %>%
mutate(expected = expectedDfTidy$n)
summaryDf2way <-
summaryDf2way %>%
mutate(stdSqDiff = (n - expected)**2 / expected)
pander(summaryDf2way)
```
<colgroup><col style="width: 15%"> <col style="width: 19%"> <col style="width: 12%"> <col style="width: 15%"> <col style="width: 15%"></colgroup>
| searched | 车手比赛 | N 号 | 预期 | 标准平方差 |
| --- | --- | --- | --- | --- |
| FALSE | 黑色 | 36244 | 36883.67 个 | 2009 年 11 月 |
| TRUE | Black | 1219 | 579.33 条 | 第 706.31 条 |
| FALSE | 白色 | 239241 | 238601.3 条 | 1.71 条 |
| TRUE | White | 3108 | 3747.67 美元 | 109.18 条 |
```r
# finally, compute chi-squared statistic by
# summing the standardized squared differences
chisq <- sum(summaryDf2way$stdSqDiff)
sprintf("Chi-squared value = %0.2f", chisq)
```
```r
## [1] "Chi-squared value = 828.30"
```
在计算了卡方统计之后,我们现在需要将其与卡方分布进行比较,以确定它与我们在无效假设下的期望相比有多极端。这种分布的自由度是![](img/206275e6a74b69f71d35c279b3e2d1c0.jpg)——因此,对于类似于这里的 2x2 表,![](img/ba8975aad6841ab43375c0463cf596e7.jpg)。这里的直觉是计算预期频率需要我们使用三个值:观察总数和两个变量的边际概率。因此,一旦计算出这些值,就只有一个数字可以自由变化,因此有一个自由度。鉴于此,我们可以计算卡方统计的 p 值:
```r
pval <- pchisq(chisq, df = 1, lower.tail = FALSE)
sprintf("p-value = %e", pval)
```
```r
## [1] "p-value = 3.795669e-182"
```
![](img/0a218a4dadec40dff9f8d2264a3c9d39.jpg)的 p 值非常小,表明如果种族和警察搜查之间真的没有关系,观察到的数据就不太可能,因此我们应该拒绝独立性的无效假设。
我们还可以使用 r 中的`chisq.test()`函数轻松执行此测试:
```r
# first need to rearrange the data into a 2x2 table
summaryDf2wayTable <-
summaryDf2way %>%
dplyr::select(-expected, -stdSqDiff) %>%
spread(searched, n) %>%
dplyr::select(-driver_race)
chisqTestResult <- chisq.test(summaryDf2wayTable, 1, correct = FALSE)
chisqTestResult
```
```r
##
## Pearson's Chi-squared test
##
## data: summaryDf2wayTable
## X-squared = 800, df = 1, p-value <2e-16
```
## 12.4 标准化残差
当我们发现卡方检验的显著效果时,这告诉我们,在无效假设下,数据是不可能的,但它并没有告诉我们 _ 数据有什么不同。为了更深入地了解数据与我们在零假设下预期的差异,我们可以检查模型的残差,该残差反映了数据(即观察到的频率)与每个单元中模型的偏差(即预期频率)。而不是查看原始残差(仅根据数据中观察的数量而变化),更常见的是查看其他 _ 标准化残差 _,其计算如下:_
![](img/fc6dfc834654c8a6d8b2323c1bd6b480.jpg)
其中![](img/31df9c730e19ca29b59dce64b99d98c1.jpg)和![](img/d8fdd0e28cfb03738fc5227885ee035a.jpg)分别是行和列的索引。我们可以为警察局的数据计算这些数据:
```r
# compute standardized residuals
summaryDf2way <-
summaryDf2way %>%
mutate(stdRes = (n - expected)/sqrt(expected))
pander(summaryDf2way)
```
<colgroup><col style="width: 15%"> <col style="width: 19%"> <col style="width: 12%"> <col style="width: 15%"> <col style="width: 16%"> <col style="width: 11%"></colgroup>
| 已搜索 | 车手比赛 | N 号 | 预期 | 标准平方差 | 标准普尔 |
| --- | --- | --- | --- | --- | --- |
| 错误的 | 黑色 | 36244 个 | 36883.67 个 | 2009 年 11 月 | -3.33 条 |
| 真的 | Black | 1219 年 | 579.33 条 | 第 706.31 条 | 26.58 美元 |
| FALSE | 白色 | 239241 个 | 238601.3 条 | 1.71 条 | 1.31 条 |
| TRUE | White | 3108 个 | 3747.67 美元 | 109.18 条 | -10.45 美元 |
这些标准化的残差可以解释为 z 分数——在这种情况下,我们看到,基于独立性,对黑人个体的搜索次数大大高于预期,而对白人个体的搜索次数大大低于预期。这为我们提供了解释显著的卡方结果所需的上下文。
## 12.5 优势比
我们还可以使用前面介绍的优势比在应急表中表示不同结果的相对可能性,以便更好地了解影响的大小。首先,我们表示每一场比赛被阻止的几率:
![](img/e6acd79788a8dbe24193b20151db6a12.jpg)
![](img/9aa2622d4176179ddfde24206732f146.jpg)![](img/be906c2bfa49b7ba6d831a99b9f47c71.jpg)
根据这一数据集,优势比显示,黑人和白人驾驶员被搜索的几率要高出 2.59 倍。
## 12.6 贝叶斯系数
![](img/ed955fccd67d36cab5ec92cfcd8c1742.jpg)
我们在关于贝叶斯统计的前一章中讨论了贝叶斯因子——你可能记得它代表了两种假设下数据的可能性比:贝叶斯因子在某种程度上类似于 p 值和影响大小,也就是说它们的解释是某种意义上的。t 主观。他们的解释有各种各样的指导方针——这是来自 Kass 和 Raftery(1995)的一个:
| 高炉 | 解释 |
| --- | --- |
| 1 到 3 | 不值得一提 |
| 3 至 20 | 积极的 |
| 20 至 150 | 坚强的 |
| 150 及以上 | 非常强壮 |
我们可以使用 BayesFactor 包中的`contingencyTableBF()`函数计算警察搜索数据的 Bayes 因子:
```r
# compute Bayes factor
# using independent multinomial sampling plan in which row totals (driver race)
# are fixed
bf <-
contingencyTableBF(as.matrix(summaryDf2wayTable),
sampleType = "indepMulti",
fixedMargin = "cols"
)
bf
```
```r
## Bayes factor analysis
## --------------
## [1] Non-indep. (a=1) : 1.8e+142 ±0%
##
## Against denominator:
## Null, independence, a = 1
## ---
## Bayes factor type: BFcontingencyTable, independent multinomial
```
这表明,在这个数据集中,支持驾驶员种族和警察搜索之间关系的证据非常有力。
## 12.7 超出 2 x 2 表的分类分析
分类分析也可以应用于应急表,其中每个变量有两个以上的类别。
例如,让我们看一下 nhanes 的数据,比较变量 _depressed_,它表示“参与者感到沮丧、沮丧或绝望的自我报告天数”。此变量编码为`None``Several``Most`。让我们测试这个变量是否与 _sleeptrouble_ 变量相关,这个变量报告个人是否向医生报告了睡眠问题。
```r
# summarize depression as a function of sleep trouble
depressedSleepTrouble <-
NHANES_adult %>%
drop_na(SleepTrouble, Depressed) %>%
count(SleepTrouble, Depressed) %>%
arrange(SleepTrouble, Depressed)
depressedSleepTroubleTable <-
depressedSleepTrouble %>%
spread(SleepTrouble, n) %>%
rename(
NoSleepTrouble = No,
YesSleepTrouble = Yes
)
pander(depressedSleepTroubleTable)
```
<colgroup><col style="width: 16%"> <col style="width: 23%"> <col style="width: 23%"></colgroup>
| 沮丧的 | 无阻力 | 是的,可重复 |
| --- | --- | --- |
| 无 | 2614 个 | 676 个 |
| 几个 | 418 个 | 249 个 |
| 大多数 | 138 个 | 145 个 |
简单地看一下这些数据,我们就可以知道这两个变量之间可能存在某种关系;特别是,尽管睡眠问题患者的总数比没有睡眠问题的患者要少很多,但对于大多数时间都处于抑郁状态的患者来说,睡眠问题患者的数量更大。比没有的要水。我们可以直接使用卡方检验对其进行量化;如果我们的数据框只包含两个变量,那么我们可以简单地将数据框作为参数提供给`chisq.test()`函数:
```r
# need to remove the column with the label names
depressedSleepTroubleTable <-
depressedSleepTroubleTable %>%
dplyr::select(-Depressed)
depressedSleepChisq <- chisq.test(depressedSleepTroubleTable)
depressedSleepChisq
```
```r
##
## Pearson's Chi-squared test
##
## data: depressedSleepTroubleTable
## X-squared = 200, df = 2, p-value <2e-16
```
这项测试表明,抑郁和睡眠问题之间有很强的关系。我们还可以计算贝叶斯因子来量化有利于替代假设的证据的强度:
```r
# compute bayes factor, using a joint multinomial sampling plan
bf <-
contingencyTableBF(
as.matrix(depressedSleepTroubleTable),
sampleType = "jointMulti"
)
bf
```
```r
## Bayes factor analysis
## --------------
## [1] Non-indep. (a=1) : 1.8e+35 ±0%
##
## Against denominator:
## Null, independence, a = 1
## ---
## Bayes factor type: BFcontingencyTable, joint multinomial
```
在这里,我们看到贝叶斯系数非常大,这表明支持抑郁和睡眠问题之间关系的证据非常有力。
## 12.8 注意辛普森悖论
上述应急表是对大量观察结果的总结,但有时会产生误导。让我们以棒球为例。下表显示了 1995-1997 年间德里克·杰特和大卫·贾斯汀的击球数据(击数/击数和平均击球数):
| 玩家 | 1995 年 | | 1996 年 | | 1997 年 | | 合并 | |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| 基特 | 12 月 48 日 | 0.250 | 183/582 年 | .314 条 | 190/654 年 | .291 条 | 385/1284 年 | **.300** |
| 大卫·正义 | 104/411 号 | **.253** | 45/140 分 | **.321** | 163/495 年 | **.329** | 312/1046 年 | .298 条 |
如果你仔细观察,你会发现有些奇怪的事情正在发生:在每一年,正义比杰特有一个更高的击球平均值,但当我们结合所有三年的数据,杰特的平均值实际上高于正义!这是一个被称为 _ 辛普森悖论 _ 的现象的例子,在这种现象中,组合数据集中的模式可能不存在于数据的任何子集中。当有另一个变量可能在不同的子集之间发生变化时,就会发生这种情况——在这种情况下,AT 蝙蝠的数量随着时间的推移而变化,1995 年司法部的击球次数更多(击球平均数较低时)。我们把它称为一个潜伏变量(htg2),每当我们检查分类数据时,注意这些变量总是很重要的。
\ No newline at end of file
# 13 建模持续关系
大多数人都熟悉 _ 相关 _ 的概念,在本章中,我们将对这个常用和误解的概念提供更正式的理解。
\ No newline at end of file
大多数人都熟悉 _ 相关 _ 的概念,在本章中,我们将对这个常用和误解的概念提供更正式的理解。
## 13.1 一个例子:仇恨犯罪和收入不平等
2017 年,网站 fivethirtyeight.com 发表了一篇题为“仇恨犯罪率上升与收入不平等(htg1)”的文章,讨论了 2016 年总统选举后仇恨犯罪率与收入不平等之间的关系。报道分析了来自联邦调查局和南方贫困法中心的仇恨犯罪数据,并据此报告:
> “我们发现,收入不平等是全美国调整人口仇恨犯罪和仇恨事件的最重要决定因素”。
此分析的数据包含在`fivethirtyeight`r 包中,这使得我们很容易访问它们。报道中的分析集中在收入不平等(定义为一个叫做 _ 基尼指数 _ 的量)与各州仇恨犯罪流行率之间的关系。
#
## 13.1.1 量化不平等:基尼指数
在我们查看报道中的分析之前,首先要了解如何使用基尼指数来量化不平等。基尼指数通常用一条曲线来定义,这条曲线描述了收入与收入水平等于或小于该水平的人口比例之间的关系,称为 _ 洛伦兹曲线 _。然而,另一种更直观的思考方式是:收入之间的相对平均绝对差异除以二(摘自[https://en.wikipedia.org/wiki/gini_coefficient](https://en.wikipedia.org/wiki/Gini_coefficient)):
![](img/1283ad31298b6a845d0b54ea5116c77d.jpg)
![Lorenz curves for A) perfect equality, B) normally distributed income, and C) high inequality (equal income except for one very wealthy individual).](img/file76.png)
图 13.1 洛伦兹曲线表示 a)完全平等,b)正态分布收入,c)高度不平等(除一个非常富有的个人外,收入相等)。
[13.1](#fig:gini0)显示了几种不同收入分配的洛伦兹曲线。左上方的面板(A)显示了一个例子,其中有 10 个人,每个人的收入完全相同。两个点之间的间隔长度相等,表明每个人在总收入中所占的份额相同。右上角的面板(B)显示了一个收入正态分布的例子。左下角的面板显示了一个不平等程度很高的例子:每个人的收入都是平等的(40000 美元),只有一个人的收入是 40000000 美元。根据美国人口普查,2010 年美国的基尼指数为 0.469,大约在我们的正态分布和最大不相等的例子之间下降了一半。
## 13.2 收入不平等是否与仇恨犯罪有关?
现在我们了解了基尼指数,我们可以看看收入不平等与仇恨犯罪率之间的关系(见图[13.2](#fig:hateCrimeGini))。
```r
hateCrimes <-
hate_crimes %>%
mutate(state_abb = state.abb[match(state,state.name)]) %>%
drop_na(avg_hatecrimes_per_100k_fbi)
hateCrimes$state_abb[hateCrimes$state=="District of Columbia"]='DC'
ggplot(hateCrimes,aes(gini_index,avg_hatecrimes_per_100k_fbi,label=state_abb)) +
geom_point() +
geom_text(aes(label=state_abb),hjust=0, vjust=0) +
theme(plot.title = element_text(size = 20, face = "bold")) +
xlab('Gini index') +
ylab('Avg hate crimes per 100K population (FBI)') +
theme(plot.margin = unit(c(1,1,1,1), "cm"))
```
![Plot of rates of hate crimes vs. Gini index.](img/file77.png)
图 13.2 仇恨犯罪率与基尼指数的关系图。
从数据来看,这两个变量之间似乎有一个正的关系。我们如何量化这种关系?
## 13.3 协方差和相关性
量化两个变量之间关系的一种方法是 _ 协方差 _。记住,单个变量的方差计算如下:
![](img/134e6de0f245b57b2f07e6a35831cb9e.jpg)
这告诉我们每个观察值与平均值相差多远。协方差告诉我们两个不同的变量在观测值之间的偏差是否存在关系。定义如下:
![](img/b268259e45c69b4bfdac3da4ed6962b5.jpg)
当 x 和 y 都高度偏离平均值时,该值将远离零;如果它们在同一方向上偏离,则协方差为正,而如果它们在相反方向上偏离,则协方差为负。让我们先看一个玩具的例子。
```r
# create data for toy example of covariance
df <-
tibble(x = c(3, 5, 8, 10, 12)) %>%
mutate(y = x + round(rnorm(n = 5, mean = 0, sd = 2))) %>%
mutate(
y_dev = y - mean(y),
x_dev = x - mean(x)
) %>%
mutate(crossproduct = y_dev * x_dev)
pander(df)
```
<colgroup><col style="width: 6%"> <col style="width: 6%"> <col style="width: 11%"> <col style="width: 11%"> <col style="width: 19%"></colgroup>
| X | 是 | Y 轴偏差 | X 轴偏差 | 叉乘 |
| --- | --- | --- | --- | --- |
| 三 | 1 个 | -6.6 条 | -4.6 节 | 30.36 天 |
| 5 个 | 3 | -4.6 | -第 2.6 条 | 11.96 年 |
| 8 个 | 8 | 0.4 倍 | 0.4 | 0.16 分 |
| 10 个 | 12 个 | 第 4.4 条 | 第 2.4 条 | 10.56 条 |
| 12 | 14 | 第 6.4 条 | 4.4 | 28.16 条 |
```r
# compute covariance
sprintf("sum of cross products = %.2f", sum(df$crossproduct))
```
```r
## [1] "sum of cross products = 81.20"
```
```r
covXY <- sum(df$crossproduct) / (nrow(df) - 1)
sprintf("covariance: %.2f", covXY)
```
```r
## [1] "covariance: 20.30"
```
我们通常不使用协方差来描述变量之间的关系,因为它随数据的总体方差水平而变化。相反,我们通常使用 _ 相关系数 _(通常在统计学家 Karl Pearson 之后称为 _Pearson 相关 _)。通过用两个变量的标准偏差缩放协方差来计算相关性:
![](img/69240baf835a9a8f0311c744556cb85f.jpg)
```r
# compute the correlation coefficient
corXY <- sum(df$crossproduct) / ((nrow(df) - 1) * sd(df$x) * sd(df$y))
sprintf("correlation coefficient = %.2f", corXY)
```
```r
## [1] "correlation coefficient = 0.99"
```
我们还可以使用 r 中的`cor()`函数轻松计算相关值:
```r
# compute r using built-in function
c <- cor(df$x, df$y)
sprintf("correlation coefficient = %.2f", c)
```
```r
## [1] "correlation coefficient = 0.99"
```
相关系数是有用的,因为它在-1 和 1 之间变化,不管数据的性质如何-事实上,我们在讨论影响大小时已经讨论过相关系数。正如我们在上一章关于影响大小的内容中看到的,1 的相关性表示一个完美的线性关系,-1 的相关性表示一个完美的负关系,0 的相关性表示没有线性关系。
我们可以计算仇恨犯罪数据的相关系数:
```r
corGiniHC <-
cor(
hateCrimes$gini_index,
hateCrimes$avg_hatecrimes_per_100k_fbi
)
sprintf('correlation coefficient = %.2f',corGiniHC)
```
```r
## [1] "correlation coefficient = 0.42"
```
#
## 13.3.1 相关性假设检验
相关值 0.42 似乎表明两个变量之间的关系相当强,但我们也可以想象,即使没有关系,这种情况也可能是偶然发生的。我们可以使用一个简单的公式来测试相关性为零的空假设,该公式允许我们将相关性值转换为 _t_ 统计:
![](img/8a7e337cf588bdff49bc2f192c7eac97.jpg)
在零假设![](img/fed676a693d5c98216a52108307827c6.jpg)下,该统计量以自由度为 t 分布。我们可以使用 r 中的`cor.test()`函数计算:
```r
# perform correlation test on hate crime data
cor.test(
hateCrimes$avg_hatecrimes_per_100k_fbi,
hateCrimes$gini_index
)
```
```r
##
## Pearson's product-moment correlation
##
## data: hateCrimes$avg_hatecrimes_per_100k_fbi and hateCrimes$gini_index
## t = 3, df = 50, p-value = 0.002
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.16 0.63
## sample estimates:
## cor
## 0.42
```
这个测试表明,R 值的可能性很低,这个极限或更高,所以我们将拒绝![](img/6937cddf1a8c4f47df75a6977e1234d1.jpg)的无效假设。注意,这个测试假设两个变量都是正态分布的。
我们也可以通过随机化来检验这一点,在随机化中,我们重复地改变其中一个变量的值并计算相关性,然后将我们观察到的相关性值与这个零分布进行比较,以确定我们观察到的值在零假设下的可能性。结果如图[13.3](#fig:shuffleCorr)所示。使用随机化计算的 p 值与 t 检验给出的答案相当相似。
```r
# compute null distribution by shuffling order of variable values
# create a function to compute the correlation on the shuffled values
shuffleCorr <- function(x, y) {
xShuffled <- sample(x)
return(cor(xShuffled, y))
}
# run this function 2500 times
shuffleDist <-
replicate(
2500,
shuffleCorr(hateCrimes$avg_hatecrimes_per_100k_fbi, hateCrimes$gini_index)
)
```
![Histogram of correlation values under the null hypothesis, obtained by shuffling values. Observed value is denoted by blue line.](img/file78.png)
图 13.3 零假设下相关值的柱状图,通过改变值获得。观测值用蓝线表示。
#
## 13.3.2 稳健相关性
在图[13.2](#fig:hateCrimeGini)中,您可能注意到了一些有点奇怪的地方——其中一个数据点(哥伦比亚特区的数据点)似乎与其他数据点非常不同。我们称之为 _ 离群值 _,标准相关系数对离群值非常敏感。例如,在图[13.4](#fig:outlierCorr)中,我们可以看到一个孤立的数据点是如何导致非常高的正相关值的,即使其他数据点之间的实际关系是完全负的。
![An simulated example of the effects of outliers on correlation. Without the outlier the remainder of the datapoints have a perfect negative correlation, but the single outlier changes the correlation value to highly positive.](img/file79.png)
图 13.4 异常值对相关性影响的模拟示例。如果没有离群值,其余数据点具有完全的负相关,但单个离群值将相关值更改为高度正相关。
解决离群值问题的一种方法是在排序后,在数据的列组上计算相关性,而不是在数据本身上计算相关性;这被称为 _ 斯皮尔曼相关性 _。图[13.4](#fig:outlierCorr)中的 Pearson 相关性为 0.83,而 Spearman 相关性为-0.45,表明等级相关性降低了异常值的影响。
我们可以使用`cor.test`函数计算仇恨犯罪数据的等级相关性:
```r
corTestSpearman <- cor.test( hateCrimes$avg_hatecrimes_per_100k_fbi,
hateCrimes$gini_index,
method = "spearman")
corTestSpearman
```
```r
##
## Spearman's rank correlation rho
##
## data: hateCrimes$avg_hatecrimes_per_100k_fbi and hateCrimes$gini_index
## S = 20000, p-value = 0.8
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.033
```
现在我们看到相关性不再显著(实际上接近于零),这表明 Fivethirtyeight 博客帖子的声明可能由于离群值的影响而不正确。
#
## 13.3.3 贝叶斯相关分析
我们也可以使用贝叶斯分析来分析五个第八个数据,这有两个优点。首先,它为我们提供了一个后验概率——在本例中,相关值超过零的概率。其次,贝叶斯估计将观察到的证据与 _ 先验 _ 相结合,从而使相关估计 _ 正则化,有效地将其拉向零。在这里,我们可以使用 bayesmed 包中的`jzs_cor`函数来计算它。_
```r
bayesCor <- jzs_cor(
hateCrimes$avg_hatecrimes_per_100k_fbi,
hateCrimes$gini_index
)
```
```r
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 50
## Unobserved stochastic nodes: 4
## Total graph size: 230
##
## Initializing model
```
```r
bayesCor
```
```r
## $Correlation
## [1] 0.41
##
## $BayesFactor
## [1] 11
##
## $PosteriorProbability
## [1] 0.92
```
请注意,使用贝叶斯方法估计的相关性略小于使用标准相关系数估计的相关性,这是由于该估计基于证据和先验的组合,从而有效地将估计缩小到反渗透。但是,请注意,贝叶斯分析对异常值不具有鲁棒性,它仍然表示有相当强的证据表明相关性大于零。
## 13.4 相关性和因果关系
当我们说一件事导致另一件事时,我们的意思是什么?关于因果关系意义的讨论在哲学上有着悠久的历史,但在统计学上,我们通常认为因果关系的一种方式是实验控制。也就是说,如果我们认为因子 x 导致因子 y,那么操纵 x 的值也应该操纵 y 的值。
在医学上,有一套被称为[_koch 的假设 _](https://en.wikipedia.org/wiki/Koch%27s_postulates)的观点,在历史上一直被用来确定一个特定的有机体是否引起疾病。基本思想是,有机体应该存在于有疾病的人身上,而不存在于没有疾病的人身上——因此,消除有机体的治疗也应该消除疾病。此外,感染有机体的人应使他们感染该疾病。巴里·马歇尔博士的工作就是一个例子,他假设胃溃疡是由一种细菌(_ 幽门螺杆菌 _)引起的。为了证明这一点,他感染了这种细菌,很快他的胃就出现了严重的炎症。然后他用抗生素治疗自己,他的胃很快就恢复了。他后来因这项工作获得了诺贝尔医学奖。
通常我们想测试因果假设,但实际上我们不能做实验,因为这是不可能的(“人类碳排放与地球气候之间的关系是什么?”)或不道德(“严重虐待对儿童大脑发育有什么影响?”)但是,我们仍然可以收集与这些问题相关的数据。例如,在后一个例子中,我们可以潜在地从遭受虐待的儿童和未遭受虐待的儿童那里收集数据,然后我们可以询问他们的大脑发育是否不同。
假设我们做了这样的分析,我们发现被虐待儿童的大脑发育比未被虐待儿童差。这是否表明虐待会导致大脑发育不良?不,当我们观察到两个变量之间的统计关联时,这两个变量中的一个肯定会引起另一个。然而,这两个变量都有可能受到第三个变量的影响;在这个例子中,虐待儿童可能与家庭压力有关,家庭压力也可能通过较少的智力投入、食物压力或许多其他可能导致大脑发育不良。大道。重点是,两个变量之间的相关性通常告诉我们有什么东西导致了其他事情,但它并不能告诉我们是什么导致了什么。正如统计学家 EdwardTufte 所说,“相关性并不意味着因果关系,但它是一个很好的提示。”
#
## 13.4.1 因果图
描述变量之间因果关系的一种有用方法是通过 _ 因果图 _,它将变量显示为圆,并将变量之间的因果关系显示为箭头。例如,图[13.5](#fig:simpleCausalGraph)显示了学习时间和我们认为应该受到影响的两个变量之间的因果关系:考试成绩和考试完成时间。
图 13.5 显示三个变量之间因果关系的图表:学习时间、考试成绩和考试结束时间。绿色箭头表示一种积极的关系(即学习时间越长,考试成绩越高),红色箭头表示一种消极的关系(即学习时间越长,考试完成越快)。
然而,事实上,对完成时间和成绩的影响并不是直接由花费在学习上的时间量造成的,而是由学生通过学习获得的知识量造成的。我们通常会说知识是一个潜在的(htg0)变量——也就是说,我们不能直接测量它,但是我们可以看到它反映在我们可以测量的变量中(比如成绩和完成时间)。图[13.6](#fig:latentCausalGraph)显示了这一点。
图 13.6 显示了与上述相同的因果关系的图,但现在也显示了使用平方框的潜在变量(知识)。
在这里,我们可以说知识(htg0)介导了学习时间和成绩/完成时间之间的关系。这意味着,如果我们能够保持知识的恒定性(例如,通过服用一种能立即引起遗忘的药物),那么学习时间的长短将不再对成绩和完成时间产生影响。
请注意,如果我们简单地测量考试成绩和完成时间,我们通常会看到他们之间的负面关系,因为完成考试最快的人通常得到最高的分数。然而,如果我们将这种相关性解释为因果关系,这将告诉我们,为了获得更好的成绩,我们实际上应该更快地完成考试!这个例子说明了从非实验数据推断因果关系是多么困难。
在统计学和机器学习领域,有一个非常活跃的研究团体,目前正在研究从非实验数据推断因果关系的时间和方式问题。然而,这些方法往往需要强有力的假设,通常必须谨慎使用。
## 13.5 阅读建议
* [朱迪亚·珀尔的《为什么》(htg1)一书-对因果推理背后的思想的极好介绍。](http://bayes.cs.ucla.edu/WHY/)
\ No newline at end of file
# 13 建模持续关系
大多数人都熟悉 _ 相关 _ 的概念,在本章中,我们将对这个常用和误解的概念提供更正式的理解。
## 13.1 一个例子:仇恨犯罪和收入不平等
2017 年,网站 fivethirtyeight.com 发表了一篇题为“仇恨犯罪率上升与收入不平等(htg1)”的文章,讨论了 2016 年总统选举后仇恨犯罪率与收入不平等之间的关系。报道分析了来自联邦调查局和南方贫困法中心的仇恨犯罪数据,并据此报告:
> “我们发现,收入不平等是全美国调整人口仇恨犯罪和仇恨事件的最重要决定因素”。
此分析的数据包含在`fivethirtyeight`r 包中,这使得我们很容易访问它们。报道中的分析集中在收入不平等(定义为一个叫做 _ 基尼指数 _ 的量)与各州仇恨犯罪流行率之间的关系。
#
## 13.1.1 量化不平等:基尼指数
在我们查看报道中的分析之前,首先要了解如何使用基尼指数来量化不平等。基尼指数通常用一条曲线来定义,这条曲线描述了收入与收入水平等于或小于该水平的人口比例之间的关系,称为 _ 洛伦兹曲线 _。然而,另一种更直观的思考方式是:收入之间的相对平均绝对差异除以二(摘自[https://en.wikipedia.org/wiki/gini_coefficient](https://en.wikipedia.org/wiki/Gini_coefficient)):
![](img/1283ad31298b6a845d0b54ea5116c77d.jpg)
![Lorenz curves for A) perfect equality, B) normally distributed income, and C) high inequality (equal income except for one very wealthy individual).](img/file76.png)
图 13.1 洛伦兹曲线表示 a)完全平等,b)正态分布收入,c)高度不平等(除一个非常富有的个人外,收入相等)。
[13.1](#fig:gini0)显示了几种不同收入分配的洛伦兹曲线。左上方的面板(A)显示了一个例子,其中有 10 个人,每个人的收入完全相同。两个点之间的间隔长度相等,表明每个人在总收入中所占的份额相同。右上角的面板(B)显示了一个收入正态分布的例子。左下角的面板显示了一个不平等程度很高的例子:每个人的收入都是平等的(40000 美元),只有一个人的收入是 40000000 美元。根据美国人口普查,2010 年美国的基尼指数为 0.469,大约在我们的正态分布和最大不相等的例子之间下降了一半。
## 13.2 收入不平等是否与仇恨犯罪有关?
现在我们了解了基尼指数,我们可以看看收入不平等与仇恨犯罪率之间的关系(见图[13.2](#fig:hateCrimeGini))。
```r
hateCrimes <-
hate_crimes %>%
mutate(state_abb = state.abb[match(state,state.name)]) %>%
drop_na(avg_hatecrimes_per_100k_fbi)
hateCrimes$state_abb[hateCrimes$state=="District of Columbia"]='DC'
ggplot(hateCrimes,aes(gini_index,avg_hatecrimes_per_100k_fbi,label=state_abb)) +
geom_point() +
geom_text(aes(label=state_abb),hjust=0, vjust=0) +
theme(plot.title = element_text(size = 20, face = "bold")) +
xlab('Gini index') +
ylab('Avg hate crimes per 100K population (FBI)') +
theme(plot.margin = unit(c(1,1,1,1), "cm"))
```
![Plot of rates of hate crimes vs. Gini index.](img/file77.png)
图 13.2 仇恨犯罪率与基尼指数的关系图。
从数据来看,这两个变量之间似乎有一个正的关系。我们如何量化这种关系?
## 13.3 协方差和相关性
量化两个变量之间关系的一种方法是 _ 协方差 _。记住,单个变量的方差计算如下:
![](img/134e6de0f245b57b2f07e6a35831cb9e.jpg)
这告诉我们每个观察值与平均值相差多远。协方差告诉我们两个不同的变量在观测值之间的偏差是否存在关系。定义如下:
![](img/b268259e45c69b4bfdac3da4ed6962b5.jpg)
当 x 和 y 都高度偏离平均值时,该值将远离零;如果它们在同一方向上偏离,则协方差为正,而如果它们在相反方向上偏离,则协方差为负。让我们先看一个玩具的例子。
```r
# create data for toy example of covariance
df <-
tibble(x = c(3, 5, 8, 10, 12)) %>%
mutate(y = x + round(rnorm(n = 5, mean = 0, sd = 2))) %>%
mutate(
y_dev = y - mean(y),
x_dev = x - mean(x)
) %>%
mutate(crossproduct = y_dev * x_dev)
pander(df)
```
<colgroup><col style="width: 6%"> <col style="width: 6%"> <col style="width: 11%"> <col style="width: 11%"> <col style="width: 19%"></colgroup>
| X | 是 | Y 轴偏差 | X 轴偏差 | 叉乘 |
| --- | --- | --- | --- | --- |
| 三 | 1 个 | -6.6 条 | -4.6 节 | 30.36 天 |
| 5 个 | 3 | -4.6 | -第 2.6 条 | 11.96 年 |
| 8 个 | 8 | 0.4 倍 | 0.4 | 0.16 分 |
| 10 个 | 12 个 | 第 4.4 条 | 第 2.4 条 | 10.56 条 |
| 12 | 14 | 第 6.4 条 | 4.4 | 28.16 条 |
```r
# compute covariance
sprintf("sum of cross products = %.2f", sum(df$crossproduct))
```
```r
## [1] "sum of cross products = 81.20"
```
```r
covXY <- sum(df$crossproduct) / (nrow(df) - 1)
sprintf("covariance: %.2f", covXY)
```
```r
## [1] "covariance: 20.30"
```
我们通常不使用协方差来描述变量之间的关系,因为它随数据的总体方差水平而变化。相反,我们通常使用 _ 相关系数 _(通常在统计学家 Karl Pearson 之后称为 _Pearson 相关 _)。通过用两个变量的标准偏差缩放协方差来计算相关性:
![](img/69240baf835a9a8f0311c744556cb85f.jpg)
```r
# compute the correlation coefficient
corXY <- sum(df$crossproduct) / ((nrow(df) - 1) * sd(df$x) * sd(df$y))
sprintf("correlation coefficient = %.2f", corXY)
```
```r
## [1] "correlation coefficient = 0.99"
```
我们还可以使用 r 中的`cor()`函数轻松计算相关值:
```r
# compute r using built-in function
c <- cor(df$x, df$y)
sprintf("correlation coefficient = %.2f", c)
```
```r
## [1] "correlation coefficient = 0.99"
```
相关系数是有用的,因为它在-1 和 1 之间变化,不管数据的性质如何-事实上,我们在讨论影响大小时已经讨论过相关系数。正如我们在上一章关于影响大小的内容中看到的,1 的相关性表示一个完美的线性关系,-1 的相关性表示一个完美的负关系,0 的相关性表示没有线性关系。
我们可以计算仇恨犯罪数据的相关系数:
```r
corGiniHC <-
cor(
hateCrimes$gini_index,
hateCrimes$avg_hatecrimes_per_100k_fbi
)
sprintf('correlation coefficient = %.2f',corGiniHC)
```
```r
## [1] "correlation coefficient = 0.42"
```
#
## 13.3.1 相关性假设检验
相关值 0.42 似乎表明两个变量之间的关系相当强,但我们也可以想象,即使没有关系,这种情况也可能是偶然发生的。我们可以使用一个简单的公式来测试相关性为零的空假设,该公式允许我们将相关性值转换为 _t_ 统计:
![](img/8a7e337cf588bdff49bc2f192c7eac97.jpg)
在零假设![](img/fed676a693d5c98216a52108307827c6.jpg)下,该统计量以自由度为 t 分布。我们可以使用 r 中的`cor.test()`函数计算:
```r
# perform correlation test on hate crime data
cor.test(
hateCrimes$avg_hatecrimes_per_100k_fbi,
hateCrimes$gini_index
)
```
```r
##
## Pearson's product-moment correlation
##
## data: hateCrimes$avg_hatecrimes_per_100k_fbi and hateCrimes$gini_index
## t = 3, df = 50, p-value = 0.002
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.16 0.63
## sample estimates:
## cor
## 0.42
```
这个测试表明,R 值的可能性很低,这个极限或更高,所以我们将拒绝![](img/6937cddf1a8c4f47df75a6977e1234d1.jpg)的无效假设。注意,这个测试假设两个变量都是正态分布的。
我们也可以通过随机化来检验这一点,在随机化中,我们重复地改变其中一个变量的值并计算相关性,然后将我们观察到的相关性值与这个零分布进行比较,以确定我们观察到的值在零假设下的可能性。结果如图[13.3](#fig:shuffleCorr)所示。使用随机化计算的 p 值与 t 检验给出的答案相当相似。
```r
# compute null distribution by shuffling order of variable values
# create a function to compute the correlation on the shuffled values
shuffleCorr <- function(x, y) {
xShuffled <- sample(x)
return(cor(xShuffled, y))
}
# run this function 2500 times
shuffleDist <-
replicate(
2500,
shuffleCorr(hateCrimes$avg_hatecrimes_per_100k_fbi, hateCrimes$gini_index)
)
```
![Histogram of correlation values under the null hypothesis, obtained by shuffling values. Observed value is denoted by blue line.](img/file78.png)
图 13.3 零假设下相关值的柱状图,通过改变值获得。观测值用蓝线表示。
#
## 13.3.2 稳健相关性
在图[13.2](#fig:hateCrimeGini)中,您可能注意到了一些有点奇怪的地方——其中一个数据点(哥伦比亚特区的数据点)似乎与其他数据点非常不同。我们称之为 _ 离群值 _,标准相关系数对离群值非常敏感。例如,在图[13.4](#fig:outlierCorr)中,我们可以看到一个孤立的数据点是如何导致非常高的正相关值的,即使其他数据点之间的实际关系是完全负的。
![An simulated example of the effects of outliers on correlation. Without the outlier the remainder of the datapoints have a perfect negative correlation, but the single outlier changes the correlation value to highly positive.](img/file79.png)
图 13.4 异常值对相关性影响的模拟示例。如果没有离群值,其余数据点具有完全的负相关,但单个离群值将相关值更改为高度正相关。
解决离群值问题的一种方法是在排序后,在数据的列组上计算相关性,而不是在数据本身上计算相关性;这被称为 _ 斯皮尔曼相关性 _。图[13.4](#fig:outlierCorr)中的 Pearson 相关性为 0.83,而 Spearman 相关性为-0.45,表明等级相关性降低了异常值的影响。
我们可以使用`cor.test`函数计算仇恨犯罪数据的等级相关性:
```r
corTestSpearman <- cor.test( hateCrimes$avg_hatecrimes_per_100k_fbi,
hateCrimes$gini_index,
method = "spearman")
corTestSpearman
```
```r
##
## Spearman's rank correlation rho
##
## data: hateCrimes$avg_hatecrimes_per_100k_fbi and hateCrimes$gini_index
## S = 20000, p-value = 0.8
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.033
```
现在我们看到相关性不再显著(实际上接近于零),这表明 Fivethirtyeight 博客帖子的声明可能由于离群值的影响而不正确。
#
## 13.3.3 贝叶斯相关分析
我们也可以使用贝叶斯分析来分析五个第八个数据,这有两个优点。首先,它为我们提供了一个后验概率——在本例中,相关值超过零的概率。其次,贝叶斯估计将观察到的证据与 _ 先验 _ 相结合,从而使相关估计 _ 正则化,有效地将其拉向零。在这里,我们可以使用 bayesmed 包中的`jzs_cor`函数来计算它。_
```r
bayesCor <- jzs_cor(
hateCrimes$avg_hatecrimes_per_100k_fbi,
hateCrimes$gini_index
)
```
```r
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 50
## Unobserved stochastic nodes: 4
## Total graph size: 230
##
## Initializing model
```
```r
bayesCor
```
```r
## $Correlation
## [1] 0.41
##
## $BayesFactor
## [1] 11
##
## $PosteriorProbability
## [1] 0.92
```
请注意,使用贝叶斯方法估计的相关性略小于使用标准相关系数估计的相关性,这是由于该估计基于证据和先验的组合,从而有效地将估计缩小到反渗透。但是,请注意,贝叶斯分析对异常值不具有鲁棒性,它仍然表示有相当强的证据表明相关性大于零。
## 13.4 相关性和因果关系
当我们说一件事导致另一件事时,我们的意思是什么?关于因果关系意义的讨论在哲学上有着悠久的历史,但在统计学上,我们通常认为因果关系的一种方式是实验控制。也就是说,如果我们认为因子 x 导致因子 y,那么操纵 x 的值也应该操纵 y 的值。
在医学上,有一套被称为[_koch 的假设 _](https://en.wikipedia.org/wiki/Koch%27s_postulates)的观点,在历史上一直被用来确定一个特定的有机体是否引起疾病。基本思想是,有机体应该存在于有疾病的人身上,而不存在于没有疾病的人身上——因此,消除有机体的治疗也应该消除疾病。此外,感染有机体的人应使他们感染该疾病。巴里·马歇尔博士的工作就是一个例子,他假设胃溃疡是由一种细菌(_ 幽门螺杆菌 _)引起的。为了证明这一点,他感染了这种细菌,很快他的胃就出现了严重的炎症。然后他用抗生素治疗自己,他的胃很快就恢复了。他后来因这项工作获得了诺贝尔医学奖。
通常我们想测试因果假设,但实际上我们不能做实验,因为这是不可能的(“人类碳排放与地球气候之间的关系是什么?”)或不道德(“严重虐待对儿童大脑发育有什么影响?”)但是,我们仍然可以收集与这些问题相关的数据。例如,在后一个例子中,我们可以潜在地从遭受虐待的儿童和未遭受虐待的儿童那里收集数据,然后我们可以询问他们的大脑发育是否不同。
假设我们做了这样的分析,我们发现被虐待儿童的大脑发育比未被虐待儿童差。这是否表明虐待会导致大脑发育不良?不,当我们观察到两个变量之间的统计关联时,这两个变量中的一个肯定会引起另一个。然而,这两个变量都有可能受到第三个变量的影响;在这个例子中,虐待儿童可能与家庭压力有关,家庭压力也可能通过较少的智力投入、食物压力或许多其他可能导致大脑发育不良。大道。重点是,两个变量之间的相关性通常告诉我们有什么东西导致了其他事情,但它并不能告诉我们是什么导致了什么。正如统计学家 EdwardTufte 所说,“相关性并不意味着因果关系,但它是一个很好的提示。”
#
## 13.4.1 因果图
描述变量之间因果关系的一种有用方法是通过 _ 因果图 _,它将变量显示为圆,并将变量之间的因果关系显示为箭头。例如,图[13.5](#fig:simpleCausalGraph)显示了学习时间和我们认为应该受到影响的两个变量之间的因果关系:考试成绩和考试完成时间。
图 13.5 显示三个变量之间因果关系的图表:学习时间、考试成绩和考试结束时间。绿色箭头表示一种积极的关系(即学习时间越长,考试成绩越高),红色箭头表示一种消极的关系(即学习时间越长,考试完成越快)。
然而,事实上,对完成时间和成绩的影响并不是直接由花费在学习上的时间量造成的,而是由学生通过学习获得的知识量造成的。我们通常会说知识是一个潜在的(htg0)变量——也就是说,我们不能直接测量它,但是我们可以看到它反映在我们可以测量的变量中(比如成绩和完成时间)。图[13.6](#fig:latentCausalGraph)显示了这一点。
图 13.6 显示了与上述相同的因果关系的图,但现在也显示了使用平方框的潜在变量(知识)。
在这里,我们可以说知识(htg0)介导了学习时间和成绩/完成时间之间的关系。这意味着,如果我们能够保持知识的恒定性(例如,通过服用一种能立即引起遗忘的药物),那么学习时间的长短将不再对成绩和完成时间产生影响。
请注意,如果我们简单地测量考试成绩和完成时间,我们通常会看到他们之间的负面关系,因为完成考试最快的人通常得到最高的分数。然而,如果我们将这种相关性解释为因果关系,这将告诉我们,为了获得更好的成绩,我们实际上应该更快地完成考试!这个例子说明了从非实验数据推断因果关系是多么困难。
在统计学和机器学习领域,有一个非常活跃的研究团体,目前正在研究从非实验数据推断因果关系的时间和方式问题。然而,这些方法往往需要强有力的假设,通常必须谨慎使用。
## 13.5 阅读建议
* [朱迪亚·珀尔的《为什么》(htg1)一书-对因果推理背后的思想的极好介绍。](http://bayes.cs.ucla.edu/WHY/)
\ No newline at end of file
此差异已折叠。
此差异已折叠。
此差异已折叠。
此差异已折叠。
此差异已折叠。
此差异已折叠。
此差异已折叠。
此差异已折叠。
此差异已折叠。
此差异已折叠。
此差异已折叠。
此差异已折叠。
此差异已折叠。
此差异已折叠。
此差异已折叠。
此差异已折叠。
此差异已折叠。
此差异已折叠。
此差异已折叠。
此差异已折叠。
此差异已折叠。
此差异已折叠。
此差异已折叠。
此差异已折叠。
Markdown is supported
0% .
You are about to add 0 people to the discussion. Proceed with caution.
先完成此消息的编辑!
想要评论请 注册