R语言数据分析与作图示例

数据分析与作图往往是一起的，分析的过程也是作图的过程，而分析的结果往往需要通过作图展示出来，作图的前提往往是对数据进行分析。

这里用到的数据为存于期考成绩.tsv。使用如下命令读入：

CJ <- read.delim("ftp://ftp.biolab.wang/data/期考成绩.tsv", fileEncoding="GBK")

这个数据共六列，前边五列是每个学生的基本信息，第6列是某次期考的成绩。共有606名学生。那么面对这样的一堆数据，如果不经过分析和可视化，很难看出有用的信息。

数据的可视化(作图)使用ggplot2包，这里先载入。往后的的所有用到相应的命令的地方都默认了已经载入了ggplot2。

library(ggplot2)

使用boxplot能够快速地查看各班成绩的分布情况，

ggplot(data = CJ) +
    geom_boxplot(mapping = aes(x=班级, y=成绩))

这里通过boxplot能够非常直观看出来，整体上1-7班的成绩比8-12班的好。更多信息可以展示出来，比如把性别区分出来，把专业表示出来。

ggplot(data = CJ) +
    geom_boxplot(mapping = aes(x=班级, y=成绩, colour=性别))

可以通过分面将不同的专业分开，以便不同专业之间的比较 (留意这里facet_grid里设置的scale和space参数，试试去掉看看效果)：

ggplot(data = CJ) +
    geom_boxplot(mapping = aes(x=班级, y=成绩, colour=性别)) +
    facet_grid(.~专业, scales = "free_x", space = "free")

boxplot展示分布细节不够，可以附加一个小提琴图(violin plot), 这样能够更清楚地看出各班的成绩分布情况。为了能够让两个图层(layer)都展示出来，这里将boxplot的宽度改小，然后是小提琴图的填充改成透明的 (也可以把boxplot放到后边，这样就是处于上层了):

ggplot(data = CJ) +
    geom_boxplot(mapping = aes(x=班级, y=成绩), width=0.125) +
    geom_violin(mapping = aes(x=班级, y=成绩), fill="transparent") 

突然，我对学霸比较感兴趣，想看看各班超过85分的都有谁。这时候我们可以使用geom_text标注出来。

ggplot(data = CJ) +
    geom_boxplot(mapping = aes(x=班级, y=成绩), width=0.125) +
    geom_violin(mapping = aes(x=班级, y=成绩), fill="transparent") +
    geom_text(data=CJ[CJ$成绩 > 85, ], mapping = aes(x=班级, y=成绩, label=姓名))

可以看到，一些标注挤在一起了，这时候可以使用ggrepel包的geom_text_repel函数替代geom_text。

ggplot(data = CJ) +
    geom_boxplot(mapping = aes(x=班级, y=成绩), width=0.125) +
    geom_violin(mapping = aes(x=班级, y=成绩), fill="transparent") +
    ggrepel::geom_text_repel(data=CJ[CJ$成绩 > 85, ], mapping = aes(x=班级, y=成绩, label=姓名)) 

另一个大家关心的问题可能是每个班有多少人挂科了，有多人优秀，比如我们这里以60、70、80为界分为挂科、合格、中等、优秀四档，那么我们看看每个班的人数分布情况如何。首先是进行统计：

CJ$档次 <- cut(CJ$成绩, 
               breaks=c(0, 60, 70, 80, 101), 
               right = FALSE, 
               labels=c("挂科", "合格", "中等", "优秀"))
 
library(plyr)
TJ <- ddply(CJ, .(班级, 档次), summarize, 人数=length(学号))

很显然，这个统计依然有几十行，一眼看过去并不能获取各个班的情况的直观信息，这里简单用geom_bar来做个图：

ggplot(data=TJ) +
    geom_bar(mapping = aes(x=班级, y=人数, fill=档次), 
        stat="identity", 
        position=position_stack(reverse=TRUE)) +
    geom_text(mapping = aes(x=班级, y=人数, group=档次, label=人数), 
        position=position_stack(vjust=0.5, reverse=TRUE))

通过这样的一个图，我们能够非常直观的对班级之间进行比较。 要排除总人数的干扰直接看各层次的比例，可以简单得改一下参数就得到这样的图:

ggplot(data=TJ) +
    geom_bar(mapping = aes(x=班级, y=人数, fill=档次), 
        stat="identity", 
        position=position_fill(reverse=TRUE)) +
    geom_text(mapping = aes(x=班级, y=人数, group=档次, label=人数), 
        position=position_fill(vjust=0.5, reverse=TRUE)) +
    ylab("人数占比")

上图还能直接从CJ通过ggplot2自己的stat得到，可以从geom出发，也可以从stat出发，两段代码如下：

ggplot(data=CJ) +
    geom_bar(mapping = aes(x=班级, fill=档次), 
        position=position_fill(reverse=TRUE)) +
    geom_text(mapping = aes(x=班级, group=档次, label=..count..), 
        stat = "count",
        position=position_fill(vjust=0.5, reverse=TRUE)) +
    ylab("人数占比")

得到的图是完全一样的。注意其中的stat=“count”和label=..count..的用法。

ggplot(data=CJ) +
    stat_count(mapping = aes(x=班级, fill=档次),
	position=position_fill(reverse=TRUE)) +
    stat_count(mapping = aes(x=班级, group=档次, label=..count..), 
        geom = "text",
        position=position_fill(vjust=0.5, reverse=TRUE)) +
    ylab("人数占比")

得到的图是完全一样的。注意其中的stat_count和geom_bar是默认对应的，因此前一个省略了。

先对各个班的平均值（MEAN）、标准误（SEM）、95%置信区间(CI)等进行计算。

CJ <- read.delim("ftp://ftp.biolab.wang/data/期考成绩.tsv", fileEncoding="GBK")
 
PJ <- ddply(CJ, .(专业, 班级), summarize, MEAN = mean(成绩), SD=sd(成绩), N=length(学号))
PJ$CI95Error <- qnorm(0.95)*PJ$SD/sqrt(PJ$N)     # 95%置信区间是均值偏差
PJ$SEM       <- PJ$SD/sqrt(PJ$N)                 # 标准误
 
PJ  # 看看得到的数据框是什么样子的

作图，可以是柱状图来比较各班的平均分，通常还要加上误差线(error bar), 误差线有SD、SEM、95%CI三种形式，都有人在用。我这里使用SEM，作图如下：

ggplot(PJ) +
    geom_bar(aes(x=班级, y=MEAN), stat = "identity") +
    geom_errorbar(aes(x=班级, ymin = MEAN-SEM, ymax=MEAN+SEM), width=0.2) +
    scale_y_continuous(breaks=seq(0, 100, 10), expand=c(0,0)) +
    xlab(NULL) + ylab("考试成绩") +
    theme_classic()

从图上能够直观看出各个班级的平均分的高低的差别，差别是否显著则需要统计检验。这里可以使用方差分析。

aovE <- aov(成绩~班级, data=CJ)
 
summary(aovE)

得到结果如下：

             Df Sum Sq Mean Sq F value              Pr(>F)    
班级         11  12332  1121.1   18.16 <0.0000000000000002 ***
Residuals   594  36680    61.8                                
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

从上边的结果的看得出来，显然班级之间是存在显著的差异的，哪些班级之间存在显著的差异，则可以使用Post Hoc检验。

P.H.test <- TukeyHSD(aovE)
P.H.result <- data.frame(P.H.test[[1]])
head(P.H.result)      # 结果为4列，最后一列的p.adj就是我们需要的矫正后的p值

可以把结果展示在图上，用geom_segment和geom_text把检验显著的p值标注出来，这样比较直观地看出哪些组间有显著的差异。

P.H.result$班级1 <- gsub("-.*", "", row.names(P.H.result))  # 从行名中提取班级信息，班级之间是用-相连的，因此预先自己处理好班级名称，避免班级名称本身带-。如果无法避免则可以自己研究怎样把两个班级名称分别提出来。
P.H.result$班级2 <- gsub(".*-", "", row.names(P.H.result))  # 这里简单的把-及前边的内容去掉，就得到后一个班级的名称
 
P.H.result.sig <- P.H.result[P.H.result$p.adj < 0.05, ]
P.H.result.sig$y <- floor(max(PJ$MEAN + PJ$SEM)) + 5*(1:nrow(P.H.result.sig))  #增加一列，用于作图
 
# 使用因子，这样比较好的控制班级的显示顺序，也使得填上上去的p值更准确地对应
BJs <- c("2001班", "2002班", "2003班", "2004班", "2005班", "2006班", "2007班", "2008班", "2009班", "2010班", "2011班", "2012班")
PJ$班级 <- factor(PJ$班级, levels=BJs)
P.H.result.sig$班级1 <- factor(P.H.result.sig$班级1, levels= BJs)
P.H.result.sig$班级2 <- factor(P.H.result.sig$班级2, levels= BJs)
 
P.H.result.sig$x <- (as.numeric(P.H.result.sig$班级1) + as.numeric(P.H.result.sig$班级2))/2   #p值放在两个班级的中间
 
ggplot() +
    geom_bar(data=PJ, aes(x=班级, y=MEAN), stat = "identity") +
    geom_errorbar(data=PJ, aes(x=班级, ymin = MEAN-SEM, ymax=MEAN+SEM), width=0.2) +
    geom_segment(data=P.H.result.sig, aes(x=班级1, xend=班级2, y=y, yend=y)) +
    geom_text(data=P.H.result.sig, aes(x=x, y=y, label=sprintf("adjusted p = %.6f", p.adj)), vjust=-0.25) +
    scale_y_continuous(limits =c(0, max(P.H.result.sig$y)+5), breaks=seq(0, 70, 10), expand=c(0,0)) +
    xlab(NULL) + ylab("考试成绩") +
    theme_classic()

最后得到下边这个图，很容易看出两两比较哪些班级之间的平均分存在差异。

本网页已被访问6407次。