####第七章####
####7.1 描述性统计分析####
#本节将关注连续型变量的中心趋势、变化性和分布形状的分析方法。
#本节使用第1章中MotorTrend杂志的mtcars数据集中的一此变量。
#我们的关注重点是每加仑汽油行驶英里数(mpg)、马力(hp）和车重（wt）。
myvars <- c("mpg", "hp", "wt")
head(mtcars[myvars])
#首先查看所有32种车型的描述性统计量，然后按变速箱类型(am)和发动机汽缸配置(vs)分组考察描述性统计量。
#变速箱类型是一个以0表示自动挡(0=automatic)、1表示手动挡（1=manual）来编码的变量。
#汽缸配置的编码为0=V-shape、1=straight。
#(`head()`函数是一个非常有用的函数，用于查看数据对象（通常是向量、矩阵、数据框等）的前几行或前几个元素。
#默认情况下，`head()`函数返回对象的前6行或元素。
#这对于快速查看大型数据集或矩阵的开头部分非常有用，而无需打印整个对象。)
##7.1.1方法云集
#代码清单7-1用summary()函数来获取描述性统计量
myvars <- c("mpg", "hp", "wt")
summary(mtcars[myvars])
#summary()函数描述了最小值、最大值、四分位数和均值。
#第五章中的apply()和sapply()函数计算任意描述性统计量
#apply()用于矩阵，sapply()用于数据框
#sapply()的格式如下：
sapply(x, FUN, options)
#x为数据框，FUN为任意一个函数，如果指定了options，则这个选项直接传递给FUN
#这里可以插入的函数有mean(),sd(),var(),max(),median(),length(),range()和quantile()
#函数fivenum()可以返回图基五数总括（最小值，下四分位数，中位数，上四分位数和最大值）
#代码清单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))
}
myvars <- c("mpg", "hp", "wt")
sapply(mtcars[myvars], mystats)
##7.1.2 更多方法
#若干用户贡献的包都提供了计算描述性统计量的函数，
#其中包括Hmisc、pastecs、psych、skimr和summytools。
#这里只显示前3个。
install.packages("Hmisc")
install.packages("pastecs")
install.packages("psych")
#代码7-3用Hmisc包中的describe()包返回变量和观测值的数量、缺失值、
#和唯一值的数目、均值、分位数、5个最大值和最小值，代码如下：
library(Hmisc)
myvars <- c("mpg", "hp", "wt")
describe(mtcars[myvars])
#pastecs包中有一个为名stat.desc()的函数，可以计算种类繁多的描述性统计量，格式如下：
stat.desc(x, basic = TRUE, desc = TRUE, norm = TRUE, p = 0.95)
#其中的x是一个数据框或时间序列。
#若basic=TRUE（默认值)，则计算其中所有值、空值、缺失值的数量，以及最小值、最大值、值域，还有总和。
#若desc=TRUE（同样也是默认值）,则计算中位数、均值、均值的标准误、平均数置信度为95%的置信区间、方差、标准差以及变异系数最后.
#若norm=TRUE（不是默认的)，则返回正态分布统计量，包括偏度和峰度（以及它们的统计显著程度）和Shapiro-Wilk正态检验结果。
#这里使用了p值来计算平均数的置信区间（默认置信度为0.95）。
#代码清单7-4提供了一个示例:
library(pastecs)
myvars <- c("mpg", "hp", "wt")
stat.desc(mtcars[myvars])
#psych包也有一个名为describe(）的函数，
#它可以计算观测值中非缺失值的数量、平均数、标准差、中位数、截尾均值、绝对中位差、最小值、最大值、值域、偏度、峰度和均值的标准误。
#代码清单7-5中提供了一个示例:
library(psych)
myvars <- c("mpg", "hp", "wt")
describe(mtcars[myvars])
#注意：在前面的示例中，psych包和Hmisc包均提供了名为describe()的函数。
#R如何知道该使用哪个呢？简言之，如代码清单7-5所示，最后载入的包优先。
#在这里，psych在Hmisc之后被载入，然后显示了一条信息。
#提示Hmisc包中的describe()函数被psych句7中的同名函数所屏蔽（masked）。
#输入describe()后，R在搜索这个函数时将首先找到psych包中的函数并执行它。
#如果我们还是想使用Hmisc包中的版本，可以输入Hmisc:describe(mt)。
#这个函数仍然在那里，我们只是需要给予R更多信息来找到它。
##7.1.3分组计算描述性统计量
#可以使用by()函数计算分组统计量，格式如下：
by(data, INDICES, FUN)
#其中的data是一个数据框或矩阵，INDICES是一个因子或因子组成的列表，定义了分组.
#FUN可以是任意一个函数，用来操作数据框中的所有列。
#代码清单7-6提供了一个示例。
dstats <- function(x)sapply(x, mystats) #这里定义了一个新的函数dstats，调用了代码7-2种的mystats函数
myvars <- c("mpg", "hp", "wt")
by(mtcars[myvars], mtcars$am, dstats) #按am进行分组
#接下来的示例（代码清单7-7）根据两个分组变量（am和vs）生成了汇总统计量，并使用自定义标签输出了各组的结果。
#另外，在计算统计量的时候忽略了缺失值。
dstats <- function(x)sapply(x, mystats, na.omit = TRUE)
myvars <- c("mpg", "hp", "wt")
by(mtcars[myvars],
list(Transmission = mtcars$am,
Engine = mtcars$vs),
FUN = dstats)
#虽然前面的示例使用的是mystats()函数，但我们也可以使用Hmisc和psych包中的describe()函数或者pastecs包中的stat.desc()函数。
#事实上，by()函数提供了通用的处理机制，可以针对分组逐一进行任意的分析操作。
##7.1.4使用dplyr济学宁交互式汇总数据
#到目前为止，我们所讨论的函数集中在为整个数据框计算综合的描述性统计量。
#但是，在交互式、探索性数据分析中，我们的目标是回答具有针对性的问题。
#在这种情况下，我们希望获取特定观测值的特定的统计量。
#3.11节介绍的dolyr包为我们提供了快速并灵活地达成此目标的工具。
#函数summarize()和summarize_all()可以用来计算任何指定的统计量，函数group_by()可以用来指定要计算这些统计量的分组。
#在此举个示例，我们使用carData包中的Salaries数据框来询问并回答一系列问题。
#这个数据集包含美国一所大学在2008年至2009年之间397名教授9个月的薪水（以美元为单位）数据。
#这些数据是一项正在进行的，旨在监测男性教师和女性教师的薪水差异项目的一部分。
#在继续操作之前，请确保已安装carData包和dplyr包:
install.packages(c("carData", "dplyr"))
library(carData)
library(dplyr)
#现在397名教授的薪水中位数和薪水范围是多少？
Salaries %>% summarize(med = median(salary),
min = min(salary),
max = max(salary))
#不同性别和级别的教授数量，薪水的中位数和薪水范围是多少？
Salaries %>%
group_by(rank, sex) %>%
summarize(n = length(salary),
med = median(salary),
min = min(salary),
max = max(salary))
#不同性别和级别的教授的平均任职年限和获得博士学位后的平均任职年限是多少？
Salaries %>%
group_by(rank, sex) %>%
select(yrs.service, yrs.since.phd) %>%
summarise_all(mean)
#函数summarize_all（)为每个非分组变量（yrs.service和yrs.since.phd）的计算汇总统计量。
#如果我们希望为每个变量计算多个统计量，则可以提供一个清单。
#例如，summarize_all(list(mean=mean，std=sd))将计算每个变量的均值和标准差。
#使用dplyr包的一个优点是结果以tibble（数据框）的形式返回，便于我们进一步分析这些汇总结果、绘制图形、重新设置汇总结果的格式并打印出来。
#dplyr包还提供了汇总数据的简易机制。
#一般来说，数据分析人员对于输出哪些描述性统计量以及如何设置它们的格式有着自己的偏好。
##7.1.5 结果的可视化
#分布特征的数值刻画的确很重要，但是这并不能代替数据的视觉呈现。
#对于定量变量，我们有直方图（6.4节）、核密度图（6.5节）、箱线图（6.6节）和点图（6.7节），它们都可以让我们洞悉那些依赖于观察一小部分描述性统计量时忽略的细节。
####7.2 频数表和列联表####
#在本节中，我们将着眼于分类变量的频数表和列联表，以及相应的独立性检验、相关性的度量、图形化展示结果的方法。
#我们除了使用基础安装中的函数，还将连带使用vcd包和gmodels包中的函数。
#在下面的示例中，假设A、B和c代表分类变量。
#本节中的数据来自vcd包中的Arthritis数据集。
#这份数据来自Kock&Edward（1988）.表示了一项风湿性关节炎新疗法的双盲临床实验的结果。
library(vcd)
head(Arthritis)
#7.2.1生成列联表
#生成列联表的方式有很多，表7-1中罗列的函数如下：
table(varl, 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 <- table(Arthritis$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是一个矩阵或者数据框
#要进行交叉的比变量出现在公式的右侧，~的右侧
#+作为分隔符
#若某个变量写在公式的左侧，则其为一个频数向量
#代码如下：
mytable <- xtabs(~ Treatment + Improved, data = Arthritis)
mytable
#可以使用margin.table()或prop.table()生成边际频数和比例：（这里的边际频数指的是按照行或列分组时的每一组的总频数，John注）
margin.table(mytable, 1)
prop.table(mytable, 1)
#其中1指代table()函数的第1个变量（即行变量）。每一行比例之和为1。
margin.table(mytable, 2)
prop.table(mytable, 2)
#这里的2代表table()函数的第2个变量（即列变量）。每一列比例之和为1。
#各单元格所占比例如下：
prop.table(mytable) #所有单元格比例之和为1.
#可以用addmargins()函数为这些表格添加边际和
addmargins(mytable)
addmargins(prop.table(mytable))
#在使用addmargins()函数时默认是为表中所有变量创建边际和。
#作为对照，便有以下代码
addmargins(prop.table(mytable,1), 2)
addmargins(prop.table(mytable,2) ,1)
#注意table()函数默认忽略缺失值（NA）。
#要在频数统计中将NA视为一个有效的类别，请定参数useNA="ifany"。
#创建二维列联表的第3种方法是使用gmodels包中的CrossTablel)函数。
#CrossTablel()函数仿照SAS中PROCFREO或SPSS中CROSSTABS的形式生成二维列联表。
#示例参见代码清单7-8。
install.packages("gmodels")
library(gmodels)
CrossTable(Arthritis$Treatment, Arthritis$Improved)
#CrossTable（）函数有很多选项，可以做许多事情：
#计算（行、列、单元格）的百分比；指定小数位数；
#进行卡方、Fisher和McNemar独立性检验；
#计算期望和（Pearson、标准化、调整的标准化）残差；
#将缺失值作为一种有效值；进行行和列标题的标注；生
#成SAS或SPSS风格的输出。
#更多详情，请参阅help(CrossTable）。
#3、多维列联表
#table()和xtabs(）都可以基于3个或更多的分类变量生成多维列联表 。
#函数margin.table()、prop.table()和addmargins()可以自然地推广到高于二维的情况 。
#另外，ftable()函数可以以一种紧凑且吸引人的方式输出多维列联表 。
#代码清单7-9提供了一个示例 。
mytable <- xtabs(~ Treatment + Sex + Improved, data = Arthritis)
mytable
ftable(mytable)
margin.table(mytable, 1)
margin.table(mytable, 2)
margin.table(mytable, 3)
ftable(prop.table(mytable, c(1, 2)))
ftable(addmargins(prop.table(mytable, c(1, 2)), 3))
#第①部分代码生成了三维分组各单元格的频数。
#这段代码同时演示了如何使用ftable()函数输出更为紧凑和吸引人的表格。
#第②部分代码为治疗方式（Treatment）、性别（Sex）和改善情况（Improved）生成了边际频数。
#由于使用公式~Treatement+Sex+Improve创建了这个表，所以Treatment需要通过下标1来引用，Sex通过下标2来引用，Improved通过下标3来引用。
#第③部分代码为治疗方式（Treatment）×改善情况（Improved）分组的边际频数，由不同性别（Sex）的单元加和而成。
#每个Treatment×sex组合中改善情况为None、Some和markea患者的比例由④给出。
#一般来说，比例将被添加到不在prop.table（）调用中的下标上（本例中是第3个下标，或称Improved）。
#在最后一个例子中可以看到这一点，我们在那里为第3个下标添加了边际和。
#如果想得到百分比而不是比例，可以将结果表格乘以100。
#例如语句：
ftable(addmargins(prop.table(mytable, c(1, 3)), 2)) * 100
#这里和书中代码略有不同，目的是为了和书中展示的相同（John注）
#7.2.2 独立性检验
#包括卡方独立检验。Fisher精确检验和Cochran-Mantel-Haenszel检验
#1、卡方独立检验
#使用chisq.test()函数对二维表的行变量和列变量进行卡方独立检验
#代码7-10如下
library(vcd)
mytable <- xtabs(~ Treatment + Improved, data = Arthritis)
chisq.test(mytable)
mytable <- xtabs(~ Improved + Sex, data = Arthritis)
chisq.test(mytable)
#df = 2是自由度（John注）
#书第152页有这段代码的详细解释，由于这部分很容易看懂，就不写了（John注）
#2、Fisher精确检验
#其格式为fisher.test(mytable)
mytable <- xtabs(~ Treatment + Improved, data = Arthritis)
fisher.test(mytable)
#3、Cochran-Mantel-Haenszel检验
#mantelhaen.test（）函数可用来进行Cochran-Mantel-Haenszel卡方检验，
#其原假设是，两个名义变量在第3个变量的每一层中都是条件独立的。
#下列代码可以检验治疗方式和改善情况在性别的每一水平下是否独立。
#此检验假设不存在三阶交互作用（治疗方式×改善情况×性别）。
mytable <- xtabs(~ Treatment + Improved + Sex, data = Arthritis)
mantelhaen.test(mytable)
#7.2.3 相关性度量
#上一节中的显著性检验评估了是否存在充分的证据以拒绝变量间相互独立的原假设。
#如果可以拒绝原假设，那么我们的兴趣就会自然而然地转向用以衡量相关性强弱的相关性度量。
#vcd包中的assocstats()函数可以用来计算二维列联表的phi系数、列联系数和Cramer'sV系数。
#代码清单7-11提供了一个示例。
library(vcd)
mytable <- xtabs(~ Treatment + Improved, data = Arthritis)
assocstats(mytable)
##7.2.4结果的可视化
#R中拥有的以可视化方式探索分类变量间关系的方法要远远超出其他多数统计软件。
#通常我们会使用条形图进行一维频数的可视化（参见6.1节）。
#vcd包中拥有优秀的、用于可视化多维数据集中分类变量间关系的函数，可以绘制相关图（参见11.3节）和马赛克图（参见11.4节）。
#最后，ca包中的对应分析函数允许使用多种几何表示（Nenadic&Greenacre，2007）可视地探索列联表中行和列之间的关系。
#对列联表的讨论暂时到此为止，我们将在第11章和第19章中探讨更多高级话题。
#下面我们来看看各种类型的相关系数。
####7.3 相关分析####
install.packages(c("psych", "ggm"))
#7.3.1相关的类型
#R可以计算多种相关系数，包括Pearson相关系数、Spearman相关系数、Kendall相关系数偏相关系数、多分格（polychoric）相关系数和多系列（polyserial）相关系数。
##1.Pearson、Spearman和Kendall相关系数
#Pearson积差相关系数衡量了两个定量变量之间的线性相关程度。
#Spearman等级相关系数则衡量分级定序变量之间的相关程度。
#Kendall's Tau相关系数也是一种非参数的等级相关度量。
#cor（）函数可以计算这3种相关系数，而cov（）函数可用来计算协方差。
#两个函数的参数有很多，其中与相关系数的计算有关的格式可以简化为：
cor(x, use = , method = )
#x:矩阵或数据框
#use:指定缺失数据的处理方式。
#可选值为all.obs（假设不存在缺失数据——遇到缺失数据时将报错）、
#everything（遇到缺失数据时，相关系数的计算结果将被设为missing）、complete.obs（行删除）
#以及pairwise.complete.obs（成对删除）
#method:指定相关系数的类型。可选值为Pearson、Spearman或Kendall
#默认参数为use = "everything", method = "pearson"
#代码清单7-12如下
states <- state.x77[, 1:6]
cov(states)
cor(states)
cor(states, method = "spearman")
#注意，在默认情况下得到的结果是一个方阵（所有变量之间两两计算相关）。
#我们同样可以计算非方形的相关矩阵。d代码如下：
x <- states[, c("Population", "Income", "Illiteracy", "HS Grad")]
y <- states[, c("Life Exp", "Murder")]
cor(x, y)
##2.偏相关
#偏相关分析是指在对其他一个或多个定量变量的影响进行控制的条件下分析两个定量变量之间的相互关系。
#我们可以使用ggm包中的pcor（）函数计算偏相关系数。
#ggm包没有被默认安装，在第一次使用之前需要先进行安装。
#函数调用格式为：
pcor(u, S)
#其中u是一个数值向量，前两个数值表示要计算相关系数的变量下标，其余的数值为条件变量（即要排除影响的变量）的下标。
#s为变量的协方差阵。
#这个示例有助于阐明用法：
library(ggm)
colnames(states)
pcor(c(1,5,2,3,6), cov(states))
#控制后面的列，对前两列分析偏相关值（John注）
##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-13中的示例:
cor.test(states[,3], states[,5])
#这段代码检验了预期寿命和谋杀率的Pearson相关系数为0的原假设。
#假设总体的相关度为0，则预计在一千万次中只会有少于一次的机会见到0.703这样大的样本相关度（即p=1.258e-08）。
#因为这种情况几乎不可能发生，所以我们可以拒绝原假设，从而支持了要研究的猜想，即预期寿命和谋杀率之间的总体相关度不为0。
#cor.test()函数每次只能检验一种相关关系
#psych包中的corr.tesr()函数可以一次做很多事情（同时用三种方法进行显著性检验）
library(psych)
corr.test(states, use = "complete")
#参数use可以选择pairwise或complete,是对缺失值处理方法选择成对删除还是行删除
#method 选项可以选择"pearson"（默认值）, "kendall", "spearman"
#7.3.3 相关关系的可视化
#在第11章进行讨论
####7.4 t检验####
#这里用MASS包中的UScrime数据集进行演示。（这里对t检验进行了简单的描述，就不写了，见书158页，John注）
#7.4.1 两独立样本的t检验
#调用格式如下：
t.test(y ~ x, data)
#其中y是一个数值型变量，x是一个二分类变量，
#调用格式或为：
t.test(y1, y2)
#其中，y1和y2位数值型向量（各组的结果变量）
#data的取值为一个包含了这些变量的矩阵或数据框
#注意：这里的t检验默认假定方差不相等，并使用Welsh的修正自由度
#可以添加一个参数，var.equal = TRUE以假定方差相等，并使用合并方差分析。
#默认的备择假设是双侧的，
#可以添加一个参数alternative = "less"或alternative = "greater"进行有方向的检验
#在下列代码中，我们使用了一个假设方差不相等的双侧检验，比较了南方（group1）和非南方（group0)各州的监禁概率：
install.packages("MASS")
library(MASS)
t.test(Prob ~ So, data = UScrime)
t.test(Prob ~ So, data = UScrime, var.equal = TRUE) #假定方差齐
#对于样本量小的t检验需要在之前做正态性检验和方差齐性检验，John注
#7.4.2 非独立样本的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))
#7.4.3 多于两组的情况
#第9章介绍ANOVA分析（方差分析）
####7.5 组间差异的非参数检验####
#当无法满足参数检验时，可选用非参数检验
#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="areater"进行有方向的检验。
#用Mann-Whitney U检验回答7.4节中关于监禁率的问题，将得到这些结果：
with(UScrime, by(Prob, So, median))
wilcox.test(Prob ~ So, data = UScrime)
#我们可以再次拒绝南方各州和非南方各州监禁率相同的假设（p<0.001）。
# Wilcoxon 符号秩检验是非独立样本t检验的一种非参数替代方法。
#它适用于两组成对数据和无法保证正态分布假设的情境。(不满足方差齐性和正态性，John注)
#调用格式与 Mann-Whitney U检验完全相同，不过还可以添加参数paired=TRUE。
#让我们用它解答7.4节中的失业率问题：
sapply(UScrime[c("U1", "U2")], median)
with(UScrime, wilcox.test(U1, U2, paired = TRUE))
#与t检验得出了相同的结论
####7.6 组间差异的可视化####
2596
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
