《R语言实战》day7-第七章基本统计分析_r语言am表示类别型的分组变量么?-CSDN博客

本文链接：https://blog.csdn.net/lim18719266191/article/details/132806634

本章内容

 描述性统计分析

 频数表和列联表

 相关系数和协方差

 t检验

 非参数统计

两两探索所选择变量之间的关系

导入问题

 各车型的油耗如何？特别是在对车型的调查中，每加仑汽油行驶英里数的分布是什么样的？（均值、标准差、中位数、值域等。）

 在进行新药实验后，用药组和安慰剂组的治疗结果（无改善、一定程度的改善、显著的改善）相比如何？实验参与者的性别是否对结果有影响？

 收入和预期寿命的相关性如何？它是否明显不为零？

 美国的某些地区是否更有可能因为你犯罪而将你监禁？不同地区的差别是否在统计上显著？

7.1 描述性统计分析

7.1.1 方法云集

myvars <- c("mpg", "hp", "wt") 
summary(mtcars[myvars])

summary（）：最小值、最大值、四分位数；均值（数量型）；频数统计（因子向量和逻辑型向量）
sapply（）：x是你的数据框（或矩阵），FUN为一个任意的函数。如果指定了options，它们将被传递给FUN。你可以在这里插入的典型函数有mean()、sd()、var()、min()、max()、median()、length()、range()和quantile()。
函数fivenum()可返回图基五数总括（Tukey’s five-number summary，即最小值、下四分位数、中位数、上四分位数和最大值）

sapply(x, FUN, options)

例子：通过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")
myvars
sapply(mtcars[myvars],mystats)

此时，用来计算的x是从mtcars数据中取的mpg、hp、wt列
FUN=用来计算的函数=前面定义好的mystats

7.1.2 更多方法

Hmisc包中的describe()函数可返回变量和观测的数量、缺失值和唯一值的数目、平均值、分位数，以及五个最大的值和五个最小的值。
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。

例子：

library(pastecs)
myvars<-c("mpg","hp","wt")
stat.desc(mtcars[myvars])

library(psych) 
myvars <- c("mpg", "hp", "wt") 
describe(mtcars[myvars])

stat.desc的结果

psych的包下载得有问题，再看看

7.1.3 分组计算描述性统计量

1. 用aggregate（）

myvars <- c("mpg", "hp", "wt") 
aggregate(mtcars[myvars], by=list(am=mtcars$am), mean) 
aggregate(mtcars[myvars], by=list(am=mtcars$am), sd)

注意list(am=mtcars$am)的使用。如果使用的是list(mtcars$am)，则am列将被标注为 Group.1而不是am。你使用这个赋值指定了一个更有帮助的列标签。
如果有多个分组变量，可以使用by=list(name1=groupvar1, name2=groupvar2, ... , nameN=groupvarN)这样的语句。遗憾的是，aggregate()仅允许在每次调用中使用平均数、标准差这样的单返回值函数。 它无法一次返回若干个统计量。

2.用by（）：可以同时返回若干个

by(data, INDICES, FUN)

其中data是一个数据框或矩阵，
INDICES是一个因子或因子组成的列表，定义了分组，
FUN是任意函数。

例子

dstats <- function(x)sapply(x, mystats) 
myvars <- c("mpg", "hp", "wt") 
by(mtcars[myvars], mtcars$am, dstats)

7.1.4 分组计算的拓展：doBy中的summaryBy（）

summaryBy(formula, data=dataframe, FUN=function)

其中的formula接受以下的格式： var1 + var2 + var3 + ... + varN ~ groupvar1 + groupvar2 + ... + groupvarN 在~左侧的变量是需要分析的数值型变量，而右侧的变量是类别型的分组变量。
function可为任何内建或用户自编的R函数。

自编函数：比如

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))
}

例子：用doBy包中的summaryBy()分组计算概述统计量

library(doBy) 
summaryBy(mpg+hp+wt~am, data=mtcars, FUN=mystats)

按am来分类，data=mtcars，使用的函数FUN=前面自设的mystats

7.2 频数表和列联表

7.2.1 生成频数表

1. 一维列联表

table()函数生成简单的频数统计表
- table()函数默认忽略缺失值（NA）。要在频数统计中将NA视为一个有效的类别，请设定参数useNA="ifany"。
用prop.table()将这些频数转化为比例值
用prop.table()将这些频数转化为比例值

2. 二维列联表

1）mytable <- table(A, B)
2）xtabs()函数还可使用公式风格的输入创建列联表，
- 格式为： mytable <- xtabs(~ A + B, data=mydata)
- 总的来说，要进行交叉分类的变量应出现在公式的右侧（即 ~符号的右方），以+作为分隔符。若某个变量写在公式的左侧，则其为一个频数向量（在数据已经被表格化时很有用）

mytable <- xtabs(~ Treatment+Improved, data=Arthritis) 
mytable

你可以使用margin.table()和prop.table()函数分别生成边际频数和比例。行和与行比例可以这样计算：
- margin.table(mytable, 1)
- prop.table(mytable, 1)
- 下标1指代table()语句中的第一个变量

addmargins()函数为这些表格添加边际和
- 默认是对所有变量（如图）
- addmargins(prop.table(mytable, 1), 2) 则仅为行
- addmargins(prop.table(mytable, 2), 1) 则仅为列

3）使用gmodels包中的CrossTable()函数

library(gmodels) 
CrossTable(Arthritis$Treatment, Arthritis$Improved)

CrossTable()函数有很多选项，可以做许多事情：

计算（行、列、单元格）的百分比；
指定小数位数；
进行卡方、Fisher和McNemar独立性检验；
计算期望和（皮尔逊、标准化、调整的标准化）残差；
将缺失值作为一种有效值；
进行行和列标题的标注；
生成SAS或SPSS风格的输出。
参阅help(CrossTable)以了解详情。

3. 多维列联表

table() 和 xtabs() 都可以基于三个或更多的类别型变量生成多维列联表。
margin.table()、prop.table()和addmargins()函数可以自然地推广到高于二维的情况。
ftable()函数可以以一种紧凑而吸引人的方式输出多维列联表。

mytable <- xtabs(~ Treatment+Sex+Improved, data=Arthritis) 
mytable 

 ftable(mytable)

margin.table(mytable, 1) # 边际频数

margin.table(mytable, c(1, 3)) #治疗情况*改善情况的边际频数

# 这里的c（1，3）就是第1、3个变量的意思

ftable(prop.table(mytable, c(1, 2)))
# 治疗情况（Treatment）× 性别（Sex）的各类改善情况比例

ftable(addmargins(prop.table(mytable, c(1, 2)), 3)) 
# 加了边际和的情况

# 如果想得到百分比而不是比例，可以将结果表格乘以100。例如：
ftable(addmargins(prop.table(mytable, c(1, 2)), 3)) * 100

7.2.2 独立性检验

卡方检验
fisher检验
Cochran-Mantel-Haenszel

1. 卡方独立性：chisq.test()

library(vcd) 
mytable <- xtabs(~Treatment+Improved, data=Arthritis) 
chisq.test(mytable)

p＜0.01：结果显示不独立

这里的p值表示从总体中抽取的样本行变量与列变量是相互独立的概率。由于概率值很小，所以你拒绝了治疗类型和治疗结果相互独立的原假设

mytable <- xtabs(~Improved+Sex, data=Arthritis) 
chisq.test(mytable)

p>0.05 独立
产生警告信息的原因是，表中的6个单元格之一（男性－一定程度上的改善）有一个小于5的值，这可能会使卡方近似无效。

2. Fisher精准检验：fisher.test()

Fisher精确检验的原假设是：边界固定的列联表中行和列是相互独立的。
调用格式为fisher.test(mytable)，其中的mytable是一个二维列联表
这里的fisher.test()函数可以在任意行列数大于等于2的二维列联表上使用，但不能用于2×2的列联表。

3. Cochran-Mantel-Haenszel检验：mantelhaen.test()

其原假设是，两个名义变量在第三个变量的每一层中都是条件独立的。
下列代码可以检验治疗情况和改善情况在性别的每一水平下是否独立。
此检验假设不存在三阶交互作用（治疗情况×改善情况×性别）。

mytable <- xtabs(~Treatment+Improved+Sex, data=Arthritis) 
mantelhaen.test(mytable)

7.2.3 相关性的度量

例子：二维列联表的相关性度量

library(vcd) 
mytable <- xtabs(~Treatment+Improved, data=Arthritis) 
assocstats(mytable)

vcd包也提供了一个kappa()函数，可以计算混淆矩阵的Cohen’s kappa值以及加权的kappa值。
（举例来说，混淆矩阵可以表示两位评判者对于一系列对象进行分类所得结果的一致程度。

7.3 相关

7.3.1 相关的类型

1. Pearson、Spearman、Kendall相关

cor()函数可以计算这三种相关系数，而cov()函数可用来计算协方差。
两个函数的参数有很多，其中与相关系数的计算有关的参数可以简化为：
- cor(x, use= , method= )

如果想选部分，计算非方形的相关矩阵，可用：

x <- states[,c("Population", "Income", "Illiteracy", "HS Grad")] 
y <- states[,c("Life Exp", "Murder")] 
cor(x,y)

上述结果并未指明相关系数是否显著不为0（即，根据样本数据是否有足够的证据得出总体相关系数不为0的结论）。由于这个原因，你需要对相关系数进行显著性检验

2. 偏相关：偏相关是指在控制一个或多个定量变量时，另外两个定量变量之间的相互关系。

你可以使用 ggm包中的pcor()函数计算偏相关系数。调用格式：

pcor(u, S)

其中的u是一个数值向量，前两个数值表示要计算相关系数的变量下标，其余的数值为条件变量（即要排除影响的变量）的下标。
S为变量的协方差阵。

library(ggm) 
colnames(states)
[1] "Population" "Income" "Illiteracy" "Life Exp" "Murder" "HS Grad" 
pcor(c(1,5,2,3,6), cov(states)) # 在控制了收入、文盲率和高中毕业率的影响时，人口和谋杀率之间的相关系数为0.346。

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）。

 cor.test(states[,3], states[,5])

遗憾的是，cor.test()每次只能检验一种相关关系。但幸运的是，psych包中提供的 corr.test()函数可以一次做更多事情。corr.test()函数可以为Pearson、Spearman或Kendall 相关计算相关矩阵和显著性水平。

library(psych) 
corr.test(states, use="complete")

参数use=的取值可为"pairwise"或"complete"（分别表示对缺失值执行成对删除或行删除）。
参数method=的取值可为"pearson"（默认值）、"spearman"或"kendall"。
这里可以看到，人口数量和高中毕业率的相关系数（–0.10）并不显著地不为0（p=0.5）

其他显著性检验：在多元正态性的假设下，psych包中的pcor.test() 函数①可以用来检验在控制一个或多个额外变量时两个变量之间的条件独立性。

pcor.test(r, q, n)

其中的r是由pcor()函数计算得到的偏相关系数，q为要控制的变量数（以数值表示位置），n为样本大小。

psych包中的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)

调用格式或为：

t.test(y1, y2)

其中的y1和y2为数值型向量（即各组的结果变量）。
可选参数data的取值为一个包含了这些变量的矩阵或数据框。
与其他多数统计软件不同的是，这里的t检验默认假定方差不相等，并使用Welsh 的修正自由度。
你可以添加一个参数var.equal=TRUE以假定方差相等，并使用合并方差估计。
默认的备择假设是双侧的（即均值不相等，但大小的方向不确定）。
你可以添加一个参数 alternative="less"或alternative="greater"来进行有方向的检验。

# 假设方差不等的双侧检验
library(MASS) 
t.test(Prob ~ So, data=UScrime)

由于结果变量是一个比例值，你可以在执行t检验之前尝试对其进行正态化变换。在本例中，所有对结果变量合适的变换（Y/1-Y、log(Y/1-Y)、arcsin(Y)、arcsin(sqrt(Y)）都会将检验引向相同的结论。

7.4.2 非独立样本的t检验

较年轻（14~24岁）男性的失业率是否比年长（35~39岁）男性的失业率更高？在这种情况下，这两组数据并不独立。你不能说亚拉巴马州的年轻男性和年长男性的失业率之间没有关系。

在两组的观测之间相关时，你获得的是一个非独立组设计（dependent groups design）。前-后测设计（pre-post design）或重复测量设计（repeated measures design）同样也会产生非独立的组。
非独立样本的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))

重点关注：mean of the differences；如果它足够大，可以保证拒绝年长的年轻男性的平均失业率相同的假设。

7.4.3 多于两组的情况：方差分析（ANOVA，需要能够假设数据是从正态总体中独立抽样而得——第九章

7.5 组间差异的非参数检验：数据无法满足t检验或ANOVA的参数假设

7.5.1 两组比较

1. 若两组数据独立，可以使用Wilcoxon秩和检验来评估观测是否是从相同的概率分布中抽得的

（更广为人知的名字是Mann-Whitney U检验）（即，在一个总体中获得更高得分的概率是否比另一个总体要大）。

调用格式为：

wilcox.test(y ~ x, data)

其中的y是数值型变量，而x是一个二分变量。

wilcox.test(y1, y2)

其中的y1和y2为各组的结果变量。

可选参数data的取值为一个包含了这些变量的矩阵或数据框。默认进行一个双侧检验。你可以添加参数exact来进行精确检验，指定alternative="less" 或alternative="greater"进行有方向的检验。

 with(UScrime, by(Prob, So, median)) 
 wilcox.test(Prob ~ So, data=UScrime)

你可以再次拒绝南方各州和非南方各州监禁率相同的假设（p<0.001）。

2. Wilcoxon符号秩检验是非独立样本t检验的一种非参数替代方法。它适用于两组成对数据和无法保证正态性假设的情境。

调用格式与Mann-Whitney U检验完全相同，不过还可以添加参数 paired=TRUE。

sapply(UScrime[c("U1","U2")], median) 
with(UScrime, wilcox.test(U1, U2, paired=TRUE))

当t检验的假设合理时，参数检验的功效更强（更容易发现存在的差异）。而非参数检验在假设非常不合理时（如对于等级有序数据）更适用。

7.5.2 多于两组的比较

Q：想比较美国四个地区（东北部、南部、中北部和西部）的文盲率，应该怎么做呢？

单向设计（one-way design），我们可以使用参数或非参数的方法来解决这个问题。

如果无法满足ANOVA设计的假设，那么可以使用非参数方法来评估组间的差异。
如果各组独立，则Kruskal-Wallis检验将是一种实用的方法。

kruskal.test(y ~ A, data)

其中的y是一个数值型结果变量，A是一个拥有两个或更多水平的分组变量（grouping variable）。

如果各组不独立（如重复测量设计或随机区组设计），那么Friedman检验会更合适。

friedman.test(y ~ A | B, data)

例子

states <- data.frame(state.region, state.x77)
kruskal.test(Illiteracy ~ state.region, data=states)

显著性检验的结果意味着美国四个地区的文盲率各不相同（p<0.001）。虽然你可以拒绝不存在差异的原假设，但这个检验并没有告诉你哪些地区显著地与其他地区不同。
要回答这个问题，你可以使用Wilcoxon检验每次比较两组数据。
一种更为优雅的方法是在控制犯第一类错误的概率（发现一个事实上并不存在的差异的概率）的前提下，执行可以同步进行的多组比较，这样可以直接完成所有组之间的成对比较。
- 用wmc（）
  - 每次用Wilcoxon检验比较两组，并通过p.adj()函数调整概率值。

例子

source("http://www.statmethods.net/RiA/wmc.txt") 
# source()函数下载并执行了定义wmc()函数的R脚本
states <- data.frame(state.region, state.x77) 

wmc(Illiteracy ~ state.region, data=states, method="holm")

函数的形式是wmc(y ~ A, data, method)，其中y是数值输出变量，A是分组变量，data是包含这些变量的数据框，method指定限制I类误差的方法。代码清单7.17使用的是基于Holm（1979）提出的调整方法，它可以很大程度地控制总体I类误差率（在一组成对比较中犯一次或多次I类错误的概率）。参阅help(p.adjust) 以查看其他可供选择的方法描述。