R语言分析单身狗和舔狗哪个更短命？

一、载入数据与概览
二、单因素方差分析
- 2.1 假设检验
  - 2.1.1 是否符合正态分布
正态分布检验
- 2.2 单因素方差检验
- 2.3 多重比较
参考资料

你猜，做单身狗和做舔狗，到底哪个命更短？
一只果蝇的一辈子平均寿命 60 天，有一个因素还会让它们命更短，那就是，能不能交配。
频繁和雄果蝇交配而受孕的雌性果蝇很容易短命，而长期单身不能交配的雄性，也会短命。如果让一只雄性和八只雌性待在一起，但这八只都不能与雄性果蝇交配，那这只雄性果蝇会比其他单身狗还要少活 20 多天……

人们如何得出这样的结论？
——通过实验和数据分析。

一、载入数据与概览

该实验者把 125 只雄果蝇分成五组。每组进行分配不同数目的雌果蝇。因为新受精的雌性通常至少两天内不会再交配，所以可以借此来控制雄果蝇能否与其交配。

G1：纯独居
G2：一夫一妻（1 virgin female fruit fly）
G3：一夫八妻（8 virgin female fruit fly）
G4：一夫一妻（1 pregnant female fruit fly，不能交配）
G5：一夫八妻（8 pregnant female fruit fly，不能交配）

数据分组情况（图片来源于Topbook）

The data set is from Hanley, J. A. and S. H. Shapiro (1994) “Sexual activity and the lifespan of male fruit flies: a dataset that gets attention.” Journal of Statistics Education 2.

查看数据，这里一共两列数据，第一列是果蝇活了多少天，第二列告诉我们它来自哪一组。
Fruitfly.csv

# 查看数据
> Fruitfly.df = read.csv("Fruitfly.csv",header = T)
> View(Fruitfly.df[,1:2])

Fruitfly.df[,1:2]

确认每组都有 25 只果蝇。

# 查看每组数目
> table(Fruitfly.df$group)
G1 G2 G3 G4 G5
25 25 25 25 25

使用箱式图快速概览各组存活天数情况。这群兄弟中有只活了不到 10 天的短命鬼，也有活了快百天的长寿老蝇。但活得太长或太短，都没有什么参考价值，盒式图就是要让我们来看盒子跟黑线的位置。

# 箱式图
> ggplot(data =Fruitfly.df,aes(x = group,y =days,fill = group))+geom_boxplot()+labs(
    x = "Group",
    y = "Life Span",
    title = "Lifespan of Each Group",
  )

箱式图：横坐标表示第一二三四五组，纵坐标表示果蝇的寿命。每个盒子中间的黑线表示中位数，上下须线代表非离常值的最大值

很明显，整体看来，第四组和第五组的果蝇更短命，尤其是第五组。它们分别是一夫一妻和一夫八妻，并且所有的雌性果蝇都不能和为爱鼓掌。不仅如此，从数据看起来，面对八只雌性果蝇却不能为爱鼓掌的第五组雄性果蝇比第四组还短命。

# 查看每组的平均数
> tapply(Fruitfly.df$days, Fruitfly.df$group, mean)
G1    G2    G3    G4    G5 
63.56 64.80 63.36 56.76 38.72

这些数据是不是已经意味着，做舔狗比做单身狗还惨呢？
——不能这么草率地下结论，我们还得进行科学的统计分析来判断。

二、单因素方差分析

首先介绍下单因素方差分析。

单因素方差分析可以看成基础统计中两样本 t 检验的一个推广，要比较试验观测值的某个因变量（称为“指标”）按照一个分组变量（称为“因素”）分组后，各组的因变量均值有无显著差异。 “单因素”是指按照一种分组方式分组。
H0 假设为各组均值相等：

$R语言分析单身狗和舔狗哪个更短命？ - 图4$

如果拒绝，希望找出哪些组的均值两两之间是有显著差异的。
模型中假定各个观测独立，误差服从正态分布，且误差的方差相等
这样的问题就称为单因素方差分析问题。

以本次研究为例，我们研究是就是雄果蝇在不同组别的寿命差异，便是一个单因素方差分析问题。

2.1 假设检验

单因素方差分析的适用条件：
（1）各研究对象独立随机抽样
（2）因变量服从正态分布
（3）方差齐性

为此我们需要检验因变量是否服从正态分布、方差齐性

2.1.1 是否符合正态分布

我们先来检验数据是否符合正态分布，有两种方法使用 qqPlot() 或 shaprio.test() 检验数据是否符合正态分布

QQ plot
QQ plot 也就是 Quantile-Quantile Plots。是通过比较两个概率分布的分位数对这两个概率分布进行比较的概率图方法。如果一组实数服从正态分布，正态概率图就是一条直线。（详见 R 语言学习笔记—单因素方差分析（avo））
```
> if(!require(car))install.packages("car")
> qqPlot(lm(days ~ group, data = Fruitfly.df),simulate=TRUE,main="Q-Q plot",labels=FALSE)
```
在这里，数据主要分布在 95% 置信区间内，假设成立，属于正态分布。
shapiro 检验 ```r

正态分布检验

shapiro.test(Fruitfly.df$days) Shapiro-Wilk normality test

data: Fruitfly.df$days W = 0.9899, p-value = 0.4946

<br />H0 假设是数据符合正态分布，p>0.05,原假设成立，所以数据符合正态分布。 
<a name="p2wU3"></a>
#### 2.1.2 是否符合方差齐性
为了后面可以进行多重比较，我们还需检验这些组间的数据方差是否齐性，也就是通常我们所说的方差齐性检验。采用 `bartlett.test()` 函数
```r
> bartlett.test(days ~ group, data = Fruitfly.df)
Bartlett test of homogeneity of variances
data:  days by group
Bartlett's K-squared = 2.4196, df = 4, p-value = 0.6591

H0 假设是这几组方差没有显著不同
这里 P>0.05,表明 H0 假设成立。

2.2 单因素方差检验

经过上面的检验之后，终于可以开始进行单因素方差检验
单因素方差分析检验组间差异性，采用 avo( )函数。

> aov_model <- aov(days ~ group, data = Fruitfly.df)
> summary(aov_model)
               Df Sum Sq Mean Sq F value   Pr(>F)  
group         4  11939  2984.8   13.61 3.52e-09 ***
Residuals   120  26314   219.3                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

H0 假设为组间均值没有差异
在这里 p=3.52e-09 远远小于 0.01，所以拒绝 H0 假设，组间有差异。
但 ANOVA 的结果只能说明 there is a significant difference among the groups，但无法得知具体是哪些组之间存在 difference。为此，需要进行多重比较。

2.3 多重比较

进行多个假设检验（如均值比较）的操作称为“多重比较”（multiple comparison，或 multiple testing）。多重比较的目的就是组间进行两两的比较。
R 软件本身携带的的 TukeyHSD( ) 可以形成组间对比，对于比较结果我们再筛选出 p adj<0.05 的组间比较。

> cis=TukeyHSD(aov_model)$group
> cis[cis[,4]<0.05,]
        diff       lwr        upr        p adj
G5-G1 -24.84 -36.44047 -13.239532 2.958395e-07
G5-G2 -26.08 -37.68047 -14.479532 7.232306e-08
G5-G3 -24.64 -36.24047 -13.039532 3.701044e-07
G5-G4 -18.04 -29.64047  -6.439532 3.240300e-04

从返回结果中我们可以看出，第五组的兄弟们比第一组平均少活了 24.84 天，通常要少活个 13.2 到 36.4 天不等。后面三行的数据告诉我们，第五组分别比第二、三、四组平均少活了 26.08 天、24.64 天、18.04 天……

现在可以得出结论，第五组的兄弟确实显著比其他组的寿命都短。

所以雄性果蝇们如果想活得长，要不就得跟雌性果蝇性福地生活在一起，要不就干脆做单身狗。如果做了一只求而不得的舔狗，结局就可能是英年早逝……

library(ggplot2)
library(tidyverse)
mean.df<-Fruitfly.df %>% group_by(group) %>% summarise(mean = mean(days, na.rm = T))
ggplot(mean.df, aes(x = group, y =mean,fill=group)) + geom_bar(stat = "identity",alpha=0.8)+ geom_text(aes(x = group, y =mean, label = mean),hjust = 1.5, color = "white")+coord_flip()+labs(
    x = "Group",
    y = "Mean",
    title = "Mean Lifespan of Each Group",
  ) 
ggsave(p,filename = "Output.png.png",width = 12,height = 9)

R语言分析单身狗和舔狗哪个更短命？ - 图6
当然，统计学是一门非常严谨的学科。这些数据只能说明在这组数据中面对八只雌性却不能为爱鼓掌的雄性更短命，要知道其它组果蝇的寿命差异，还得需要更多的数据和更多的分析。

参考资料

单身狗和舔狗哪个更短命？ (qq.com)
Hanley, J. A. and S. H. Shapiro (1994) “Sexual activity and the lifespan of male fruit flies: a dataset that gets attention.” Journal of Statistics Education 2.
Lab 9 Handout Solution (umn.edu)
R 语言学习笔记—单因素方差分析（avo） - 知乎 (zhihu.com)