刘小泽写于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】 输出控制
-
调节
- 根据名字判断选项作用:
color:控制颜色 font：字体 select：与选择有关 lty：这是一类简写【line type的简写】 lwd：line width method：算法- 选项接受哪些参数
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
