将军笔记之基本统计
八一建军节是中国人民解放军建军纪念日,每年的八月一日举行,由中国人民革命军事委员会设立,为纪念中国工农红军成立的节日。
值此佳节,则能不心情澎湃,于是,我决定好好学生信。可是,我想有一手的好资料啊,作为学生党,我一穷二白,三丑四没基础。
怎么办呢?嗷~~某天在生信菜鸟团的群里,发现了这个人。他给我分享了好多实用资料,然后,我自学,然后就被大神看上,当了小编。自此就可以直接请教大神问题,感觉自己起飞了。
扫他!
扫他!
扫他!
说正事!
1. 描述性统计分析:
使用的数据:汽车的mpg(每加仑汽油行驶英里数),hp(马力),wt(车重):
左边就是我们要使用的数据的前几行。总共32行。获取过程没啥好说的。
各种描述统计方法
summary()函数:
summary()函数给出了最大最小值,四分位数以及均值;
如果是逻辑向量或因子向量,summary()函数会给出频数统计。
使用sapply函数计算任意描述性统计量:sapply(x,FUN,options)。x是数据框或矩阵,FUN为一任意函数,options被传递给fun。此处FUN可以是mean, sd, var, min, max, median, length, range, quantile, fivenum()(fivenum返回图基五数总括(Tukey's five-number summary,音译的啊?))
基础安装中没有偏度和峰度的计算函数。下面是自己写的一个包含偏度和峰度的函数,并用sapply函数作用于各个变量:(sapply这不是用于数据框和矩阵的么,怎么前面说是用于列表的。只有lapply是用于列表的吧?)
函数的编辑部分没啥,最难的部分就是偏度和峰度的公式了;
注意R计算时的广播功能:x-m时,返回的是x中的每一项减m后的向量。
注意sapply的作用和结果:把参数中的函数应用于各列。
Hmisc包中的describe()函数可以对数据进行描述:
载入Hmisc包前需要载入其他几个包,按要求操作即可。(其实不用的,提示的意思是说R已经自动载入了需要的包,不需要手动载入。好傻啊...)
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值计算平均数的置信区间。
载入pastecs包之前需要载入boot包。
函数调用如上所述,没啥可说的。
psych包里也有一个名为describe()的函数,计算非缺失值的数量,均值,标准差,中位数,截尾均值,绝对中位数,最小值,最大值,值域,偏度,峰度,均值的标准误:
已经载入的其他包中含有与psych包中名字相同的函数,所以载入psych包后,其他包的同名函数被屏蔽掉,只保留最后一个载入的包中的函数名。
如果想使用之前载入的包中的函数,需要指定包的名字:Hmisc::describe(mtcars[vars])
分组计算描述性统计量
复习:aggregate()函数的用法:
选定数据集,指定要折叠的变量(就是根据其分组的变量),以及指定应用于分组的函数。要折叠的变量用列表来给出,可以指定折叠后的名字。自行联想python中字典的键和值。这种方法只能返回一个统计量。
要想一次对分组返回多个统计量,要是用by()函数,格式为:by(data,INDECES,FUN),data是数据框或矩阵,INDECES是因子或因子组成的列表,用于定义分组;FUN是任意函数。
按书上的代码写,跑不通。原因不明。第一行函数的定义就很奇怪,没有return(),用的也不是花括号。
doBy包和psych包也提供了分组计算描述性统计量的函数。doBy中的summaryBy()的使用格式是:summaryBy(formula, data=dataframe,FUN=function), fomular格式为:var1+var2+...+varn~groupvar1+groupvar2+...+groupvarN,var是需要分析的数值型变量,groupvar是类别型分组变量。
psych包中的describe.by()函数可计算和describe相同的描述性统计量:
describe.by()函数不允许指定任意函数。
reshape包用于分组描述:
这组语句能正常运行,看来函数定义没问题。
melt函数中,measure.vars=表示对这些数值型数据做描述,id.vars指定分组变量(分别指定测量变量和分组变量)
cast函数中variable仅表示重铸后数据框中的变量variable。melt后得到的数据框中有一个变量名为variable:
cast函数中的点号是什么鬼...(表示除了~左边的变量以外剩下的其他所有变量)
2. 频数表和列联表
使用的数据如下所示:
一维列联表
注意with函数的用法,复习一下。
可以使用prop.table函数将频数转化为比例值:
注意prop.table()函数的参数是一维列联表。这个函数再乘以100,就得到了百分比。
二维列联表
如果要用table函数创建二维列联表,则table(A,B),A是行变量,B是列变量
xtabs函数:mytable <-xtabs(~A+B, data=mydata),A和B是要交叉的变量,~左侧为一频数向量,在数据已经被表格化时有用。
生成Treatment和Improved两个变量的交叉列联表;生成的表,行为Treatment,列为Improved;
margin.table()函数用于生成行/列汇总;第二个参数等于1时,返回行(第一个变量)的和;第二个参数等于2时,返回列(第二个变量)的和
prop.table()函数返回比值。第二个参数等于1时返回的比值以行之和为总体;第二个参数等于2时,返回的比例以列之和为总体。
如果不加第二个参数,则返回的比例以所有值的和为总体。
使用addmargins()函数给各行各列加汇总和,第二个参数=1增加行,第二个参数=2增加列
(注意这里的“增加行”“增加列”指的是“增加一行”“增加一列”的意思,而不是“增加行汇总”“增加列汇总”的意思)
对比例二维表也是一样使用addmargins()加汇总,行汇总列汇总都可以,行列都汇总也可以,把所有单元格之和作为整体1。
但是注意:addmargins(prop.table (mytable,1),1)和addmargins (prop.table(mytable, 2),2)这两个表述是没有意义的。
table()函数默认忽略NA值。如果要不忽略NA,使用参数useNA=“ifany”。
还可以使用gmodels包中的CrossTable()函数创建二维列联表。如下所示:
并不是完全理解输出的都是什么...CrossTable()函数可以:计算行列单元格百分比,指定小数位数,卡方、Fisher、McNemar独立性检验,计算期望和皮尔逊、标准化、调整的标准化残差,将缺失值作为有效值,标注行列标题等。现用现查吧。
多维列联表
前面用于创建和加工二维列联表的函数:table(), xtabs(), margin.table(), prop.table(), addmargins()函数都可以推广到多维列联表。另外ftable()函数可以以紧凑的形式输出多维列联表。例子:
xtabs()生成多维列联表的方式是这样的。生成的表示这样的。
ftable()函数输出的多维列联表是这
样的,我觉得这个看着紧凑点。注意ftable()函数不是生成了多维列联表,而是把已有的多维列联表以一种更加紧凑的方式输出了。
使用margin.table()函数分别对多维列联表的三个维度(变量)进行汇总。也可以对任意两个或几个变量进行汇总统计。
此时多维列联表坍缩为低维列联表。
先用prop.table()函数返回各组频数占总体的比例,总体是根据第一第二个变量分组后每组的总和。然后用ftable()函数输出成紧凑点的形式。
这个函数嵌套的层数比较多:先是用prop.table()函数把mytable变成比例形式,分组同上;然后用addmargins()函数添加汇总,汇总是根据第三个变量汇总的;然后用ftable()函数输出成紧凑的形式。
独立性检验
卡方独立性检验:对二维表的行变量和列变量进行,chisq.test(mytable)
没啥可说的。
Fisher精确检验:原假设:边界固定的列联表中的行和列是相互独立的。
书上说:fisher.test()函数可以在任意行列数大于等于2的二维列联表上使用,但不能用于2*2列联表。不知道是什么意思?
Cochran-Mantel-Haenszel检验:原假设:两个名义变量在第三个变量的每一层中都是条件独立的(哪两个变量?)
用法没什么可说的,意义比较重要。根据原假设,这里p值这么小,表示的意义应该是:两个名义变量并不是在第三个变量的每一层中都独立,亦即,至少在第三个变量的一层上,两个名义变量不独立。(不能指定哪两个变量是前两个哪一个变量是第三个吗?)
相关性度量
如果经过上步检验,变量不独立,则评估变量相关性。vcd包中的assocstats()函数可以用来计算二维列联表的phi系数,列联系数和Cramer‘s V系数(这都是些啥...)。
值越大,相关性越强。vcd包提供了一个kappa()函数,可以计算混淆矩阵的Cohen's kappa值以及加权的kappa值(混淆矩阵可以表示两位评判者对于一系列对象进行分类所得结果的一致程度)(这都说了些啥???)
结果可视化:
条形图:一维频数可视化;
vcd包:可视化多维数据集 中 类别型变量间关系的函数,可以绘制马赛克图和关联图;
car包中的对应分析函数允许使用多种几何表示可视地探索列联表中行和列之间的关系。
将表转化为扁平格式:
用下面这个函数将列联表转为原始数据:
> table2flat <- function(mytalbe){
+ df <- as.data.frame(mytable)
+ row <- dim(df)[1]
+ cols <- dim(df)[2]
+ x <- NULL
+ for (i in 1:row){
+ for (j in 1:df$Freq[i]){
+ rows <- df[i,c(1:(cols-1))]
+ x <- rbind(x,rows)
+ }
+ }
+ row.names(x) <- c(1:dim(x)[1])
+ return(x)
+ }
此函数可以接受一个表格并返回一个扁平格式的数据框。
比如遇到这样的表格:
无改善 一定程度改善 显著改善
安慰剂 29 7 7
药物 13 17 21
要将上述表格扁平化,代码如下:
费这么大劲,其实所做的工作就是把各种变量组合起来,然后看每种组合的出现次数是多少,是多少,组合变量就出现多少次。
Pearson积差相关系数衡量了两个定量变量之间的线性相关程度;Spearman等级相关系数则衡量分级定序变量之间的相关程度;Kendall’s Tau相关系数也是非参数的等级相关度量。cor()函数可以计算这三种相关系数,格式cor(x,use= , method= ),其中x是矩阵或数据框,use用于指定缺失数据的处理方式,有all.obs(假设不存在缺失数据,存在时将报错),everything(遇到缺失数据时,相关系数的计算结果是missing),complete.obs(行删除),pairwise.complete.obs(成对删除)。method指定相关系数类型,有pearson, spearmen和kendall。默认参数是use="everything", method="pearson"。
cov()函数可用来计算协方差。示例如下:
cov()函数计算了方差和协方差。
cor(states)计算了Pearson积差相关系数
cor(states,method="spearman")函数计算了Spearman等级相关系数。
上述函数得到的是两两相关的方阵。也可以指定待分析变量:(原来是通过控制数据来控制的计算哪些变量之间的相关性,而不是通过函数本身控制的)
首先提取states数据集中的变量,然后对两个“变量集”做相关。注意,此时州名是rowname,自动出现在提取的“变量集”中,
偏相关:在控制一个或多个定量变量时,另外两个定量变量之间的相互关系。使用ggm包中的pcor()函数计算偏相关系数。函数调用格式:pcor(u,S)。u是数值向量,前两个数值表示要计算相关系数的变量下标,其余的数值为条件变量下标(就是要控制的变量的下标);S为变量的协方差阵(协方差阵可以用cov()函数计算)。例子:
在控制了收入、文盲率、高中毕业率影响后,人口和谋杀率之间的相关系数为0.346。
polycor包中的hetcor()函数可以计算一种混合的相关矩阵,包括Pearson积差相关系数,数值型变量和有序变量之间的多系列相关系数,有序变量之间的多分格相关系数以及二分变量之间的四份相关系数。多系列、多分格和四份相关系数都假设有序变量或二分变量有潜在的正态分布导出。
相关性的显著性检验:原假设为变量间不相关(总体相关系数为0),使用cor.test()函数对单个的相关系数进行检验。cor.test(x,y,alternative= ,method= ),其中xy是要检验相关性的变量,alternative指定双侧或单侧检验,取值为"two.side", "less", "greater";method指定要计算的相关类型,“pearson”, "spearman", "kendall"。例子:
cor.test()每次只能检验一种 相关关系。psych包中提供的corr.test()函数可以计算相关矩阵和显著性水平。例子:
use=的取值可为pairwise或complete,分别表示对缺失值执行成对删除或行删除。method=可以是pearson,spearman,kendall。
ggm包中的pcor.test()函数可以用来检验在控制一个或多个额外变量时两个变量之间的条件独立性。使用格式:pcor.test(r, q, n),其中r是pcro()函数计算得到的偏相关系数,q为要控制的变量(用数值表示位置),n为样本大小。
psych包中的r.test()函数提供了多种显著性检验方法,用到再查。
t检验:数据要服从正态分布
独立样本t检验
t.test(y~x, data),y是数值变量,x是二分变量;
t.test(y1, y2),y1,y2均为数值向量
上述函数默认方差不相等,并使用Welsh的修正自由度。添加参数var.equal=TRUE指定方差相等,并使用合并方差估计。默认的备择假设是双侧的,可以使用参数alternative=“less”或“greater” 进行单侧检验。例子:
非独立样本t检验:假定组间差异服从正态分布
对非独立样本做检验时,只需在t.test()函数中加入选项paired=TRUE即可:
组间差异的非参数检验:数据存在严重偏倚,或呈现有序关系时使用。
若两组数据独立,可以使用Wilcoxon秩和检验(其实就是Mann-Whitney U检验)来判断两样本是否来自同一总体。调用格式:wilcox.test(y~x,data),或者wilcox.test(y1,y2)。默认双侧检验,alternative=“less”或“greater”进行单侧检验;添加exact参数进行精确检验。
复习一下by()函数的用法:第一个参数:要处理的变量;第二个参数,要分组的变量;第三个参数,用于处理的函数。
wilcox符号秩检验用于非独立样本,适用于两组成对数据和无法保证正态性假设的数据,只需在wilcox.test()函数中添加paired=TRUE选项即可:(wilcoxon秩和检验用于独立样本的非参数检验;wilcoxon符号秩检验用于非独立样本的非参数检验)
多组间的比较:比较不同组间数值差异,称为单向设计(one-way design)。如果无法满足ANOVA设计的假设(数据不是从正态总体中独立抽样而得),那么使用非参方法。
各组独立:Kruskal-Wallis检验:kruskal.test(y~A, data),y是数值型变量,A是分组变量
首先把states所属地区的信息加入到各州信息中,生成新的数据集states;
然后根据state.region分组统计不同地区文盲率差异是否显著。
npmc包中的npmc()函数可以进行多组比较:npmc(var,class):然而npmc包无法在最新版本的R中使用了。网上有修改方法(http://www.w2bc.com/article/52835),不过,算了吧,先不弄了。
各组不独立:Friedman检验:friedman.test(y~A|B,data),y是数值变量,A是分组变量,B是一个用以认定匹配观测的区组变量(blocking variable)。这部分书上没有再讲。
(多组比较检验,样本独立使用kruskal.test(),样本不独立(区组设计?)则使用Friedman检验)
小结一下:
做描述统计的几个函数:
summary(dataframe);
sapply(data, FUN, options) data是数据框或矩阵,FUN函数作用于各列,options是函数的选项;
Hmisc::discribe(data)
pastecs::stat.desc(data,basic=TRUE,desc=TRUE,norm=FALSE,p=0.95)
如果basic=TRUE,计算所有值、空值、缺失值的数量,最小值,最大值,值域,总和;
desc=TRUE,计算中位数,平均值,均值的标准误,均值95%置信区间,方差,标准差,变异系数;
norm=TRUE,则返回正态分布统计量,偏度峰度及其统计显著程度,shapiro-wilk正态检验;
p值计算平均数的置信区间。
psych::describe(data)
aggregate(data, by=list(), FUN)
by(data, group, FUN) group: 分类变量
doBy::summaryBy(v1+v2+v3~g1+g2+g3, data=, FUN=) v, 数值变量;g, 分组变量;
psych::describe.by(var, group) 返回一些统计量;
d <- reshape::melt(data, measure.vars=, id.vars=) measure.vars, 数值变量;id.vars, 分类变量。id.vars组合成一个个唯一的标签;然后一列是variable变量,measure.vars中的变量名就是variable中的值;再后面一列就是前面三列限定下的值;reshape::cast(d, v1+v2~., FUN) d就是melt得到的数据框;v1, v2行变量,其余(.)列变量,用FUN函数处理(v1v2限定出的样本肯定不止一个;多个样本值使用FUN函数处理)
列联表:
table(A,B) 两个变量产生列联表;
xtabs(~A+B, data=dataframe) 产生列联表
prop.talbe(table, n) 将列联表转为比例;n=1时以行为100%;n=2以列为100%;没有n时以总数为100%;
margin.talbe(table, n) 返回对table的汇总,n=1时行汇总;n=2列汇总;没有n时对行列都汇总;
addmargins(table, n) 对table返回汇总。这个不仅仅返回汇总,连表本身也返回。
ftable(table) table由xtabs(~V1+V2)产生,将三维列表紧凑输出;
chisq.test(table) 卡方检验
fisher.test(table) Fisher检验;
mantelhaen.test(table) 两个名义变量在第三个变量的每层都是条件独立的;
相关性:
cor(data,use=, method=) use=all.obs(假设不存在缺失数据,存在时将报错)/everything(遇到缺失数据时,相关系数的计算结果是missing)/complete.obs(行删除)/pairwise.complete.obs(成对删除)。method=pearson/spearmen/kendall。默认参数是use="everything", method="pearson"。
ggm::pcor(c(n1,n2,n3...),cov(dataframe))在控制一个或多个定量变量时,另外两个定量变量之间的相互关系。前两个数值表示要计算相关系数的变量的下标(在dataframe中),其余的数值为条件变量下标(就是要控制的变量的下标);cov(dataframe)为变量的协方差阵(协方差阵可以用cov()函数计算);
polycor::hetcor()
cor.test(x, y, alternative=, method=) alternative="two.side"/"less"/"greater";method=“pearson”/"spearman"/"kendall"。
psych::corr.test(dataframe, use=, method=) use="complete"(行删除)/"pairwise"(成对删除)。method=pearson/spearman/kendall。
ggm::pcor.test()
psych::r.test()
样本是否来源于同一总体检验
t.test(y~x, data) x二分类变量;
t.test(y1, y2, var.equal=TRUE, alternative=, paired=TRUE) var.equal=TRUE指定方差相等;alternative="less"/"greater"; paired=TRUE配对;
wilcox.test(y~x, data)
wilcox.test(y1, y2, alternative="less"/"greater", exact=TRUE)
kruskal.test(y~x, data) 多组独立样本检验;
friedman.test(y~A|B, data) 区组设计,B是区组;
统计知识需要加强!
快快留言!!!
往期精彩回顾