R语言学习笔记（16）-R函数(2)
2021-02-04 本文已影响0人
Akuooo
六、独立性检验函数
1.独立性检验
独立性检验是根据频数信息判断两类因子彼此相关或相互独立的假设检验。所谓独立性就是指变量之间是独立的，没有关系。
三种方法：卡方检验、Fisher检验、Cochran-Mantel-Haenszel检验
- 假设检验（Hypothesis Testing）
假设检验是数理统计学中根据一定假设条件由样本推理总体的一种方法。
原假设——没有发生；
备择假设——发生了；
具体做法：根据问题的需要对所研究的总体作某种假设，记作H0；选取合适的统计量，这个统计量的选取要使得在假设H0成立时，其分布为已知；由实测的样本，计算出统计量的值，并根据预先给定的显著性水平进行检验，作为拒绝或接受假设H0的判断。
- p-value
就是probability的值，是通过计算得到的概率值，在原假设为真时，得到最大的或者超出所得到的额检验统计量值的概率
一般p值定位到0.05，当p<0.05拒绝原假设；p>0.05，不拒绝原假设
记忆：p值越小，假设越不靠谱，拒绝原假设；p值越大，假设越靠谱 - 卡方检验
#依旧使用Arthritis的Treatment和Improved来判定是否相互独立
#若独立则说明，改善与药物治疗无关，即药物无效果
#table()函数计算频数
> mytable <- table(Arthritis$Treatment,Arthritis$Improved)
#运用chisq.test检验
> chisq.test(mytable)
Pearson's Chi-squared test
data: mytable
X-squared = 13.055, df = 2, p-value = 0.001463
# p值小于0.05，说明,这两个变量不是独立的，而是相关联的。
#再来看Improved和sex是否有关
> mytable <- table(Arthritis$Sex,Arthritis$Improved)
> chisq.test(mytable)
Pearson's Chi-squared test
data: mytable
X-squared = 4.8407, df = 2, p-value = 0.08889
Warning message:
In chisq.test(mytable) : Chi-squared近似算法有可能不准
#可以看出，p值0.08大于0.05，说明两者无关，两变量独立，即药物无性别差异
#这些待检测变量无顺序差别
- Fisher精确检验
#fisher.test()函数进行检验
#原理略有不同
> mytable <- xtabs(~Treatment+Improved, data=Arthritis)
> fisher.test(mytable)
Fisher's Exact Test for Count Data
data: mytable
p-value = 0.001393
alternative hypothesis: two.sided
#p值为0.0013，与卡方检验的结果略有不同
- Cochran-Mantel-Haenszel检验
#使用mantelhaen.test()函数
#假设原理：两个名义变量在第三个变量每一层中都是条件独立的
#故需要三个变量
#该顺序验证的是，在性别层面上，treatment和Improved是否相互独立
> mytable <- xtabs(~Treatment+Improved+Sex, data=Arthritis)
> mantelhaen.test(mytable)
Cochran-Mantel-Haenszel test
data: mytable
Cochran-Mantel-Haenszel M^2 = 14.632, df = 2, p-value = 0.0006647
#p值很小，说明在性别层面上，药物质量和改善情况在每一个水平上不独立
#若改变变量顺序，则结果是有差别的
> mytable <- xtabs(~Treatment+Sex+Improved, data=Arthritis)
> mantelhaen.test(mytable)
Mantel-Haenszel chi-squared test with continuity correction
data: mytable
Mantel-Haenszel X-squared = 2.0863, df = 1, p-value = 0.1486
alternative hypothesis: true common odds ratio is not equal to 1
95 percent confidence interval:
0.8566711 8.0070521
sample estimates:
common odds ratio
2.619048
七、相关性分析函数
- 相关性分析
相关性分析是指对两个或多个具备相关性的变量元素进行分析，从而衡量两个变量因素的相关密切程度。相关性的元素之间需要存在一定的联系或者概率才可以进行相关性分析，也就是变量之间是否有关系。
-
相关系数
首先是验证两个变量是否相互独立，若不独立，就存在相关性。相关可以是正相关或负相关，需要根据相关系数来确定。
相关系数可以用来描述变量之间的关系，用正负号表示是正相关还是负相关，相关系数的大小表示相互关系的强弱。 -
相关性衡量指标
Pearson相关系数、Spearman相关系数、Kendall相关系数、偏相关系数、多分格（polychoric）相关系数和多系列相关（polyserial）系数 -
如何计算相关系数
（1）cor()函数
计算相关性系数都是用同一种函数：cor()，可以计算三种相关系数（Pearson、Spearman&Kendall）
①参数：
method：指定使用哪种算法，默认是Pearson相关系数
use：如何对待缺失值，不处理or删除。可选，everything、all.obs、complete.obs、na.or.complete、pairwise.complete.obs
②演示（state.x77矩阵数据集）
#推测谋杀率与哪些因素有关
> cor(state.x77)
Population Income Illiteracy Life Exp Murder HS Grad Frost Area
Population 1.00000000 0.2082276 0.10762237 -0.06805195 0.3436428 -0.09848975 -0.3321525 0.02254384
Income 0.20822756 1.0000000 -0.43707519 0.34025534 -0.2300776 0.61993232 0.2262822 0.36331544
Illiteracy 0.10762237 -0.4370752 1.00000000 -0.58847793 0.7029752 -0.65718861 -0.6719470 0.07726113
Life Exp -0.06805195 0.3402553 -0.58847793 1.00000000 -0.7808458 0.58221620 0.2620680 -0.10733194
Murder 0.34364275 -0.2300776 0.70297520 -0.78084575 1.0000000 -0.48797102 -0.5388834 0.22839021
HS Grad -0.09848975 0.6199323 -0.65718861 0.58221620 -0.4879710 1.00000000 0.3667797 0.33354187
Frost -0.33215245 0.2262822 -0.67194697 0.26206801 -0.5388834 0.36677970 1.0000000 0.05922910
Area 0.02254384 0.3633154 0.07726113 -0.10733194 0.2283902 0.33354187 0.0592291 1.00000000
#解读数据，相关系数0~1之间
#文盲率与谋杀率有很强的正相关关系；收入高、高中毕业率与谋杀率呈负相关；谋杀率与地区的相关性很低
（2） cov()
计算协方差，衡量两个变量的总体误差，计算正相关时需要用到协方差的结果。
> cov(state.x77)
Population Income Illiteracy Life Exp Murder HS Grad Frost Area
Population 19931683.7588 571229.7796 292.8679592 -4.078425e+02 5663.523714 -3551.509551 -77081.97265 8.587917e+06
Income 571229.7796 377573.3061 -163.7020408 2.806632e+02 -521.894286 3076.768980 7227.60408 1.904901e+07
Illiteracy 292.8680 -163.7020 0.3715306 -4.815122e-01 1.581776 -3.235469 -21.29000 4.018337e+03
Life Exp -407.8425 280.6632 -0.4815122 1.802020e+00 -3.869480 6.312685 18.28678 -1.229410e+04
Murder 5663.5237 -521.8943 1.5817755 -3.869480e+00 13.627465 -14.549616 -103.40600 7.194043e+04
HS Grad -3551.5096 3076.7690 -3.2354694 6.312685e+00 -14.549616 65.237894 153.99216 2.298732e+05
Frost -77081.9727 7227.6041 -21.2900000 1.828678e+01 -103.406000 153.992163 2702.00857 2.627039e+05
Area 8587916.9494 19049013.7510 4018.3371429 -1.229410e+04 71940.429959 229873.192816 262703.89306 7.280748e+09
#反应出来的问题与cor()一致
#若只想比较某几个变量
> x <- state.x77[,c(1,2,3,6)]
> y <- state.x77[,c(4,5)]
> head(x)
Population Income Illiteracy HS Grad
Alabama 3615 3624 2.1 41.3
Alaska 365 6315 1.5 66.7
Arizona 2212 4530 1.8 58.1
Arkansas 2110 3378 1.9 39.9
California 21198 5114 1.1 62.6
Colorado 2541 4884 0.7 63.9
> head(y)
Life Exp Murder
Alabama 69.05 15.1
Alaska 69.31 11.3
Arizona 70.55 7.8
Arkansas 70.66 10.1
California 71.71 10.3
Colorado 72.06 6.8
> cor(x,y)
Life Exp Murder
Population -0.06805195 0.3436428
Income 0.34025534 -0.2300776
Illiteracy -0.58847793 0.7029752
HS Grad 0.58221620 -0.4879710
（3）计算其他相关系数，可通过R扩展包实现
ggm包中的pcor()可计算偏相关系数
pcor()需要输入两个重要参数
u：数值向量，前两个数值表示要计算相关系数的下标，其余数值为条件变量的下标
s：cov计算出来的协方差结果
> colnames(state.x77)
[1] "Population" "Income" "Illiteracy" "Life Exp" "Murder" "HS Grad" "Frost"
[8] "Area"
#若想控制收入水平、文盲率和高中毕业率的影响，看人口和谋杀率之间的关系
#人口是第一列，谋杀是第五列
> pcor(c(1,5,2,3,6),cov(state.x77))
[1] 0.3462724
八、相关性检验函数
进行相关性系数计算后，如何进行统计学显著性检验？
计算出来的相关性系数到底能否说明问题，还需要进行量化，需要p-value衡量
- cor.test()函数 相关性系数的检验
（1）参数
x，y：需要检测的相关性变量
alternative：指定是进行双侧检验还是单侧检验，取值为“two.sided”（同时检测政府相关性）、“greater”（单独正相关）、“less”（负相关）
method：Pearson、Spearman、Kendall
（2）演示
> cor.test(state.x77[,3],state.x77[,5])
Pearson's product-moment correlation
data: state.x77[, 3] and state.x77[, 5]
t = 6.8479, df = 48, p-value = 1.258e-08
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.5279280 0.8207295
sample estimates:
cor
0.7029752
#p值远小于0.05，拒绝原假设，说明两个变量之间相关系数不为零
#confidence interval：置信区间
#cor相关系数
（3）置信区间
置信区间：指由样本统计量所构造的总体参数的估计区间。在统计学中，一个概率样本的置信区间是对这个样本的某个总体参数的区间估计。置信区间展现的是这个参数的真实值有一定概率落在测量结果的周围的程度。置信区间给出的是被测量参数的测量值的可信程度。
比如，在一次大选中某人支持率为55%，而置信水平0.95以上的置信区间是（50%，60%），那么他的真实支持率有95%的几率落在百分之五十和百分之六十之间，因此他的真实支持率不足一半的可能性小于5%。
简单来说就是，光给出概率还不行，还要给出这个概率的范围
- psych包中的corr.test()函数
可一次进行多个变量之间的计算
> corr.test(state.x77)
Call:corr.test(x = state.x77)
Correlation matrix
Population Income Illiteracy Life Exp Murder HS Grad Frost Area
Population 1.00 0.21 0.11 -0.07 0.34 -0.10 -0.33 0.02
Income 0.21 1.00 -0.44 0.34 -0.23 0.62 0.23 0.36
Illiteracy 0.11 -0.44 1.00 -0.59 0.70 -0.66 -0.67 0.08
Life Exp -0.07 0.34 -0.59 1.00 -0.78 0.58 0.26 -0.11
Murder 0.34 -0.23 0.70 -0.78 1.00 -0.49 -0.54 0.23
HS Grad -0.10 0.62 -0.66 0.58 -0.49 1.00 0.37 0.33
Frost -0.33 0.23 -0.67 0.26 -0.54 0.37 1.00 0.06
Area 0.02 0.36 0.08 -0.11 0.23 0.33 0.06 1.00
Sample Size
[1] 50
Probability values (Entries above the diagonal are adjusted for multiple tests.)
Population Income Illiteracy Life Exp Murder HS Grad Frost Area
Population 0.00 1.00 1.00 1.00 0.23 1.00 0.25 1.00
Income 0.15 0.00 0.03 0.23 1.00 0.00 1.00 0.16
Illiteracy 0.46 0.00 0.00 0.00 0.00 0.00 0.00 1.00
Life Exp 0.64 0.02 0.00 0.00 0.00 0.00 0.79 1.00
Murder 0.01 0.11 0.00 0.00 0.00 0.01 0.00 1.00
HS Grad 0.50 0.00 0.00 0.00 0.00 0.00 0.16 0.25
Frost 0.02 0.11 0.00 0.07 0.00 0.01 0.00 1.00
Area 0.88 0.01 0.59 0.46 0.11 0.02 0.68 0.00
To see confidence intervals of the correlations, print with the short=FALSE option
- ggm包中的pcop.test()函数
可进行偏相关的检验
#需要三个选项
#首先是pcor要检验的偏相关系数
#然后是要控制的变量数
#最后是样本数
> x <- pcor(c(1,5,2,3,6),cov(state.x77))
> pcor.test(x,3,50)
$tval##学生t检验
[1] 2.476049
$df##自由度
[1] 45
$pvalue##p值
[1] 0.01711252
4.分组数据的相关性检验
（1）两组：这种检验可以使用t检验。
t为student学生的意思；t检验是运用t分布理论推论差异发生的概率，从而比较两个平均数的差异是否显著。
主要用于样本含量较小（＜30）的总体标准差未知的正态分布数据
在这里可使用MASS包中的UScrime数据集，包含了1960年美国47个洲的刑罚制度对犯罪率影响的信息。
#"~"分组变量，将南北分成两组做t检验
> t,test(Prob ~ So, data = UScrime)
Welch Two Sample t-test
data: Prob by So
t = -3.8954, df = 24.925, p-value = 0.0006506
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
#p值小于0.05，可以拒绝南北各洲拥有相同监禁概率的假设
（2）若想在多于两个组之间进行比较
如果数据满足正态分布，可以使用方差分析；
如果不满足正太分布，就需要使用非参数检验的方法
非参数检验，Nonparametric tests，在总体方差未知或知道甚少的情况下，利用样本数据对总体分布形态等进行推断的方法。由于非参数检验方法在推断过程中不涉及有关总体分布的参数，因而得名为“非参数”检验
参数检验，parametric tests，是在总体分布形式已知的情况下，对总体分布的参数如均值、方差进行推断的方法，也就是数据分布已知，比如满足正态分布。
但在数据分析过程中，由于种种原因，往往无法对数据总体的分布情况进行简单的假定，这时参数检验方法就无法适用了。
就需要使用非参数检验方法，它正是基于这种考虑，在总体方差未知或知道很少的情况下，利用样本数据对总体分布形态进行推断的方法。非参数检验放在在推断过程中不涉及有关总体分布的参数，因而被称为非参数检验
九、绘图函数
R语言四大作图系统
-
基础绘图系统；
基础绘图函数都在graphics包中，R启动时已经默认加载了这个包
graphics
#用此命令展示R能够作的一些图
#敲回车继续
> demo(graphics)
image.png
image.png
image.png
> help(package = "graphics")
graphics
(1)高级绘图：一步到位，可直接绘制出图（搭好框架）
(2)低级绘图：不能单独使用，必须在高级绘图产生图形的基础上，对图形进行调整，比如加一条线，加上标题文字等。（精雕细琢）
（3）输入数据格式
散点图：x和y两个坐标数据
直方图：因子
热力图：数据矩阵
最好能够看到哪个图形，就能猜出是用哪个函数做出来的
（4）plot()函数
①输入数据为向量
以women数据集为例
#plot函数可以接受单独的一个数值向量
> plot(women$height)
#也可以接受两个数值向量
> plot(women$height,women$weight)
height.png
height+weight.png
②输入数据为因子→直方图
以mtcars数据集为例
> plot(as.factor(mtcars$cyl))
as.factor(mtcars$cyl)
③因子+数值→箱线图
> plot(as.factor(mtcars$cyl),mtcars$carb)
> plot(mtcars$carb,as.factor(mtcars$cyl))
因子,数值箱线图.png
反过来
④两个都是因子→脊柱图
> plot(as.factor(mtcars$cyl),as.factor(mtcars$carb))
两个因子.png
⑤接受公式
如包含“~”结构，输出二者关系图
> plot(women$height ~ women$weight)
image.png
> fit <- lm(height ~ weight, data = women)
> fit
Call:
lm(formula = height ~ weight, data = women)
Coefficients:
(Intercept) weight
25.7235 0.2872
> plot(fit)
#会生成四幅图
S3系统：是R中面向对象的编程概念，R中每个对象可以添加很多属性，S3包括属性、泛型函数和方法三个内容。
plot之所以支持多种数据结构格式，原因在于plot支持多种属性的数据格式
> methods(plots)
methods(plot)
[1] plot.acf* plot.data.frame* plot.decomposed.ts* plot.default
[5] plot.dendrogram* plot.density* plot.ecdf plot.factor*
[9] plot.formula* plot.function plot.hclust* plot.histogram*
[13] plot.HoltWinters* plot.isoreg* plot.lm* plot.medpolish*
[17] plot.mlm* plot.ppr* plot.prcomp* plot.princomp*
[21] plot.profile.nls* plot.raster* plot.spec* plot.stepfun
[25] plot.stl* plot.table* plot.ts plot.tskernel*
[29] plot.TukeyHSD*
see '?methods' for accessing help and source code
#里面有很多plot开头的函数相当于plot函数是一个家族，包含很多个plot开头的函数。
#在进行plot会制作图时，就会支持很多种数据格式
（5）par()，par是parameter的简称，用来对绘图参数进行设置，有很多选项参数
> par()
$xlog
[1] FALSE
$ylog
[1] FALSE
$adj
[1] 0.5
$ann
[1] TRUE
$ask
[1] FALSE
$bg
[1] "white"
$bty
[1] "o"
$cex
[1] 1
$cex.axis
[1] 1
$cex.lab
[1] 1
$cex.main
[1] 1.2
$cex.sub
[1] 1
$cin
[1] 0.15 0.20
$col
[1] "black"
$col.axis
[1] "black"
$col.lab
[1] "black"
$col.main
[1] "black"
$col.sub
[1] "black"
$cra
[1] 14.4 19.2
$crt
[1] 0
$csi
[1] 0.2
$cxy
[1] 1.845331 0.860517
$din
[1] 5.541667 5.354167
$err
[1] 0
$family
[1] ""
$fg
[1] "black"
$fig
[1] 0 1 0 1
$fin
[1] 5.541667 5.354167
$font
[1] 1
$font.axis
[1] 1
$font.lab
[1] 1
$font.main
[1] 2
$font.sub
[1] 1
$lab
[1] 5 5 7
$las
[1] 0
$lend
[1] "round"
$lheight
[1] 1
$ljoin
[1] "round"
$lmitre
[1] 10
$lty
[1] "solid"
$lwd
[1] 1
$mai
[1] 1.02 0.82 0.82 0.42
$mar
[1] 5.1 4.1 4.1 2.1
$mex
[1] 1
$mfcol
[1] 1 1
$mfg
[1] 1 1 1 1
$mfrow
[1] 1 1
$mgp
[1] 3 1 0
$mkh
[1] 0.001
$new
[1] FALSE
$oma
[1] 0 0 0 0
$omd
[1] 0 1 0 1
$omi
[1] 0 0 0 0
$page
[1] TRUE
$pch
[1] 1
$pin
[1] 4.301667 3.514167
$plt
[1] 0.1479699 0.9242105 0.1905058 0.8468482
$ps
[1] 12
$pty
[1] "m"
$smo
[1] 1
$srt
[1] 0
$tck
[1] NA
$tcl
[1] -0.5
$usr
[1] 113.04 165.96 57.44 72.56
$xaxp
[1] 120 160 4
$xaxs
[1] "r"
$xaxt
[1] "s"
$xpd
[1] FALSE
$yaxp
[1] 58 72 7
$yaxs
[1] "r"
$yaxt
[1] "s"
$ylbias
[1] 0.2
#这些都是绘图的默认设置
par()中的参数可以单独设置，对整个绘图全局起作用；也可以用在每个函数中。
例如，在绘制mtcars因子时，给条形图加颜色，可以设置cyl选项
> plot(as.factor(mtcars$cyl),col=c("red","green","blue"))
image.png
- lattice包；
- ggplot2包；
- grid包
十、自定义函数
编写函数的目的就是为了减少重复代码的书写，从而让R脚本更为简洁高效，增加可读性
在R中，不敲函数的括号就可以显示其代码；但有些函数封装起来了，这种方法就看不了
- R自定义函数
（1）函数名称
函数命令与功能相关；
可以是字母和数字的组合，但必须是字母开头
（2）函数声明
利用function函数来声明
myfun <- function(选项参数)
{
函数体#最重要部分，包括很多逻辑判断和循环
}
（3）函数参数
（4）函数体
- 编写一个简单的函数
（1）计算偏度与峰度函数
相关概念:
偏度（skewness），是统计数据分布偏斜方向和程度的度量，是统计数据分布对称程度的数字特征
峰度（peakedness；kurtosis）又称峰态系数。表征概率密度分布曲线在平均值处峰值高低的特征数。
编写函数：
> 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))
+ }
> x <- 1:100
> mystats(x)
n mean stdev skew kurtosis
100.00000 50.50000 29.01149 45.22571 38.31910
- 循环与向量化操作
很多情况下都需要使用循环语句，尤其是在自定义函数中。
R向量化操作，在函数内部是通过循环来实现的，也就是说本来需要使用循环的地方，在函数内部已经帮你实现好了。如果是自定义函数，还是需要使用循环控制结构的。
R中有一般现代化编程语言中都有的控制结构，控制循环语句与其他语言也差不多（if条件判断，for循环，while循环，switch语句等）
（1）完整循环
①条件判断，是真或者假
②用于循环执行的结构
③表达式
（2）R中的for循环
> for (i in 1:10) {print ("Hello,World")}
[1] "Hello,World"
[1] "Hello,World"
[1] "Hello,World"
[1] "Hello,World"
[1] "Hello,World"
[1] "Hello,World"
[1] "Hello,World"
[1] "Hello,World"
[1] "Hello,World"
[1] "Hello,World"
（2）while循环
> i=1;while(i <= 10) {print ("Hello,World");i=i+1;}
[1] "Hello,World"
[1] "Hello,World"
[1] "Hello,World"
[1] "Hello,World"
[1] "Hello,World"
[1] "Hello,World"
[1] "Hello,World"
[1] "Hello,World"
[1] "Hello,World"
[1] "Hello,World"
> i=1;while(i <= 10) {print ("Hello,World");i=i+2;}
[1] "Hello,World"
[1] "Hello,World"
[1] "Hello,World"
[1] "Hello,World"
[1] "Hello,World"
假设不加i=i+1,i=i+2，则这个循环会一直持续下去，可以按Esc键退出
故在使用while循环时，要注意条件表达式
（3）if else结构
> score=70;if (score >60 ) {print ("Passwd") } else {print ("Failed")}
[1] "Passwd"
#简写版本
> ifelse( score >60,print ("Passwd"),print ("Failed"))
[1] "Passwd"
[1] "Passwd"
如果score是一个向量，操作起来会更容易，可以对向量中每一个值进行判断
（4）switch语句
如果条件有多种情况，不是简单的真or假，就可以使用switch语句
典型例子：判断成绩
#选择type，若是计算平均值则mean，……
> centre <- function(x, type) {
+ switch(type,
+ mean = mean(x),
+ median = median(x),
+ trimmed = mean(x, trim = .1))
+ }
> x <- rcauchy(10)
> centre(x, "mean")
[1] -0.09539974
> centre(x, "median")
[1] 0.1717119
> centre(x, "trimmed")
[1] 0.009479639
