5991 1

[程序分享] R语言自定义相同指标内多组数据对比:单因素方差分析函数 [推广有奖]

  • 5关注
  • 8粉丝

博士生

38%

还不是VIP/贵宾

-

威望
0
论坛币
3053 个
通用积分
73.0212
学术水平
37 点
热心指数
36 点
信用等级
23 点
经验
33617 点
帖子
298
精华
0
在线时间
141 小时
注册时间
2012-11-23
最后登录
2022-9-18

+2 论坛币
k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群
赵安豆老师微信:zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

感谢您参与论坛问题回答

经管之家送您两个论坛币!

+2 论坛币

自定义相同指标内多组数据对比:单因素方差分析函数

注意:无论以哪种形式作为数据输入,符号~右边表示分类变量,符合左边表示因变量,即要研究的指标数据

aov.func<-function(data,alpha=0.05,plot.logic=T,p.adjust.mod="holm"){

+ error<-FALSE

+ data[,1]<-as.factor(data[,1])

+ group.levels<-as.numeric(levels(data[,1]))#分类变量的水平取值

+ for(i in group.levels){

+ if(shapiro.test(data[which(data[,1]==i),2])$p.value<alpha){+ print(paste("ERROR:",i,"组数据不服从正态分布"))

#验证不同水平下的指标是否服从正态分布,只要有一组不服从正态分布,就带引第i组不服从正态分布

+ error<-TRUE

+ }

+ }

+ if(error){

+ return()

+ }

+ else{

+ print("符合正态性前提")

+ }

+ #检验不同水平下指标的方差齐次性

+ if(bartlett.test(data[,2]~data[,1])$p.value<alpha){

+ print("ERROR:符合方差齐次性前提")

+ }

#验证不同水平下指标数据的方差是否相同,如果p-value大于alpha表示服从,小于则表示不服从

+ #绘制箱线图

+ if(plot.logic){

+ boxplot(data[,2]~data[,1])

+ }

+ #方差分析

+ sol<-aov(data[,2]~data[,1])

+ p.value<-summary(sol)[[1]][,5][1]

+ tab<-matrix(NA,nrow=3,ncol=5)

+ dimnames(tab)<-list(c("分类变量A","随机性误差","总和"),c("自由度","差异平方和","均方","统计量F","P-value"))

+ tab[1:2,1]<-summary(sol)[[1]][,1];tab[3,1]<-sum(tab[1:2,1])#自由度

+ tab[1:2,2]<-summary(sol)[[1]][,2];tab[3,2]<-sum(tab[1:2,2])#差异平方和

+ tab[1:2,3]<-summary(sol)[[1]][,3];#均方

+ tab[1:2,4]<-summary(sol)[[1]][,4];#统计量F

+ tab[1,5]<-p.value;#p-value

+ print("====方差分析====")

+ print(tab)

+ #多重T检验

+ if(p.value<alpha){

+ print(paste(p.value,"<",alpha,":各水平下的指标有明显差别"))

+ sol.t<-pairwise.t.test(data[,2],data[,1],p.adjust.method=p.adjust.mod)

+ tab<-sol.t[[3]]

+ tab[which((sol.t[[3]]>alpha)==TRUE)]<-"无差别"

+ tab[which((sol.t[[3]]>alpha)==FALSE)]<-"显著差别"

+ print("====多重T检验====")

+ print(tab)

+ }

+ else{

+ print(paste(p.value,">",alpha,":个水平下的指标无明显差别"))

+

+ }

+ return(sol)

+ }

> 自定义函数aov.func的参数说明:

data:可以是matrixdata.frame两种对象,第一列表示不同水平下的取值,第二列是对应的指标数据,例如:kpi<-c(rnorm(20,80,5),rnorm(25,195,5),rnorm(23,256,5)),group<-c(rep(1,20),rep(2,25),rep(3,23))

data<-data.frame(group,kpi),

Alpha表示显著性水平,默认为0.05

Plot.logic如果去T默认,时绘制data的箱线图,如果去F时不进行绘图。

P.adjust.mod设置在左多重T检验时,修正p-value的值。

例:数据为某电商2个年龄段的日客单价数据。

>  price1<-c(79,79,87,79,71,84,82,82,85,82,88,81,71,76,81,75,72,83,76)

>price2<-c(193,192,191,181,191,192,187,191,196,196,196,196,192,194,192,189,196,202,197,206,199,191,194,195,196)

> date<-rep(1:2,c(19,25))

> price<-c(price1,price2)

> data<-data.frame(group=date,x=price)

> data

   group   x

1      1  79

2      1  79

3      1  87

4      1  79

5      1  71

6      1  84

7      1  82

8      1  82

9      1  85

10     1  82

11     1  88

12     1  81

13     1  71

14     1  76

15     1  81

16     1  75

17     1  72

18     1  83

19     1  76

20     2 193

21     2 192

22     2 191

23     2 181

24     2 191

25     2 192

26     2 187

27     2 191

28     2 196

29     2 196

30     2 196

31     2 196

32     2 192

33     2 194

34     2 192

35     2 189

36     2 196

37     2 202

38     2 197

39     2 206

40     2 199

41     2 191

42     2 194

43     2 195

44     2 196

> sol<-aov.func(data)

[1] "符合正态性前提"

[1] "====方差分析===="

           自由度 差异平方和         均方  统计量F      P-value

分类变量A       1 140712.579 140712.57895 5780.327 1.289547e-46

随机性误差     42   1022.421     24.34336       NA           NA

总和           43 141735.000           NA       NA           NA

[1] "1.28954695480505e-46 < 0.05 :各水平下的指标有明显差别"

[1] "====多重T检验===="

  1         

2 "显著差别"

>


二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

关键词:单因素方差分析 方差分析 R语言 自定义 单因素 因变量 内多

已有 1 人评分经验 论坛币 学术水平 热心指数 收起 理由
dxystata + 100 + 20 + 1 + 1 精彩帖子

总评分: 经验 + 100  论坛币 + 20  学术水平 + 1  热心指数 + 1   查看全部评分

沙发
liquidambar 发表于 2018-1-12 11:08:55 |只看作者 |坛友微信交流群
summary(aov(x~group,dat=data))
已有 1 人评分论坛币 收起 理由
admin_kefu + 10 热心帮助其他会员

总评分: 论坛币 + 10   查看全部评分

使用道具

您需要登录后才可以回帖 登录 | 我要注册

本版微信群
加好友,备注cda
拉您进交流群

京ICP备16021002-2号 京B2-20170662号 京公网安备 11010802022788号 论坛法律顾问:王进律师 知识产权保护声明   免责及隐私声明

GMT+8, 2024-4-20 00:26