写R脚本：以对多分组数据进行正态分布检验为例

Posted 2021-02-13 yanlingbin

问题提出

正态分布检验一次只能检验一个分组，如果有多组数据需要检验，则需要运行多次

解决思路

• 使用循环命令可以实现按一定规则计算。
• 如果以后也经常需要使用，写成脚本调用更方便些，需要使用的时候直接调用即可。
• 脚本针对的场景相对直接使用循环命令更广泛写，如果仅使用循环命令，很多参数与类型定义直接使用数据对应的即可，因此写脚本难度相对高一些，但是设计出可以通用的脚本也是小小的成就。

设计框架

函数名与参数

从分组是否符合生态分布需要输入数据出发，定义的函数设计以下内容：

• 函数名，这个无所谓，只要不是数字开头，不含符号字符即可。先做个简单的脚本，只用Shapiro-Wilk法进行检验，因为要检验多个，就假设multi，关键词之间，可以用“.”隔开，函数名shapiro.test.multi。
• 数据集，数据类型多样，通常是数据框。
• 分组变量名，需要获取分组变量不重复的值，分组的数据可能是文本型也可能是数值，通常是文本型。所以先只定义文本型的。
• 计算值变量名，根据不同的分组名提取此值为向量，此数据必须为数值。

输出数据框

从输出来看，需要得到一张表，这个表要包含计算得到的重要信息:

• 分组名
• 统计量，对于Shapiro-Wilk来说，即W
• 概率值，即p.value
• 分布的判断结果，如果p.value > 0.05，则判断符合正态分布

定义用于存储数据的数据框

• 数据框应包含输出数据框所需变量信息
• 定义为空的数据框以便填入数据

使用管道函数简洁代码

• 需要加载管道函数的程序包magrittr，使用require()，判断是否已经加载，如未加载则加载。

shapiro.test.multi <- function(   #定义函数名
  data,   #定义函数第一个参数
  value,  #第2参数
  group)  #第3参数
  {       #开始计算
  
  require(magrittr)   #按需要加载管道函数包
  
  table(data[,group]) %>%   #提取分组信息，此处即为统计group中值出现的次数，达到了去重的目的
    data.frame(.) -> a1   #将提取信息从table格式转为数据库data.frame并存为a1，这样才能提取其中一列转为向量
  
  a2 <- as.vector(a1[,1])  #将a1数据的第一列转为向量，这一列即为不重复的分组信息
    
  data = data.frame(group = data[,group],  #对数据集进行关键变量的提取，提取分组变量为data数据集的group变量
                    value = data[,value])  #提取计算值为data数据集的value
  
  test.result <- data.frame(No=0,        #行号
                            Name=0,      #分组名
                            W=0,         #W值
                            p.value=0,   #p值
                            norm.test=0) #检测结果

  for (i in (1:length(a2))){     #定义for循环计算，从1到a2的长度值这样一个区间，a2的长度即为分组的数量
    subset(data,                 #指定取数据集    换行写代码使层次更清晰
           group == a2[i],       #定义去子集的条件，“==”为判断
           select = value) %>%   #定义需要取集的变量/列，“=”为定义
      .[,1] %>%                  #  "."定义计算结果放置的位置
      shapiro.test(.) -> t.r     #进行正态检验分布并存储为t.r
    
    test.result[i,1] = i              #存储组序号
    test.result[i,2] = a2[i]          #存储分组名
    test.result[i,3] = t.r$statistic  #存储W统计量
    test.result[i,4] = t.r$p.value    #存储计算的p值
    
    if      #if判断
    (t.r$p.value > 0.05)           #判断条件
      test.result[i,5] = "Norm"    #通过判断后的命令
      else 
        test.result[i,5] = "Other_situation"  #未通过判断执行的命令
    } #结束循环计算
  
   test.result[nrow( test.result)+1,1] = "Test Method:"  #给数据框加上检验正态分布方法信息，在最后一行之后加上一行，在第1列放入次数据
   test.result[nrow( test.result),2] = "Shapiro-Wilk"    #同上行，存在第二列
   
  test.result  #显示用于存储计算结果的数据框
}              #脚本结束

此处以iris数据为例，计算测量的三个物种的花瓣长度值是否都符合正态分布假设。

shapiro.test.multi(iris,value = "Sepal.Length",group = "Species")

运行结果

-	No	Name	W	p.value	norm.test
1	1	setosa	0.9776985	0.4595132	Norm
2	2	versicolor	0.9778357	0.4647370	Norm
3	3	virginica	0.9711794	0.2583147	Norm
4	Test Method:	Shapiro-Wilk	NA	NA	NA