R函数了解一下？🤓

刘小泽写于18.8.29-30

R函数是什么？

其实就是对一些编程语言的封装，编写函数可以减少重复代码的书写，让R脚本可以简洁运行，并且可以快速进行传播。

一般一个函数可以完成一项特定功能，统计中的函数比如sum、var、mean等，都有自己的工作范围

查看函数源代码

没有被封装的：敲函数名不加括号就可以查看

自定义的函数

• 函数名称：与功能相关，必须字母开头
• 函数声明：function声明，就是告诉R里面的是整套函数
• 函数参数
• 函数体：最后需要用return 返回结果，也就是我们用这个函数的目的

myfunc <- function (选项参数) {
    函数体
}
# 不同选项之间用逗号分隔，选项参数之间用等号，可以只写选项，但如果有参数的话需要标明参数默认值
# 大括号中的函数体是最重要的部分，会包含许多逻辑判断、循环等

一定会用到的循环判断

常用的格式

• if条件判断，如score=70;ifelse(score>60, print("pass"), print("fail"))
• for循环，如for (i in 1:10) {print ("hello, world")}
• while循环，如i=1;while (i <= 10) {print("hello, world"); i=i+1}
• switch语句

循环语句包括？

• 条件判断T/F？
• 执行循环的结构
• 表达式

输入数据类型

简单句几个例子

• 向量：sum、mean、sd、range、median、sort、order等
• 矩阵或数据框：cbind、rbind等
• 数值矩阵：heatmap等

R初始加载的函数有1200多个，很难也没必要全部记住，只需要了解他们的规律，看他们的命令基本就能明白意思，具体的使用就需要知道函数的选项、参数使用。选项负责解决做还是不做的问题，参数解决做的话做多少的问题。有时候为什么会感觉：明明一件感觉很简单的工作，脑子里想的很明白，第一步第二步，然后实际用R却无从下手，最后默默用了excel…就是因为对选项参数不了解，也就是利器虽有，不在你手～

选项参数

与linux及其他的编程语言类似

• 输入控制

  一般放在第一位，告诉我们函数能接受哪种类型数据

  file:后接一个文件【比如read.table、read.csv等】；
  data：输入一个数据框；
  x：单独的一个对象，一般是向量，但也可以是矩阵/列表；
  x，y：需要两个输入变量；
  x，y，z：需要三个输入变量；
  formula：公式【如：lm等统计函数】
  na.rm：去掉缺失值
  ...：没有输入限制【如：data.frame】
  

• 输出控制

• 调节

  1. 根据名字判断选项作用:

  color:控制颜色
  font：字体
  select：与选择有关
  lty：这是一类简写【line type的简写】
  lwd：line width
  method：算法
  

  1. 选项接受哪些参数

  main:设置标题用的，因此要求字符串【不能是向量】
  na.rm：只能是T/F这样的逻辑值
  axis：其中的side控制坐标轴方向【方向只有四个，取值只能是1-4】
  fig：图形的位置【需要给出四个元素的向量，定位图片四个角】
  

数学统计相关函数

数据统计是数据分析的基础，包括了sum、mean、方差var、标准差sd、对数log、指数exp、绝对值abs等

如果再结合apply家族函数，就能对矩阵、数据框的行/列进行统计

概率分布

概率分布

这些概率分布前面加上d、p、q、r就构成了函数

其中，d表示密度，p表示分布（你会想：分布不是distribution吗？没错，但是密度density出现的比较早，d被抢占了）、q表示分位数（quantile）、r表示随机数（random）

比如像rnorm、runif、dexp等

单是正态分布的四种情况dorm、pnorm、rnorm、qnorm，除了都有的mean、sd设置以外，其余的参数设置还是不尽相同

生成随机数runif

随机生成0-1的数，runif(n) n是个数

如果想生成1以上的数，可以用runif(20)*10就生成了20个大于1的个位数；但是这样做是很麻烦的。因此再学习runif的帮助信息，发现有min、max的设置，因此直接runif(20,min=1,max=10)

随机数的特点就是“随机”，每次的数肯定不同，但有时为了重复验证，需要同样的“随机数”，因此，可以在runif之前用set.seed() 设定一下，括号中给一个任意数字，就可以绑定这一组随机数到我们这里的seed，以后再想看这次的随机数，就再次输入set.seed(之前的随机数)，再运行即可

描述性统计—获得数据的简历

• summary 运行就能对数据集得到详细统计【最值、四分位数、数值型变量均值（或因子/逻辑向量的频数统计）】

• fivenum 直接返回最值、上/下分位数、中位数

• Hmisc::describe 返回观测值、缺失值、唯一值的数目、均值、分位数、5个最大值、5个最小值

• pastecs::stat.desc 设置basic=TRUE，返回基本的统计（观测值、空值、最值、范围、总数）；设置desc=TRUE，返回描述性统计（中位数、平均值、标准误SE、置信区间CI、p值、方差var、标准偏差dev、方差系数coef.var）；设置norm=T，返回标准分布统计（skew、kurtosis、normtest）

• psych::describe 设置trim=0.1就能去除最低/最高10%的部分

• 【对数据分组描述】MASS::aggregate ，将分组信息通过一个列表指定出来

  >aggregate(mtcars[c("mpg","disp","hp")],by=list(Cylinder=mtcars$cyl),mean)
  #这里就是根据mtcars中的cyl气缸数，对mpg、disp、hp进行统计，并计算他们的平均数
    Cylinder      mpg     disp        hp
  1        4 26.66364 105.1364  82.63636
  2        6 19.74286 183.3143 122.28571
  3        8 15.10000 353.1000 209.21429
  #结果发现8缸的车比4缸马力（hp）大了2.5倍左右，但是相同油量它行驶的路程只有4缸的0.57倍
  

• doBy::summaryBy

  > summaryBy(mpg+disp+hp ~ cyl, data=mtcars, FUN=mean)
    cyl mpg.mean disp.mean   hp.mean
  1   4 26.66364  105.1364  82.63636
  2   6 19.74286  183.3143 122.28571
  3   8 15.10000  353.1000 209.21429
  

  但是，这个函数一次只能返回一个统计值，如果想一次统计既分了组，又做多个统计，那么可以用psych::describe.by

• psych::describe.by

  > describeBy(mtcars,list(cylinder=mtcars$cyl))
  # 优点是：可以进行分组；缺点就是：不能使用自定义的函数
  

如果两个包中的函数同时被载入，比如这里的describe，那么后导入的包函数会覆盖前面的

频数统计～哦应该是统计出现次数

比如想比较不同组之间的差异，就需要利用频数；但分组怎么来的呢？借助因子！

如何分组？

如果一个数据本身是因子的话，就可以直接分组，例如mtcars中的cyl气缸数统计值

##如果原数据就是因子型
> mtcars$cyl = as.factor(mtcars$cyl) #将cyl【转为因子】
> split(mtcars,mtcars$cyl) #split函数根据cyl因子【进行分组】

##如果原数据不是因子型
#比如mtcars中的mpg，【利用cut函数】（他的介绍就是Convert Numeric to Factor）
> cut(mtcars$mpg,c(seq(10,40,10)))#指定切割区间，10-40范围中，每隔10切一次，然后把mpg的原始数据归类进去
#将刚才cut的结果，【结合table进行频数统计】，就实现了分组
> tmp5 = cut(mtcars$mpg,c(seq(10,40,10))) %>% table()
.
(10,20] (20,30] (30,40] 
     18      10       4 

分组后进行统计？

#有了table计算出来的频数，如果要计算频率分布，可以用prop.table
> prop.table(tmp5)                
.
(10,20] (20,30] (30,40] 
  56.25   31.25   12.50 
#计算百分比，只需要在后面乘以100即可
#######################以上是一维的表格统计#####################
#先生成二维数据框，可以使用table或者xtabs
> library(vcd)                 
> table(Arthritis$Treatment,Arthritis$Improved)
         
          None Some Marked
  Placebo   29    7      7
  Treated   13    7     21
#或者也可以
> with(Arthritis,{table(Treatment,Improved)})
         Improved
Treatment None Some Marked
  Placebo   29    7      7
  Treated   13    7     21  
#最方便应该是xtabs
> tmp6 = xtabs(~Treatment+Improved, Arthritis)
         Improved
Treatment None Some Marked
  Placebo   29    7      7
  Treated   13    7     21

#再对列联表进行统计频数【margin.table】
> margin.table(tmp6,1) #其中1指行，表示对行进行频数统计；2指列
Treatment
Placebo Treated 
     43      41  
#统计频率【prop.table】
> prop.table(tmp6,1)
         Improved
Treatment      None      Some    Marked
  Placebo 0.6744186 0.1627907 0.1627907
  Treated 0.3170732 0.1707317 0.5121951
#将计算的边际频数添加到二联表中【addmargins】可以只添加行或者列
> addmargins(tmp6,2)
         Improved
Treatment None Some Marked Sum
  Placebo   29    7      7  43
  Treated   13    7     21  41  
#要做三维的列联表，只需要使用【xtabs+ftable】
> xtabs(~Treatment+Improved+Sex, Arthritis) %>% ftable()
                   Sex Female Male
Treatment Improved                
Placebo   None             19   10
          Some              7    0
          Marked            6    1
Treated   None              6    7
          Some              5    2
          Marked           16    5                 

独立性检验是什么？🤒

前提要有频数统计信息，根据它判断两类因子彼此之间有多大的缘分（相关性）；独立就是说二者无缘（唉！）比如一种气味只吸引男性，那么对这种气味来讲，男女之间就是独立的

先要知道假设检验(Hypothesis testing)

它根据一定的假设条件，由样本估计总体。设定假设：
1.原假设（两个因子无关=》没有发生臆想的事）；2.被择假设（两个因子有关=〉真的发生了预测的事😮）。这两种假设在逻辑上互补

https://www.douban.com/group/topic/116857671/

p-value：（Probability）当原假设为真时，得到的最大的检验概率。
一般阈值⍺定为5%。【不是绝对，如果数据噪音比较大，可以稍微放宽限制，设为10%；如果想再精细化数据进行检测，就可以设为1%】
当得到的p值小于这个⍺值时，拒绝原假设。一句话，p值越小，越能拒绝原假设

卡方检验chiseq.test

根本思想就是在于比较理论频数和实际频数的吻合程度

## tmp6是上面得到的improved与treatment的频数关系表
> chisq.test(tmp6)

    Pearson's Chi-squared test

data:  tmp6
X-squared = 13.055, df = 2, p-value = 0.001463
#p值明显小于5%，因此treatment与improved是有关的

## 再看一下improved与sex是否有联系
> tmp7 = xtabs(~Improved+Sex, Arthritis)
> chisq.test(tmp7)

    Pearson's Chi-squared test

data:  tmp7
X-squared = 4.8407, df = 2, p-value = 0.08889
#p值大于5%，因此二者没有关系

Fisher精确检验

它比较适用于小样本，仅能用于二维列联表中，行和列相互独立。计算每一种可能的概率，然后合计假设检验情况下的概率

Cochran-Mantel-Haenszel检验

原理：需要三个变量，假设：前两个变量分别与第三个变量相比，都是独立的
【注意】三个变量的顺序对于结果的影响很重要

#建立一个Improved和Treatment在Sex基础上的列联表
> tmp8 = xtabs(~Improved+Treatment+Sex, Arthritis)
> mantelhaen.test(tmp8)

    Cochran-Mantel-Haenszel test

data:  tmp8
Cochran-Mantel-Haenszel M^2 = 14.632, df = 2, p-value = 0.0006647 
#p值小于5%，说明Improved+Treatment与Sex相比都是独立的

#再建立一个Treatment和Sex在Improved基础上的列联表
> xtabs(~Treatment+Sex+Improved,Arthritis) %>% mantelhaen.test()

    Mantel-Haenszel chi-squared test with continuity correction

data:  .
Mantel-Haenszel X-squared = 2.0863, df = 1, p-value = 0.1486
#结果表示Treatment和Sex不全是独立的

相关性分析-测测缘分几何

衡量指标

Pearson相关系数、Spearman、Kendall、偏相关、多分格（polychoric）、多系列（polyserial）相关系数

相关函数1 cor()

可以指定method=Person/Spearman/Kendall。默认Person

结果产生一个对角矩阵，相关系数0-1之间，正负号表示正相关还是负相关

相关函数2 cov()

计算协方差，即两个变量的总体误差。计算偏相关时需要用到协方差结果

它反映的结果和cor（）一致

如果直接使用cor（）或cov（），得到的结果是所有样本的相关性，如果想要检测某几个样本的相关性，可以这样：

> a1=state.x77[,c(2,3,6)]
> b1=state.x77[,c(5,7)]
> cor(a1,b1)
               Murder      Frost
Income     -0.2300776  0.2262822
Illiteracy  0.7029752 -0.6719470
HS Grad    -0.4879710  0.3667797

除了Person/Spearman/Kendall这三个，还想计算其他的怎么办？【用扩展包】

比如计算偏相关：用pcor(u, S)

第一个参数u：一组数字组成的向量（前两个数字是必需的计算相关的列号，其余数字是作为辅助、参考的列号）；

第二个参数S：cor（）计算出的协方差结果

# 先看下state.x77的列名
> colnames(state.x77)
[1] "Population" "Income"     "Illiteracy" "Life Exp"   "Murder"     "HS Grad"    "Frost"     
[8] "Area" 
#想着重看一下Illiteracy和Murder之间的相关性，并且采用Population、Income、HS Grad作为参考
> pcor(c(3,5,1,2,6),cov(state.x77))
[1] 0.5989741 #相关性还比较高

相关性检验-分析完了，不得检验看看？

比如上面👆计算出来的相关性有的在0.6/0.7以上，我们感觉比较高，那么到底是真高还是我们的幻觉？这里就需要一步验证

其实我们上面独立性检验也是需要分析（算出p值）+验证（与⍺比较），只不过过程简单就一步到位

检验某一特定组的相关性

使用cor.test()，需要四个较为重要的选项参数

• x、y是需要检测相关性的变量（要求长度一致）
• alternative：可选"two.sided", "greater" or "less"，双侧检验还是单侧
  分别表示同时检测正负相关、正相关、负相关
• method：与cor函数一样，Person/Spearman/Kendall三个选项【可缩写】

> cor.test(state.x77[,3],state.x77[,5],alternative = "two.sided",method="p")

    Pearson's product-moment correlation

data:  state.x77[, 3] and state.x77[, 5]
t = 6.8479, df = 48, p-value = 1.258e-08 #p远小于5%，可以拒绝原假设
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval: #置信区间
 0.5279280 0.8207295
sample estimates:
      cor 
0.7029752  #也给出了相关系数

cor.test()除了计算p、cor以外，还给出了置信区间（95 percent confidence interval）

置信区间：算出一个计算的概率发生的范围，判断我们算的这个概率是否真实有效

比如：2018年9月新iPhone的最大6.5寸叫什么名字？综合考虑各种曝光率推测，叫XL的可能性为85%，而置信水平0.95以上的置信区间算出来是80%-90%，那么这个名字“XL”的真实度有95%以上的几率落在80%-90%之间。之前得到的85%落在这个范围内，因此之前推测的85%可能性为假的概率只有不到5%

同时检测多组相关性

cor.test（）一次只能检验一组变量，如果有多组就比较麻烦，可以使用psych::corr.test【不难发现，它多了一个r，这个就是递归recursive的意思】

使用很简单，直接corr.test(矩阵) 即可，不仅可以计算相关系数，还给出了检测值【校正过的】

偏相关检验

可以用ggm::pcor.test()

#先利用pcor计算偏相关系数
p1 = pcor(c(3,5,1,2,6),cov(state.x77))
#需要三个选项参数: 
# r：pcor计算的偏相关系数
# q：需要控制的变量数（就是除了必需的那两个以外还剩几个）这里是1，2，6这三个
# n：样本数 【用nrow统计一下行数】这里是50

> pcor.test(p1,3,50)
$tval #t检验
[1] 5.01773
$df #自由度
[1] 45
$pvalue 
[1] 8.671525e-06

分组数据相关性检验

经常会用到对两个组进行的检验，比如处理组与对照组

可以使用t检验也叫做student (Student's t test)【为什么叫学生，这里有一个小故事可以去查一下】，是用t分布理论来推论差异发生的概率，从而比较两个平均数的差异是否显著。主要用于样本含量较小（例如n<30），总体标准差σ未知的正态分布

两个组之间比较 t.test

#格式是t.test(y~x)
#y是数字型变量，x是因子型变量
> t.test(Prob ~ So, UScrime)

    Welch Two Sample t-test

data:  Prob by So
t = -3.8954, df = 24.925, p-value = 0.0006506 #p值很小，可以推翻原假设（So代表的南方与非南方地区的Prop值有相关性）
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.03852569 -0.01187439
sample estimates:
mean in group 0 mean in group 1 
     0.03851265      0.06371269 

多个组之间比较

• 数据满足正态分布，可以用方差分析

• 不满足正态分布，可以用非参数检验

  非参数检验：Nonparametric tests，当总体方差未知或知道较少的情况下，用样本数据对总体的分布进行推断

  参数检验：总体分布已知（比如满足正态分布），对总体分布的参数如均值、方差等进行推断

数据类型	非参数检验	参数检验
单样本	Wilcoxon符号秩和检验 Mann–Whitney U 检验	单样本t检验
配对样本	Wilcoxon配对符号秩和检验	配对t检验
两独立样本	Wilcoxon秩和检验	两独立样本t检验
多组独立样本	Kruskal-Wallis H 检验	单因素方差分析
随机区组	Friedman检验	两因素无交互作用方差分析

做个图？

图片最大的价值在于促使我们发现从未预料到的事情 - John Tukey

四个作图系统

• 基础绘图：ls("package:graphics") 显示内置的作图函数，demo("graphics") 展示
• lattice包、
• ggplot2包、
• grid包

难点--输入什么样的数据

散点图要求x、y坐标；

直方图（图形不是连贯的，是分组的）需要因子型变量；

热图的各个色块组成就像一个矩阵，每个色块就是一个数值，颜色深的色块代表数值越大，因此它需要矩阵

箱线图需要一个因子型和一个数值型

不同的数据对应不同的图

【简单了解下就好】plot函数如何识别不同类型的数据？

好像同一组数据，但是不同的处理方法都能出图，那么函数到底识别哪种数据格式呢？

因为它属于S3系统【面向对象】

• 属性：R中每个对象（变量、函数等）可以添加很多属性
• 泛型函数
• 方法

因此，看似plot是万能的，但其实它的背后是一大堆成型的智囊团（函数），遇到什么数据就调取哪个函数，类似的还有summary、print

plot做宏观布局，par函数在plot基础上进行微调，par有72变,输入names(par())可以查看全部；
当然现在作图基本都用ggplot2，但是plot的思想还是要学习的

---

欢迎关注我们的公众号～_～　　
我们是两个农转生信的小硕，打造生信星球，想让它成为一个不拽术语、通俗易懂的生信知识平台。需要帮助或提出意见请后台留言或发送邮件到Bioplanet520@outlook.com

Welcome to our bioinfoplanet!