《R语言实战》第7章

# 第七章 基本统计分析
# 本章内容
#   描述性统计分析 
#   频数表和列联表 
#   相关系数和协方差 
#   t检验
#   非参数统计

# 7.1 描述性统计分析
# 本节中，我们将关注分析连续型变量的中心趋势、变化性和分布形状的方法。为了便于说明， 我们将使用第1章中Motor Trend 杂志的车辆路试(mtcars)数据集。我们的关注焦点是每加仑 汽油行驶英里数(mpg)、马力(hp)和车重(wt)。
vars <- c("mpg", "hp", "wt")
head(mtcars[vars])
# 我们将首先查看所有32种车型的描述性统计量，然后按照变速箱类型(am)和汽缸数(cyl)考察描述性统计量。变速箱类型是一个以0表示自动挡、1表示手动挡来编码的二分变量，而汽缸 数可为4、5或6。

# 7.1.1 方法云集
# 对于基础安装，你可以使用summary()函数来获取描述性统计量。
# 代码清单7-1 通过summary()计算描述性统计量
summary(mtcars[vars])
# summary()函数提供了最小值、最大值、四分位数和数值型变量的均值，以及因子向量和逻 辑型向量的频数统计。
# 你可以使用第5章中的apply()函数或sapply()函数计算所选择的任意描 述性统计量。对于sapply()函数，其使用格式为:
# sapply(x, function, options)
# 其中的x是你的数据框(或矩阵)，FUN为一个任意的函数。如果指定了options，它们将被传递 给FUN。你可以在这里插入的典型函数有mean、sd、var、min、max、median、length、range 和quantile。函数fivenum()可返回图基五数总括(Tukey’s five-number summary，即最小值、 下四分位数、中位数、上四分位数和最大值)。
# 代码清单7-2 通过sapply()计算描述性统计量
mystats <- function(x, na.omit = FALSE) {
  if (na.omit)
    x <- x[!is.na(x)]
  m <- mean(x)
  n <- length(x)
  s <- sd(x)
  skew <- sum((x - m) ^ 3 / s ^ 3) / n
  kurt <- sum((x - m) ^ 4 / s ^ 4) / n - 3
  return(c(n = n, mean = m, stdev = s, skew = skew, kurtosis = kurt))
}
# 调用上面的函数
sapply(mtcars[vars], mystats)
# 对于样本中的车型，每加仑汽油行驶英里数的平均值为20.1，标准差为6.0。分布呈现右偏(偏 度+0.61)，并且较正态分布稍平(峰度-0.37)。

# 扩展
# 若干用户贡献包都提供了计算描述性统计量的函数，其中包括Hmisc、pastecs和psych。
# Hmisc包中的describe()函数可返回变量和观测的数量、缺失值和唯一值的数目、平均值、 分位数，以及五个最大的值和五个最小的值。
library(Hmisc)
describe(mtcars[vars])
# pastecs包中有一个名为stat.desc()的函数，它可以计算种类繁多的描述性统计量。使用 格式为:
# stat.desc(x, basic = TRUE, desc = TRUE, norm = FALSE, p = 0.95)
# 其中的x是一个数据框或时间序列。若basic=TRUE(默认值)，则计算其中所有值、空值、缺失 值的数量，以及最小值、最大值、值域，还有总和。若desc=TRUE(同样也是默认值)，则计算 中位数、平均数、平均数的标准误、
#  平均数置信度为95%的置信区间、方差、标准差以及变异系 数。最后，若norm=TRUE(不是默认的)，则返回正态分布统计量，包括偏度和峰度(以及它们 的统计显著程度)和Shapiro–Wilk正态检验结果。这里使用了p值来计算平均数的置信区间(默 认置信度为0.95)。
install.packages("pastecs")
# 代码清单7-4 通过pastecs包中的stat.desc()函数计算描述性统计量
library(pastecs)
stat.desc(mtcars[vars])
# psych包也拥有一个名为describe()的函数，它可以计算非缺失值的数量、 平均数、标准差、中位数、截尾均值、绝对中位差、最小值、最大值、值域、偏度、峰度和平均 值的标准误。
install.packages("psych")
library(psych)
describe(mtcars[vars])
# 在前面的示例中，psych包和Hmisc包均提供了名为describe()的函数。R如何知道该 使用哪个呢?简言之，如代码清单7-5所示，最后载入的程序包优先。在这里，psych在 Hmisc之后被载入，然后显示了一条信息，提示Hmisc包中的describe()函数被psych 包中的同名函数所屏蔽(masked)。键入describe()后，R在搜索这个函数时将首先找 到psych包中的函数并执行它。如果你想改而使用Hmisc包中的版本，可以键入 7 Hmisc::describe(mt)。

# 7.1.2 分组计算描述性统计量
# 在比较多组个体或观测时，关注的焦点经常是各组的描述性统计信息，而不是样本整体的描 述性统计信息。
# 代码清单7-6 使用aggregate()分组获取描述性统计量
aggregate(mtcars[vars], by = list(am = mtcars$am), mean)
aggregate(mtcars[vars], by = list(am = mtcars$am), sd)
# 注意list(am=mtcars$am)的使用。如果使用的是list(mtcars$am)，则am列将被标注为 Group.1而不是am。你使用这个赋值指定了一个更有帮助的列标签。如果有多个分组变量，可以使用by=list(name1=groupvar1, name2=groupvar2, ... , groupvarN)这样的语句。
# 遗憾的是，aggregate()仅允许在每次调用中使用平均数、标准差这样的单返回值函数。
# 它无法一次返回若干个统计量。要完成这项任务，可以使用by()函数。格式为:
# by(data, INDICES, FUN)
# 其中data是一个数据框或矩阵，INDICES是一个因子或因子组成的列表，定义了分组，FUN是任 意函数。
# 代码清单7-7 使用by()分组计算描述性统计量
dstats <- function(x)(c(meann = mean(x), sd = sd(x)))
# 以下语句会报错： (list) object cannot be coerced to type 'double'
by(mtcars[vars], mtcars$am, dstats)
# 使用以下语句解决上面的异常
by(as.numeric(unlist(mtcars["mpg"])), mtcars$am, dstats)
by(as.numeric(unlist(mtcars["hp"])), mtcars$am, dstats)
by(as.numeric(unlist(mtcars["wt"])), mtcars$am, dstats)

# 扩展
# doBy包和psych包也提供了分组计算描述性统计量的函数。同样地，它们未随基本安装发布， 必须在首次使用前进行安装。doBy包中summaryBy()函数的使用格式为:
# summaryBy(formula, data = dataframe, FUN = function)
# 其中的formula接受以下的格式:
# var1 + var2 + var3 + ... + varN ~ groupvar1 + groupvar2 +... + groupvarN
# 在~左侧的变量是需要分析的数值型变量，而右侧的变量是类别型的分组变量。function 可为任何内建或用户自编的R函数。
install.packages("doBy")
# 代码清单7-8 使用doBy包中的summaryBy()分组计算概述统计量
library(doBy)
summaryBy(mpg + hp + wt ~ am, data = mtcars, FUN = mystats)
# 与前面的示例不同，describe.by()函数不允许指定任意函数，所以它的普适性较低。若 存在一个以上的分组变量，你可以使用list(groupvar1, groupvar2, ... , groupvarN) 来表示它们。但这仅在分组变量交叉后不出现空白单元时有效。
# 最后，可以使用5.6.3节中描述的reshape包灵活地按组导出描述性统计量。首先，使用:
# dfm <- melt(dataframe, measure.vars = y, id.vars = g)
# 融合数据框。其中的dataframe包含着数据，y是一个向量，指明了要进行概述的数值型变量(默 认使用所有变量)，而g是由一个或多个分组变量组成的向量。然后使用:
# cast(dfm, groupvar1, groupvar2 + ... + variable - ., FUN)
# 重铸数据。分组变量以+号分隔，这里的variable只取其字面含义(即仅表示重铸后数据框中的变量variable)，而FUN是一个任意函数。
library(reshape)
dstats <- function(x) (c(n = length(x), mean = mean(x), sd = sd(x)))
dfm <- melt(mtcars, measure.vars = c("mpg", "hp", "wt"))
cast(dfm, am + cyl + variable ~ ., dstats)


# 7.2 频数表和列联表
# 在本节中，我们将着眼于类别型变量的频数表和列联表，以及相应的独立性检验、相关性的 度量、图形化展示结果的方法。我们除了使用基础安装中的函数，还将连带使用vcd包和gmodels 包中的函数。
# 本节中的数据来自vcd包中的Arthritis数据集。这份数据来自Kock & Edward(1988)，表 示了一项风湿性关节炎新疗法的双盲临床实验的结果。前几个观测是这样的:
library(vcd)
head(Arthritis)
# 治疗情况(安慰剂治疗、用药治疗)、性别(男性、女性)和改善情况(无改善、一定程度的改善、显著改善)均为类别型因子1。下一节中，我们将使用此数据创建频数表和列联表(交叉的分类)。

# 7.2.1 生成频数表
#       表7-1 用于创建和处理列联表的函数
#       函数                              描述
# table(var1, var2, ..., varN)      使用 N 个类别型变量(因子)创建一个 N 维列联表
# xtabs(formula, data)              根据一个公式和一个矩阵或数据框创建一个 N 维列联表
# prop.table(table, margins)        依margins定义的边际列表将表中条目表示为分数形式
# margin.table(table, margins)      依margins定义的边际列表计算表中条目的和
# addmargins(table, margins)        将概述边margins(默认是求和结果)放入表中
# ftable(table)                     创建一个紧凑的“平铺”式列联表

# 1. 一维列联表
# 可以使用table()函数生成简单的频数统计表。示例如下:
mytable <- with(Arthritis, table(Improved))
mytable
# 可以用prop.table()将这些频数转化为比例值:
prop.table(mytable)
# 或使用prop.table()*100转化为百分比:
prop.table(mytable) * 100

# 2. 二维列联表
# 对于二维列联表，table()函数的使用格式为:
# mytable <- table(A, B)
# 其中的A是行变量，B是列变量。除此之外，xtabs()函数还可使用公式风格的输入创建列联表， 格式为:
# mytable <- xtabs(~ A + B, data = mydata)
# 其中的mydata是一个矩阵或数据框。总的来说，要进行交叉分类的变量应出现在公式的右侧(即~ 符号的右方)，以+作为分隔符。若某个变量写在公式的左侧，则其为一个频数向量(在数据已经 被表格化时很有用)。
# 对于Arthritis数据，有:
mytable <- xtabs(~ Treatment + Improved , data = Arthritis)
mytable
# 你可以使用margin.table()和prop.table()函数分别生成边际频数和比例。行和与行比 例可以这样计算:
# 行和
margin.table(mytable, 1)
# 列和
margin.table(mytable, 2)
# 行比
prop.table(mytable, 1)
# 列比
prop.table(mytable, 2)
# 以上四个语句中的下标1指代table()语句中的第一个变量，下标2指代table()语句中的第二个变量。观察表格可以发现，与接受安慰剂的个体中有显著改善的16%相比，接受治疗的个体中的51%的个体病情有了显著的改善。
# 各单元格所占比例可用如下语句获取:
prop.table(mytable)
# 你可以使用addmargins()函数为这些表格添加边际和。例如，以下代码添加了各行的和与各列的和:
addmargins(mytable)
addmargins(prop.table(mytable))
# 在使用addmargins()时，默认行为是为表中所有的变量创建边际和。作为对照:
addmargins(prop.table(mytable, 1), 2)
# 仅添加了各行的和。类似地，
addmargins(prop.table(mytable, 2), 1)
# 添加了各列的和。在表中可以看到，有显著改善患者中的25%是接受安慰剂治疗的。
# table()函数默认忽略缺失值(NA)。要在频数统计中将NA视为一个有效的类别，请设定参数useNA="ifany"
# 使用gmodels包中的CrossTable()函数是创建二维列联表的第三种方法。CrossTable() 函数仿照SAS中PROC FREQ或SPSS中CROSSTABS的形式生成二维列联表。示例参阅代码清单 7-11。
# 代码清单7-11 使用CrossTable生成二维列联表
install.packages("gmodels")
library(gmodels)
CrossTable(Arthritis$Treatment, Arthritis$Improved)
# CrossTable()函数有很多选项，可以做许多事情:计算(行、列、单元格)的百分比;指 定小数位数;进行卡方、Fisher和McNemar独立性检验;计算期望和(皮尔逊、标准化、调整的 标准化)残差;将缺失值作为一种有效值;进行行和列标题的标注;生成SAS或SPSS风格的输出。 参阅help(CrossTable)以了解详情。

# 如果有两个以上的类别型变量，那么你就是在处理多维列联表。我们将在下面考虑这种情况。
# 3. 多维列联表
# table()和xtabs()都可以基于三个或更多的类别型变量生成多维列联表。margin.table()、
#  prop.table()和addmargins()函数可以自然地推广到高于二维的情况。另外，ftable()函 数可以以一种紧凑而吸引人的方式输出多维列联表。
# 代码清单7-12 三维列联表
# 各单元格的频数
mytable <- xtabs(~ Treatment + Sex + Improved, data = Arthritis)
mytable
ftable(mytable)
# 边际频数
margin.table(mytable, 1)
margin.table(mytable, 2)
margin.table(mytable, 3)
# 治疗情况(Treatment) x 改善情况(Improved)的边际频数
margin.table(mytable, c(1, 3))
# 治疗情况(Treatment) x 性别(Sex)的边际频数
ftable(prop.table(mytable, c(1, 2)))
# 为第三个下标添加了边际和
ftable(addmargins(prop.table(mytable, c(1, 2)), 3))
# 如果想得到百分比而不是比例，可以将结果表格乘以100。例如:
ftable(addmargins(prop.table(mytable, c(1, 2)), 3)) * 100

# 7.2.2 独立性检验
# 1. 卡方独立性检验
# 你可以使用chisq.test()函数对二维表的行变量和列变量进行卡方独立性检验。
# 代码清单7-13 卡方独立性检验
library(vcd)
# (1)治疗情况和改善情况不独立
mytable <- xtabs(~ Treatment + Improved, data = Arthritis)
chisq.test(mytable)
# (2)性别和改善情况独立
mytable <- xtabs(~ Treatment + Sex, data = Arthritis)
chisq.test(mytable)
# 在结果(1)中，患者接受的治疗和改善的水平看上去存在着某种关系(p < 0.01)。而患者性别和改善情况之间却不存在关系(p > 0.05)(2)
# 这里的p值表示从总体中抽取的样本行变量与列变量是相互独立的概率。由于(1)的概率值很小，所以你拒绝了治疗类型和治疗结果相互独立的原假 设。由于(2)的概率不够小，故没有足够的理由说明治疗结果和性别之间是不独立的。代码清单7-13 中产生警告信息的原因是，表中的6个单元格之一(男性 - 一定程度上的改善)有一个小于5的值， 这可能会使卡方近似无效。

# 2. Fisher精确检验
# 可以使用fisher.test()函数进行Fisher精确检验。Fisher精确检验的原假设是:边界固定 的列联表中行和列是相互独立的。其调用格式为fisher.test(mytable)，其中的mytable是 一个二维列联表。示例如下:
mytable <- xtabs(~ Treatment + Improved, data = Arthritis)
fisher.test(mytable)
# 与许多统计软件不同的是，这里的fisher.test()函数可以在任意行列数大于等于2的二维 列联表上使用，但不能用于2×2的列联表。

# 3. Cochran-Mantel—Haenszel检验
# mantelhaen.test()函数可用来进行Cochran—Mantel—Haenszel卡方检验，其原假设是，两 个名义变量在第三个变量的每一层中都是条件独立的。下列代码可以检验治疗情况和改善情况在 性别的每一水平下是否独立。此检验假设不存在三阶交互作用(治疗情况×改善情况×性别)。
mytable <- xtabs(~ Treatment + Improved + Sex, data = Arthritis)
mantelhaen.test(mytable)
# 结果表明，患者接受的治疗与得到的改善在性别的每一水平下并不独立(即，分性别来看， 用药治疗的患者较接受安慰剂的患者有了更多的改善)。

# 7.2.3 相关性的度量
# 上一节中的显著性检验评估了是否存在充分的证据以拒绝变量间相互独立的原假设。如果可以拒绝原假设，那么你的兴趣就会自然而然地转向用以衡量相关性强弱的相关性度量。vcd包中的assocstats()函数可以用来计算二维列联表的phi系数、列联系数和Cramer’s V系数。
# 代码清单7-14 二维列联表的相关性度量
library(vcd)
mytable <- xtabs(~ Treatment + Improved, data = Arthritis)
assocstats(mytable)
# 总体来说，较大的值意味着较强的相关性。vcd包也提供了一个kappa()函数，可以计算混淆矩阵的Cohen’s kappa值以及加权的kappa值。(举例来说，混淆矩阵可以表示两位评判者对于一 系列对象进行分类所得结果的一致程度。)


# 7.2.4 结果的可视化
# R中拥有远远超出其他多数统计软件的、可视地探索类别型变量间关系的方法。通常，我们 会使用条形图进行一维频数的可视化(参见6.1节)。vcd包中拥有优秀的、用于可视化多维数据 集中类别型变量间关系的函数，
#  可以绘制马赛克图和关联图(参见11.4节)。最后，ca包中的对 应分析函数允许使用多种几何表示(Nenadic & Greenacre，2007)可视地探索列联表中行和列之 间的关系。

# 7.2.5 将表转换为扁平格式
# 代码清单7-15 通过table2flat将表转换为扁平格式
table2flat <- function(mytable) {
  df <- as.data.frame(mytable)
  rows <- dim(df)[1]
  cols <- dim(df)[2]
  x <- NULL
  for (i in 1:rows) {
    for (j in 1:df$Freq[i]) {
      row <- df[i, c(1:(cols - 1))]
      x <- rbind(x, row)
    }
  }
  row.names(x) <- c(1:dim(x)[1])
  return(x)
}
# 调用
table2flat(mytable)
# 代码清单7-16 使用table2flat()函数转换已发表的数据
treatment <- rep(c("Placebo", "Treated"), times = 3)
improved <- rep(c("None", "Some", "Marked"), each = 2)
Freq <- c(29, 13, 7, 17, 7, 21)
mytable <- as.data.frame(cbind(treatment, improved, Freq)) 
mydata <- table2flat(mytable)
head(mydata)



# 7.3 相关
# 相关系数可以用来描述定量变量之间的关系。相关系数的符号()表明关系的方向(正相 关或负相关)，其值的大小表示关系的强弱程度(完全不相关时为0，完全相关时为1)。
# 本节中，我们将关注多种相关系数和相关性的显著性检验。我们将使用R基础安装中的 state.x77数据集，它提供了美国50个州在1977年的人口、收入、文盲率、预期寿命、谋杀率和 高中毕业率数据。
#  数据集中还收录了气温和土地面积数据，但为了节约空间，这里将其丢弃。你 可以使用help(state.x77)了解数据集的更多信息。除了基础安装以外，我们还将使用psych 和ggm包。
# 7.3.1 相关的类型
# 1. Pearson、Spearman和Kendall相关
# Pearson积差相关系数衡量了两个定量变量之间的线性相关程度。Spearman等级相关系数则衡 量分级定序变量之间的相关程度。Kendall’s Tau相关系数也是一种非参数的等级相关度量。
# cor()函数可以计算这三种相关系数，而cov()函数可用来计算协方差。两个函数的参数有 很多，其中与相关系数的计算有关的参数可以简化为:
# cor(x, use= , method = )
#         表7-3 cor和cov的参数
# 参数          描述
# x         矩阵或数据框
# use       指定缺失数据的处理方式。可选的方式为all.obs(假设不存在缺失数据——遇到缺失数据时将报
#           错)、everything(遇到缺失数据时，相关系数的计算结果将被设为missing)、complete.obs 
#           (行删除)以及 pairwise.complete.obs(成对删除，pairwise deletion)
# method    指定相关系数的类型。可选类型为pearson、spearman或kendall
# 默认参数为use="everything"和method="pearson"。
# 代码清单7-17 协方差和相关系数
states <- state.x77[, 1:6]
cov(states)
cor(states)
cor(states, method = "spearman")
# 以上首个语句计算了方差和协方差，第二个语句则计算了Pearson积差相关系数，而第三个语句计算 了Spearman等级相关系数。举例来说，我们可以看到收入和高中毕业率之间存在很强的正相关， 而文盲率和预期寿命之间存在很强的负相关。
# 请注意，在默认情况下得到的结果是一个方阵(所有变量之间两两计算相关)。你同样可以 计算非方形的相关矩阵。观察以下示例:
x <- states[, c("Population", "Income", "Illiteracy", "HS Grad")]
y <- states[, c("Life Exp", "Murder")]
cor(x, y)
# 当你对某一组变量与另外一组变量之间的关系感兴趣时，cor()函数的这种用法是非常实用 的。注意，上述结果并未指明相关系数是否显著不为0(即，根据样本数据是否有足够的证据得 出总体相关系数不为0的结论)。由于这个原因，你需要对相关系数进行显著性检验(在7.3.2节中 阐述)。

# 2. 偏相关
# 偏相关是指在控制一个或多个定量变量时，另外两个定量变量之间的相互关系。你可以使用 ggm包中的pcor()函数计算偏相关系数。ggm包没有被默认安装，在第一次使用之前需要先进行 安装。函数调用格式为:
# pcor(u, s)
# 其中的u是一个数值向量，前两个数值表示要计算相关系数的变量下标，其余的数值为条件变量 (即要排除影响的变量)的下标。S为变量的协方差阵。这个示例有助于阐明用法:
install.packages("ggm")
library(ggm)
# 在控制了收入、文盲率和高中毕业率时人口和谋杀率的偏相关系数
pcor(c(1, 5, 2, 3, 6), cov(states))
# 本例中，在控制了收入、文盲率和高中毕业率的影响时，人口和谋杀率之间的相关系数为 0.346。偏相关系数常用于社会科学的研究中。

# 3. 其他类型的相关
# polycor包中的hetcor()函数可以计算一种混合的相关矩阵，其中包括数值型变量的 Pearson积差相关系数、数值型变量和有序变量之间的多系列相关系数、
# 有序变量之间的多分格相 关系数以及二分变量之间的四分相关系数。多系列、多分格和四分相关系数都假设有序变量或二 分变量由潜在的正态分布导出。请参考此程序包所附文档以了解更多。

# 7.3.2 相关性的显著性检验
# 在计算好相关系数以后，如何对它们进行统计显著性检验呢?常用的原假设为变量间不相关 (即总体的相关系数为0)。你可以使用cor.test()函数对单个的Pearson、Spearman和Kendall相关系数进行检验。简化后的使用格式为:
# cor.test(x, y, alternative = , method = )
# 其中的x和y为要检验相关性的变量，alternative则用来指定进行双侧检验或单侧检验(取值 为"two.side"、"less"或"greater")，而method用以指定要计算的相关类型("pearson"、 "kendall"或"spearman")。当研究的假设为总体的相关系数小于0时，请使用alternative= "less"。
# 在研究的假设为总体的相关系数大于0时，应使用alternative="greater"。在默认情况下，假设为alternative="two.side"(总体相关系数不等于0)。参考代码清单7-18中的示例。
# 代码清单7-18 检验某种相关系数的显著性
cor.test(states[, 3], states[, 5])
# 这段代码检验了预期寿命和谋杀率的Pearson相关系数为0的原假设。假设总体的相关度为0， 则预计在一千万次中只会有少于一次的机会见到0.703这样大的样本相关度(即p = 1.258e08)。 6 由于这种情况几乎不可能发生，所以你可以拒绝原假设，从而支持了要研究的猜想，即预期寿命 和谋杀率之间的总体相关度不为0。
# 遗憾的是，cor.test每次只能检验一种相关关系。但幸运的是，psych包中提供的 corr.test()函数可以一次做更多事情。corr.test()函数可以为Pearson、Spearman或Kendall 相关计算相关矩阵和显著性水平。
# 代码清单7-19 通过corr.test计算相关矩阵并进行显著性检验
library(psych)
corr.test(states, use = "complete")
# 参数use=的取值可为"pairwise"或"complete"(分别表示对缺失值执行成对删除或行删 除)。参数method=的取值可为"pearson"(默认值)、"spearman"或"kendall"。这里可以看到，人口数量和高中毕业率的相关系数(0.10)并不显著地不为0(p = 0.5)。

# 其他显著性检验
# 在7.4.1节中，我们关注了偏相关系数。在多元正态性的假设下，psych包中的pcor.test() 函数(这里可能是作者的笔误，函数pcor.test事实上包含于ggm包中)可以用来检验在控制一个或多个额外变量时两个变量之间的条件独立性。使用格式为:
# poor.test(r, q, n)
# 其中的r是由pcor()函数计算得到的偏相关系数，q为要控制的变量数(以数值表示位置)，n为样本大小。
# 在结束这个话题之前应当指出的是，psych包中的r.test()函数提供了多种实用的显著性 检验方法。此函数可用来检验:
# 某种相关系数的显著性;
# 两个独立相关系数的差异是否显著;
# 两个基于一个共享变量得到的非独立相关系数的差异是否显著;
# 两个基于完全不同的变量得到的非独立相关系数的差异是否显著。
# 参阅help(r.test)以了解详情。

# 7.3.3 相关关系的可视化
# 以相关系数表示的二元关系可以通过散点图和散点图矩阵进行可视化，而相关图 (correlogram)则为以一种有意义的方式比较大量的相关系数提供了一种独特而强大的方法。



# 7.4 t检验
# 7.4.1 独立样本的t检验
# 如果你在美国的南方犯罪，是否更有可能被判监禁?我们比较的对象是南方和非南方各州，因变量为监禁的概率。一个针对两组的独立样本t检验可以用于检验两个总体的均值相等的假设。这里假设两组数据是独立的，并且是从正态总体中抽得。检验的调用格式为:
# t.test(y ~ x, data)
# 其中的y是一个数值型变量，x是一个二分变量。调用格式或为:
# t.test(y1, y2)
# 其中的y1和y2为数值型向量(即各组的结果变量)。可选参数data的取值为一个包含了这些 变量的矩阵或数据框。与其他多数统计软件不同的是，这里的t检验默认假定方差不相等，并使 用Welsh的修正*度。
#   你可以添加一个参数var.equal=TRUE以假定方差相等，并使用合并方 差估计。默认的备择假设是双侧的(即均值不相等，但大小的方向不确定)。你可以添加一个参 数alternative="less"或alternative="greater"来进行有方向的检验。
# 在下列代码中，我们使用了一个假设方差不等的双侧检验，比较了南方(group 1)和非南 方(group 0)各州的监禁概率:
library(MASS)
t.test(Prob ~ So, data = UScrime)
# 你可以拒绝南方各州和非南方各州拥有相同监禁概率的假设(p < .001)。
# 由于结果变量是一个比例值，你可以在执行t检验之前尝试对其进行正态化变换。在本例 中，所有对结果变量合适的变换(Y/1-Y、log(Y/1-Y)、arcsin(Y)、arcsin(sqrt(Y)) 都会将检验引向相同的结论。数据变换详述于第8章。

# 7.4.2 非独立样本的t检验
# 非独立样本的t检验假定组间的差异呈正态分布。对于本例，检验的调用格式为:
# t.test(y1, y2, paired = TRUE)
# 其中的y1和y2为两个非独立组的数值向量。结果如下:
library(MASS)
sapply(UScrime[c("U1", "U2")], function(x)(c(mean = mean(x), sd = sd(x))))
with(UScrime, t.test(U1, U2, paired = TRUE))
# 差异的均值(61.5)足够大，可以保证拒绝年长和年轻男性的平均失业率相同的假设。 年轻男性的失业率更高。事实上，若总体均值相等，获取一个差异如此大的样本的概率小于 0.000 000 000 000 000 22(即2.2e16)。

# 7.4.3 多于两组的情况
# 如果想在多于两个的组之间进行比较，应该怎么做?如果能够假设数据是从正态总体中独立 抽样而得的，那么你可以使用方差分析(ANOVA)。ANOVA是一套覆盖了许多实验设计和准实验 设计的综合方法。就这一点来说，它的内容值得单列一章。你可以随时离开本节转而阅读第9章。



# 7.5 组间差异的非参数检验
# 如果数据无法满足t检验或ANOVA的参数假设，可以转而使用非参数方法。举例来说，若结果变量在本质上就严重偏倚或呈现有序关系，那么你可能会希望使用本节中的方法。
# 7.5.1 两组的比较
# 若两组数据独立，可以使用Wilcoxon秩和检验(更广为人知的名字是Mann–Whitney U检验) 来评估观测是否是从相同的概率分布中抽得的(即，在一个总体中获得更高得分的概率是否比另 一个总体要大)。调用格式为:
# wilcox.test(y ~ x, data)
# 其中的y是数值型变量，而x是一个二分变量。调用格式或为:
# wilcox.test(y1, y2)
# 其中的y1和y2为各组的结果变量。可选参数data的取值为一个包含了这些变量的矩阵或数据框。默 认进行一个双侧检验。你可以添加参数exact来进行精确检验，指定alternative="less"或 alternative="greater"进行有方向的检验。
# 如果你使用Mann–Whitney U检验回答上一节中关于监禁率的问题，将得到这些结果:
with(UScrime, by(Prob, So, median))
wilcox.test(Prob ~ So, data = UScrime)
# 你可以再次拒绝南方各州和非南方各州监禁率相同的假设(p < 0.001)。
# Wilcoxon符号秩检验是非独立样本t检验的一种非参数替代方法。它适用于两组成对数据和 无法保证正态性假设的情境。调用格式与Mann–Whitney U检验完全相同，不过还可以添加参数 paired=TRUE。让我们用它解答上一节中的失业率问题:
sapply(UScrime[c("U1", "U2")], median)
with(UScrime, wilcox.test(U1, U2, paired = TRUE))
# 你再次得到了与配对t检验相同的结论。
# 在本例中，含参的t检验和与其作用相同的非参数检验得到了相同的结论。当t检验的假设合 理时，参数检验的功效更强(更容易发现存在的差异)。而非参数检验在假设非常不合理时(如对于等级有序数据)更适用。


# 7.5.2 多于两组的比较
# 如果无法满足ANOVA设计的假设，那么可以使用非参数方法来评估组间的差异。如果各组独立，则Kruskal—Wallis检验将是一种实用的方法。如果各组不独立(如重复测量设计或随机区 组设计)，那么Friedman检验会更合适。Kruskal–Wallis检验的调用格式为:
kruskal.test(y ~ A, data)
# 其中的y是一个数值型结果变量，A是一个拥有两个或更多水平的分组变量(grouping variable)。 (若有两个水平，则它与Mann–Whitney U检验等价。)而Friedman检验的调用格式为:
# friedman.test(y ~ A | B, data)
# 其中的y是数值型结果变量，A是一个分组变量，而B是一个用以认定匹配观测的区组变量(blocking variable)。在以上两例中，data皆为可选参数，它指定了包含这些变量的矩阵或数据框。
# 让我们利用Kruskal–Wallis检验回答文盲率的问题。首先，你必须将地区的名称添加到数据 集中。这些信息包含在随R基础安装分发的state.region数据集中。
states <- as.data.frame(cbind(state.region, state.x77))
# 现在就可以进行检验了:
kruskal.test(Illiteracy ~ state.region, data = states)
# 显著性检验的结果意味着美国四个地区的文盲率各不相同(p <0.001)。
# 虽然你可以拒绝不存在差异的原假设，但这个检验并没有告诉你哪些地区显著地与其他地区 不同。要回答这个问题，你可以使用Mann–Whitney U检验每次比较两组数据。一种更为优雅的 方法是在控制犯第一类错误的概率(发现一个事实上并不存在的差异的概率)的前提下，执行可 以同步进行的多组比较，这样可以直接完成所有组之间的成对比较。npmc包提供了所需要的非 参数多组比较程序。
# 代码清单7-20 非参数多组比较
class <- state.region
var <- state.x77[, c("Illiteracy")]
mydata <- as.data.frame(cbind(class, var))
rm(class, var)
# 安装npmc包报错： package ‘npmc’ is not available (for R version 3.5.2)
install.packages("npmc")
library(npmc)
summary(npmc(mydata), type = "BF")
# 调用了npmc的语句生成了六对统计比较结果(东北部对南部、东北部对中北部、东北部对。可以从双侧的p值(p.value.2s)看 处可以看到南部的文西部、南部对中北部、南部对西部，以及中北部对西部) 出南部与其他三个地区显著不同，而其他三个地区之间并没有什么不同。在 盲率中间值更高。注意，npmc在计算积分时使用了随机数，所以每次计算的结果会有轻微的不同



# 7.6 组间差异的可视化
# 在7.4节和7.5节中，我们关注了进行组间比较的统计方法。使用视觉直观地检查组间差异，同 样是全面的数据分析策略中的一个重要组成部分。它允许你评估差异的量级、甄别出任何会影响结 果的分布特征(如偏倚、双峰或离群点)并衡量检验假定的合理程度。R中提供了许多比较组间数 11 据的图形方法，其中包括6.5节中讲解的箱线图(简单箱线图、含凹槽的箱线图、小提琴图)、6.4.1 节中叠加的核密度图，以及在第9章中讨论的评估检验假定的图形方法。