jamesbang

V1

2022/08/27阅读:20主题:雁栖湖

🤔 多组比较及绘图

1写在前面

写毕业课题统计时编写的一段代码,大量数据很快就可以统计出结果并作用,方便的很。 统计使用的是r基础stat包,绘图使用的ggplot2包。都是很常见的,网上教程也很多。

2示例数据

示例为利用excel随机生成的一列数字

y <- read.table("clipboard", header = F)

该法是直接访问的剪贴板,可以用read.xlsxread.tableread.csv等函数读取已经整理好的数据。 分组信息也可读入。

3输入分组信息

a1 <- factor(c(rep(c('模型组','低剂组','高剂组'),
each = 5)))
a2 <- factor(c(rep(c('空白组','模型组','治疗组'),
time = 3)))

eachtime参数的不同,可以自己运行一下。

4计算各组均值和标准差

mean <- tapply(y$V1,a1,mean)
sd <- tapply(y$V1,a1,sd)
data <- data.frame(mean,sd)
#排序
data$group <- factor(row.names(data), levels = c('模型组','低剂组','高剂组'))

5正态性检验

tapply(y$V1,a1,chisq.test)

6方差齐性检验

bartlett.test(y$V1,a1)

7单因素方差分析

result.aov<-aov (y$V1~a1)
summary(result.aov)
# 两两比较
pairwise.t.test(y$V1,a1,p.adjust.method = 'bonferroni')

8另一种多重比较方法

TukeyHSD(result.aov)
plot(TukeyHSD(result.aov))

9非参数检验

kruskal.test(y$V1,a1)
#两两比较
pairwise.wilcox.test(y$V1,a1,p.adjust.method = "bonferroni")

10使用sink()函数获取结果

sink("calc.txt")
print(data)
print(tapply(y$V1,a1,chisq.test))
print(bartlett.test(y$V1,a1))
#多组数据符合正态分布,且方差齐则使用单因素方差分析。
result.aov<-aov (y$V1~a1)
print(summary(result.aov))
print(pairwise.t.test(y$V1,a1,p.adjust.method = 'bonferroni'))
#多组数据不符合正态分布;或虽符合正态分布,但方差不齐,则使用非参数检验。
print(kruskal.test(y$V1,a1))
print(pairwise.wilcox.test(y$V1,a1,p.adjust.method = 'bonferroni'))
sink()

输出结果

      mean       sd  group
低剂组 14.2 2.167948 低剂组
高剂组 25.0 1.581139 高剂组
模型组 7.2 3.420526 模型组
$低剂组

Chi-squared test for given probabilities

data: X[[i]]
X-squared = 1.3239, df = 4, p-value = 0.8573


$高剂组

Chi-squared test for given probabilities

data: X[[i]]
X-squared = 0.4, df = 4, p-value = 0.9825


$模型组

Chi-squared test for given probabilities

data: X[[i]]
X-squared = 6.5, df = 4, p-value = 0.1648



Bartlett test of homogeneity of variances

data: y$V1 and a1
Bartlett's K-squared = 2.1535, df = 2, p-value = 0.3407

Df Sum Sq Mean Sq F value Pr(>F)
f 2 804.1 402.1 63.82 4.03e-07 ***
Residuals 12 75.6 6.3
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Pairwise comparisons using t tests with pooled SD

data: y$V1 and a1

低剂组 高剂组
高剂组 5.7e-05 -
模型组 0.0026 3.1e-07

P value adjustment method: bonferroni
Kruskal-Wallis rank sum test

data: y$V1 and a1
Kruskal-Wallis chi-squared = 12.545, df = 2, p-value = 0.001888

Pairwise comparisons using Wilcoxon rank sum test with continuity correction

data: y$V1 and a1

低剂组 高剂组
高剂组 0.036 -
模型组 0.035 0.036

P value adjustment method: bonferroni

11绘制条形图加误差线

library(ggplot2)
p<-ggplot(data,aes(group,mean))+
geom_col(
fill = c('white','black','gray80'),
col ='black',
width = 0.6,
position = position_dodge(0.8),
lwd = 1.5,
)+
theme_classic(
base_size = 32,
base_family = 'serif'#字体
)+
scale_y_continuous(
expand = c(0,0)
)+
scale_x_discrete(
expand = c(0.2,0.2)
)+
geom_blank(aes(y= (mean+sd)*1.2))+
geom_errorbar(
aes(ymin = mean, ymax = mean+sd),
width = 0.3,lwd = 1.5
)+
xlab("这是x轴")+
ylab("这是y轴")
p

12把图保存下来

tiff('barplot.tif')
p
dev.off()

13小提琴图,加误差线,不要图例

小提琴图和箱线图用到的是所有数据,需要构建包含所有数据的表格。

data2 <- as.data.frame(y$V1)
data2$fenzu <- a1
#注意两个表中的同样数据列名一致,要不会报错。
colnames(data2) <-c('mean','group')
#排序
data2$group <- factor(data2$group, levels = c('模型组','低剂组','高剂组'))

p2<-ggplot(data = data2,aes(group,mean,fill=group))+
geom_violin(trim=FALSE,color="black",
lwd = 1.5
)+
theme_classic(
base_size = 32,
base_family = 'serif'
)+
geom_blank(data = data, aes(y= (mean+sd)*1.2))+
geom_point(data = data, aes(group,mean),pch = 18, size = 6)+
geom_errorbar(data = data,
aes(ymin = mean-sd, ymax = mean+sd),
width = 0.1,lwd = 1.5
)+
xlab("这里是x轴")+
ylab("这里是y轴")+
theme(legend.position="none")
p2

14小提琴图加箱线图

p3<- ggplot(data2, aes(group, mean ,fill = group))+
geom_violin(trim=FALSE,color="black",
lwd = 1.5
)+
geom_boxplot(width=0.2,position=position_dodge(0.9),lwd= 1.5)+
theme_classic(
base_size = 32,
base_family = 'serif'
)+
geom_blank(data = data, aes(y= (mean+sd)*1.2))+
xlab("这里是x轴")+
ylab("这里是y轴")+
theme(legend.position="none")
p3

香蕉
香蕉
最后祝大家早日不卷!~

点个在看吧各位~ ✐.ɴɪᴄᴇ ᴅᴀʏ 〰

分类:

后端

标签:

后端

作者介绍

jamesbang
V1