单变量数据分析简介

单变量数据分析简介

Posted by omicsPie on 2018-06-07

英国著名物理学家, 原子核物理学之父欧内斯特·卢瑟福曾经说过:“如果你的实验结果需要一个统计员,那你就需要重新设计一个实验了。”but啊,今时不同往日,我们现今的实验检测技术精度越来越高,微小的数据差异都能被捕捉到,而这些数据差异的规律也就很难仅仅通过简单计算与“普通观察法”被发现。So,现今实验不仅仅是需要统计学,简直是与统计相爱相杀。

啰嗦这么多,小编主要是想唤起实验小姐姐们对统计(做统计的我)的重视,既然目的已达到,那我们就进入今天的主题,简单介绍一下单变量的统计分析方法。

统计分析包括统计描述和统计推断两部分,统计描述是指通过对样本数据进行特征描述来反应总体的特征,描述方式主要包括绘制统计图表、计算数据分布特征(如均值、众数、中位数)等方式。

在R中随机生成20个正态分布的数据,数据服从N(0,22)

1
data <- rnorm(20,0,2^2)

要计算数据data常用的数据分布特征可以用如下几个语句:

1
2
3
4
5
6
mean(data)##data的均值
median(data)##data的中位数
sd(data)##data的标准差
var(data)##data的方差
quantile(data,c(0.1,0.25,0.8))
##data10%,25%,80%,分位数

要绘制反应数据特征的图表常用如下几个语句:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
hist(data,breaks = 8,

col = "pink",

main = "直方图",

xlab = "分组",

ylab = "频数")

##直方图



plot(density(data), main = "核密度图")
polygon(density(data), col = "pink", border = "green")
rug(data,col = "red")
##核密度图
boxplot(data, main = "箱线图", col = "pink",border = "green")
##箱线图

统计描述是统计推断的基础,描述之后就是统计推断了,单变量推论统计的目的一般是通过样本数据对总体的状况进行推断,主要从参数估计和假设检验两个方面进行。

参数估计通过计算样本数据的参数来估计总体的参数,又分为点估计和区间估计两种方法。点估计就是直接将样本参数作为总体参数的估计值,而区间估计是按照预先给定的概率 所计算出来的一个包含总体参数的区间。点估计的R函数即上述统计描述中计算均值、中位数等所用的语句,区间估计则要根据所估计的统计量的分布来求相应的分位数,例如最简单的正态分布数据,在置信水平下,其均值的置信区间是

其中x(上面有划线)是样本均值,n是样本量。

假设检验通过检验样本数据是否符合预先假设的分布来推断总体是否符合预先假设的分布,常用的假设检验分参数检验和非参数检验。参数检验是通过检验样本参数所代表的未知总体的参数和已知总体的参数是否有差异,有点绕哈,其实就是想检验我们所要研究的总体是不是和我们预先假设的分布是同一个分布,常用的参数检验有t检验,ANOVA(方差分析)等。因为参数检对数据要求较高,一般要求样本数据服从正态分布且其所代表的总体间方差齐,当不符合要求时,如果仍然使用参数检验,其统计推断的结果就不够可信了,此时,我们就需要用到非参数检验,常用的非参数检验有Wilcoxon秩和检验,Kruskal-Wallis检验等。

现在在R中随机生成3组50个正态分布的数据

1
2
3
data1 <- rnorm(50,0,2^2)
data2 <- rnorm(50,2,2^2)
data3 <- rnorm(50,4,2^2)

t检验可以完成独立样本均数的检验及两个样本均数的检验,如当我们直接运行语句t.test(data1)时,就是检验data1的均数(data1的均数代表我们所研究的总体的均数)和假设总体的均数0之间是否有差异,这里置信水平默认0.95,也可以通过函数中的conf.level参数设置。而当我们运行t.test(data1,data2)时,就是检验data1的均数和data2的均数之间是否有差异,也就是检验两个数据所代表的总体均数是否有差异。 t检验适用的是两个样本之间的比较,而在日常实验中,实验往往会有多个分组,因此要比较不同分组之间的数据结果是否有差异就要用到方差分析。这里,我们用aov()函数来比较data1,data2和data3这3组数据之间是否有差异。

1
2
3
4
5
6
7
8
9
10
11
12
13
all.data=c(data1,data2,data3)
all.data <- c(data1,
data2,
data3)
#将数据组合为长向量
group.name <- c(
"data1",
"data2",
"data3")
group <- rep(groupname,each = 50)
##设置每个数据对应的分组
result <- aov(all.data~group,data=NULL)
##数据代入函数进行检验

类似参数检验,非参数检验中Wilcoxon秩和检验主要用于两组之间的比较,而Kruskal-Wallis检验则适用于多组的比较。例如,若要检验data1和data2是否是从具有相同概率分布的总体中抽样所得,可以调用函数wilcox.test(data1,data2)。若要比较data1,data2和data3这3组数据是否是从具有相同概率分布的总体中抽样所得,其函数调用和aov()类似,如下所示:

1
2
3
4
5
6
7
8
all.data=c(data1,data2,data3)
##将数据组合为长向量
groupname <- c("data1","data2","data3")
group <- rep(groupname,each = 50)
##设置每个数据对应的分组alldata <- data.frame(group,
all.data=c(data1,data2,data3))
result <- kruskal.test(all.data~group,data= alldata)
##数据代入函数

进行检验
统计分析博大精深,单单是上述的ANOVA就要cover许许多多实验设计,而不同的实验设计其数据分析也是不尽相同,这里小编只肤浅地介绍了一些最简单最常用的单变量分析及其相应的R函数调用,后续的文章也会陆续深入介绍每一种分析方法,希望能对大家有所帮助。

1. 数据结构

向量 (vector):R语言中的战斗机

向量类型是R语言的核心,很多运算都涉及到向量。我们先来一个最简单的向量。

直接在控制台(Rstudio左下角的那个地方)中输入:

1
x <- 3x

x就是一个向量,一维的,3就是元素,只是只有一个,光杆司令。

1
2
3
4
5
6
y <- c(1, 2, 3)
y
length(y)
y[1]
y[3]
x[1]

“<-”是R语言中赋予的符号,在这个公式中,意味着我们把3这个值赋给x这个变量名,这样x就是3,但是3不一定是x。

我们再创建一个长度不是1的向量,

1
2
3
4
5
6
y <- c(1, 2, 3)
y
length(y)
y[1]
y[3]
x[1]

y就是一个长度为3的向量。length函数是用来查看向量长度的函数。对于y向量来说,他的元素就是1,2,3。向量中的元素是有顺序的,第一个从1开始计算,比如我想要查看向量y的第一个元素,那么就使用y[1]。y[1]代表的就是向量y的第一个元素,那么同样的y的第三个元素就是y[3]。

矩阵(matrix):

R中的矩阵概念跟数学中完全一样。矩阵可以理解为二维的向量,也就是包括行数和列数。下面我们举个栗子:

1
2
m <- rbind(c(1,2), c(2,3))
m

rbind是一个函数,我们简单介绍一下。rbind直接从字面上理解是row combine,也就是行组合。什么意思呢?就是把几个向量按照行从上到下组合起来。用图可以很好地理解:(灵魂画师上线。。。)

rbind函数就是把向量1和2然后按照行组合起来,成为一个新的矩阵。从这上面可以看到必须满足两个条件,向量1和2必须长度一致。同样还有一个函数,cbind,大家应该猜到了,就是column combine,也就是把几个向量按照列从左到右组合。
比如

1
2
m2 <- cbind(c(1,2,3,4), c(2,3,4,5))
dim(m2)

dim是查看矩阵的行数和列数。那么,怎么查看矩阵的某个元素呢?和向量一样,需要使用[]来看,但是要看单个元素,得使用行数和列数来进行确定,比如:

1
m2[1,2]

就是指m2的第1行和第2列的元素是什么东东。
那么,能够直接从矩阵中提取子向量和子矩阵呢?当然阔以:

1
2
3
4
5
6
7
8
m2[2,]
m2[,1]
m2[c(1,2),]
m2[2,]#就是指m2的第2行的所有元素

m2[,1]#就是指m2的第1列的所有元素m2[c(1,2),]

#就是指把m2的第1和第2行的矩阵

注意:#也就是井号,是注释的意思,也就是#后面的代码是不会运行的,所以如果你想要对代码进行注释说明,那么就可以像上面一样,加一个#号,然后在后面写下注释的内容。

数据框(data.frame)

数据框和矩阵是非常类似的,但是完全不一样的,就是矩阵中的元素类型必须一致,而数据框中的元素类型则可以不同。
比如我们有3个同学的语文,数学成绩。第一列是同学的名字(字符型,character),第二列是语文成绩(数字型,numeric),第三列是数学成绩(数字型,numeric)。这时候就得使用数据框了。

1
2
3
4
d <- data.frame(name = c("Shen", "Li", "Tu"),
Yuwen = c(90, 95, 100),
Shuxue = c(85,80,90))
d

使用函数data.frame建立数据框。name,Yuwen,Shuxue是建立好之后的数据框d的列名。
前面介绍了R语言中的最为常见的几种数据结构,差不多可以满足我们后面的需求了。后面再遇到新的东西,我们再介绍。

2. 函数

和其他语言一样,函数是R语言中的战斗机。和数学概念一样,y = f(x),也就是输入变量x,然后经过函数f的处理,得到结果y。
函数也是一个对象,比如:

1
2
3
fun <- function(a, b) {
a + b
}

fun就是一个最为简单的函数。他要做的就是把a和b加起来,输出a和b之和。在这里面,a和b是形参(形式参数),不是具体的参数。只是用来把位置占起来。那么我们使用fun来做一个简单的计算:

1
fun(1,2)

这时候1和2就是实参(实际参数),也就是把1和2带进去,计算1和2的和。很多函数(也就是功能)都是别人已经写好的,比如求和sum,求根sqrt等等。 函数写好之后,很多都在不同的包(package)中,下面就介绍一下R语言中的包。

3. 包(package)

R包,类似C、Python中库的概念,指包含特定领域的函数、数据、文档等的集合。通过调用包,可以直接使用包中现成的数据、函数等,使开发方便快捷高效。

  • R的强大在于包含了各种各样的包,使用包非常有利于便捷开发。
  • 一些功能在现有的包中并不存在,需要自己实现,实现后通过打包方便代码的复用。
  • 每个包涵括一个领域相关的函数数据文档等,通过包可以有效地组织代码结构,有利于开发.

在R中,所有的函数都是封装好,放在包中的,包是R语言互相交流的最好方法。下面简单介绍一下包的安装,使用,等我们后面会介绍怎么创建,分享自己的R包。
如何找到自己的需要的R包
大家要相信很多方法都已经是存在的了,所以你能想到的功能,很多已经被别人做过了,已经有了现成的工具可以使用,也就是R包。因此,如果想实现一个功能,首先要想到的就是是不是已经存在了这样的工具,而不是自己去写。那么如何找到自己所需要的R包呢,这时候就得靠搜索引擎了(谷歌一下,你就知道)。当然,在我大天朝,谷歌早已经404了,如果你能够使用,或者能够翻墙接触到这个万恶的网站,那是最好的,如果不行,那就用必应(必应内心OS:终于想我到了啊。。。),最次,使用百度(百度:嗯?)。

如何安装R包

R包一般会公开在三个地方:

  • CRAN:The Comprehensive R Archive Network,这是R core team的R包官方存放地方。安装放在这些的R的网址是https://mirrors.tongji.edu.cn/CRAN/,怎么安装这些包呢?使用下面代码即可:
1
install.packages("package.name")#将package.name替换为自己的要安装的包的名字即可。
  • Bioconductor:这个是大部分的生信相关的包存放的地方,网址https://www.bioconductor.org/,安装使用代码如下:
1
2
3
4
## try http:// if https:// URLs are not supported
source("https://bioconductor.org/biocLite.R")
biocLite("package.name")
#将package.name替换为自己的要安装的包的名字即可。
  • Github:CRAN和bioconductor都需要经过一定的审核,才可以发布,因此另外一个最为常用的代码托管网站,github,就更为常用友好了,只要把自己写的代码和包放在github上,就可以安装R包了。安装github上的R包,需要使用到R语言大神Hadley Wickham写的包,devtools了。因此,首选需要安装这个包,代码如下:
1
2
3
4
install.packages("devtools")
library(devtools)
install_github("user.name/package.name")
#其中user.name是github作者的用户名。

如何获得R包及函数的帮助文档

安装包之后,如要观察这个包里面所有的函数,以及这个R的说明,比如对于我们刚才安装的sxtTools这个包。我们就可以使用下面代码看这个包里面都含有哪些函数:

1
help(package = "sxtTools")

如果想看某个具体的函数的信息,可以使用下面的代码:

1
?installBioc

每个函数都有详细的说明,可以通过看这些详细的官方文档来学习怎么使用,如果想更便捷,那就直接上~google必应(万不得已,用百度。。。)。

下期预告
下期我们开始学习怎么读取,输出数据,敬请期待~~~