《生物统计学》上机内容（三）

教学内容：

本节课主要介绍如何运用R语言进行单因素方差检验。
以下内容仅用于北京林业大学生态与自然保护学院《生物统计学》的第三节上机内容。

参考资料:

https://rcompanion.org/handbook/
《R语言与统计分析》（汤银才著）

数据下载链接：

lesson_one_way_aov.csv

R语言教程:

# 清理内存
rm(list = ls())
gc()

##          used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 461149 24.7     995742 53.2   643651 34.4
## Vcells 828881  6.4    8388608 64.0  1649943 12.6

# 加载 packages
# 如果packages加载没成功
# install.packages("agricolae")
# install.packages("car")
# install.packages("psych")
# install.packages("skimr")
library(agricolae)

## Warning: 程辑包'agricolae'是用R版本4.1.1 来建造的

library(car)

## 载入需要的程辑包：carData

library(psych)

## Warning: 程辑包'psych'是用R版本4.1.1 来建造的

##
## 载入程辑包：'psych'

## The following object is masked from 'package:car':
##
##     logit

library(skimr)

## Warning: 程辑包'skimr'是用R版本4.1.1 来建造的

# 加载工作区域
# 需要在setwd里面设置自己的文件路径
setwd("F:/2018年生物统计课程/数据")
getwd()

## [1] "F:/2018年生物统计课程/数据"

单因素方差检验 One-way ANOVA

练习题1：用6种培养液培养红苜蓿，每一种培养液做5次重复，测定5盆苜蓿的含氮量，结果如下表（单位: mg）。问用6种不同培养液培养的红苜蓿含氮量差异是否显著？

Table 1. Cultivation type

	tp_01	tp_02	tp_03	tp_04	tp_05	tp_06
1	19.4	17.7	17	20.7	14.3	17.3
2	32.6	24.8	19.4	21	14.4	19.4
3	27	27.9	19.1	20.5	11.8	19.1
4	32.1	25.2	11.9	18.8	11.6	16.9
5	33	24.3	15.8	18.6	14.2	20.8

# 加载数据
n_content.data <- read.csv("F:/2018年生物统计课程/数据/lesson_one_way_aov.csv")

# 将treatment变为处理因子
n_content.data$treatment <- as.factor(n_content.data$treatment)

# 对数据进行了log转化，以满足方差齐性和正态性
# 同时，将数据保存到n_content.log
n_content.data$n_content.log <- log10(n_content.data$n_content)

# 统计描述：summary功能
# 查看数据格式
head(n_content.data) # 展示数据前六行

##   treatment n_content n_content.log
## 1     tp_01      19.4      1.287802
## 2     tp_01      32.6      1.513218
## 3     tp_01      27.0      1.431364
## 4     tp_01      32.1      1.506505
## 5     tp_01      33.0      1.518514
## 6     tp_02      17.7      1.247973

str(n_content.data) # 展示数据格式

## 'data.frame':    30 obs. of  3 variables:
##  $ treatment    : Factor w/ 6 levels "tp_01","tp_02",..: 1 1 1 1 1 2 2 2 2 2 ...
##  $ n_content    : num  19.4 32.6 27 32.1 33 17.7 24.8 27.9 25.2 24.3 ...
##  $ n_content.log: num  1.29 1.51 1.43 1.51 1.52 ...

# n_content.data$n_content.log的数据汇总
summary(n_content.data$n_content.log)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
##   1.064   1.229   1.284   1.288   1.370   1.519

# 统计描述的其他选择1，skimr包
# install.packages("skimr")
# library(skimr)
skim(n_content.data)

Data summary

Name	n_content.data
Number of rows	30
Number of columns	3
————————
Column type frequency:
factor	1
numeric	2
————————
Group variables	None

Variable type: factor

skim_variable	n_missing	complete_rate	ordered	n_unique	top_counts
treatment	0	1	FALSE	6	tp: 5, tp: 5, tp: 5, tp: 5

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
n_content	0	1	20.22	5.90	11.60	16.92	19.25	23.48	33.00	▅▇▃▃▂
n_content.log	0	1	1.29	0.12	1.06	1.23	1.28	1.37	1.52	▂▃▇▂▃

# 统计描述的其他选择2，psych包
# 参考：https://www.cnblogs.com/Clactionxy/p/9865196.html
# library(psych)
des.result <- describeBy(n_content.data$n_content.log,
                         n_content.data$treatment,
                         mat= T,
                         digits = 3)

des.result

##     item group1 vars n  mean    sd median trimmed   mad   min   max range
## X11    1  tp_01    1 5 1.451 0.098  1.507   1.451 0.018 1.288 1.519 0.231
## X12    2  tp_02    1 5 1.375 0.075  1.394   1.375 0.013 1.248 1.446 0.198
## X13    3  tp_03    1 5 1.215 0.086  1.230   1.215 0.075 1.076 1.288 0.212
## X14    4  tp_04    1 5 1.299 0.025  1.312   1.299 0.016 1.270 1.322 0.053
## X15    5  tp_05    1 5 1.120 0.048  1.152   1.120 0.009 1.064 1.158 0.094
## X16    6  tp_06    1 5 1.271 0.037  1.281   1.271 0.055 1.228 1.318 0.090
##       skew kurtosis    se
## X11 -0.781   -1.358 0.044
## X12 -0.803   -1.162 0.033
## X13 -0.632   -1.457 0.038
## X14 -0.252   -2.210 0.011
## X15 -0.294   -2.238 0.021
## X16  0.005   -1.997 0.017

方差分析过程

# 对数据进行线性拟合
#library(car)
fit.lm <- lm(n_content.log ~ treatment, data = n_content.data)

# 方差齐性检验
# 使用bartlett.test进行方差齐性检验
leveneTest(n_content.log ~ treatment, data = n_content.data)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  5  0.4761 0.7904
##       24

# 使用qqPlot对数据进行qq图绘制
qqPlot(fit.lm, # 拟合方程
       simulate = TRUE, # 分布模拟？
       labels = TRUE, # 对异常值进行标记？
       col = "red", # 点的颜色
       main = "qqplot", # 总标题
       ylab = "Studentized residual" # y轴标题
       )

## [1]  1 14

# 分组进行qq图绘制
#  分割一个1行6列的绘图空间
op <- par(mfcol = c(1,6))

# 自定义qq.plot
qq.plot <- function(x){qqnorm(x); qqline(x, col = "red")}

# 使用tapply进行批量处理qq.plot图
tapply(n_content.data$n_content.log, n_content.data$treatment, qq.plot)

## $tp_01
## NULL
##
## $tp_02
## NULL
##
## $tp_03
## NULL
##
## $tp_04
## NULL
##
## $tp_05
## NULL
##
## $tp_06
## NULL

# 关闭绘图空间
par(op)

# 正态性检验
# 当分析小于50行的小样本数据时，用Shapiro-Wilk 检验
# 当分析大于50行的大样本数据时，用Kolmogorov-Smirnov检验
# shapiro.test(x)
# tapply 进行批量处理shapiro.test
tapply(n_content.data$n_content.log, n_content.data$treatment, shapiro.test)

## $tp_01
##
##  Shapiro-Wilk normality test
##
## data:  X[[i]]
## W = 0.77768, p-value = 0.05264
##
##
## $tp_02
##
##  Shapiro-Wilk normality test
##
## data:  X[[i]]
## W = 0.81895, p-value = 0.1146
##
##
## $tp_03
##
##  Shapiro-Wilk normality test
##
## data:  X[[i]]
## W = 0.87484, p-value = 0.2866
##
##
## $tp_04
##
##  Shapiro-Wilk normality test
##
## data:  X[[i]]
## W = 0.8223, p-value = 0.1216
##
##
## $tp_05
##
##  Shapiro-Wilk normality test
##
## data:  X[[i]]
## W = 0.74387, p-value = 0.02609
##
##
## $tp_06
##
##  Shapiro-Wilk normality test
##
## data:  X[[i]]
## W = 0.934, p-value = 0.6239

# step1：先使用lm进行线性拟合
# fit.lm <- lm(n_content.log ~ treatment, data = n_content.data)

# step2: 再使用car包里面Anova进行方差分析
# SPSS默认的方差分析类型为III
# library(car)
Anova(fit.lm, type = "III")

## Anova Table (Type III tests)
##
## Response: n_content.log
##              Sum Sq Df  F value    Pr(>F)
## (Intercept) 10.5340  1 2348.945 < 2.2e-16 ***
## treatment    0.3408  5   15.198 9.111e-07 ***
## Residuals    0.1076 24
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# step3: 使用Tukey检验进行多重比较
# install.packages("agricolae")
# https://myaseen208.github.io/agricolae/
# library(agricolae)
comparison <- HSD.test(fit.lm, # 拟合的模型
                       "treatment", # 处理因子
                       group = TRUE,
                       main = "Tukey comparison")

# 展示Tukey分析的结果
comparison

## $statistics
##       MSerror Df    Mean       CV       MSD
##   0.004484557 24 1.28849 5.197311 0.1309542
##
## $parameters
##    test    name.t ntr StudentizedRange alpha
##   Tukey treatment   6         4.372651  0.05
##
## $means
##       n_content.log        std r      Min      Max      Q25      Q50      Q75
## tp_01      1.451480 0.09814285 5 1.287802 1.518514 1.431364 1.506505 1.513218
## tp_02      1.375007 0.07468031 5 1.247973 1.445604 1.385606 1.394452 1.401401
## tp_03      1.214698 0.08604207 5 1.075547 1.287802 1.198657 1.230449 1.281033
## tp_04      1.298723 0.02487983 5 1.269513 1.322219 1.274158 1.311754 1.315970
## tp_05      1.120465 0.04785924 5 1.064458 1.158362 1.071882 1.152288 1.155336
## tp_06      1.270566 0.03722129 5 1.227887 1.318063 1.238046 1.281033 1.287802
##
## $comparison
## NULL
##
## $groups
##       n_content.log groups
## tp_01      1.451480      a
## tp_02      1.375007     ab
## tp_04      1.298723     bc
## tp_06      1.270566     bc
## tp_03      1.214698     cd
## tp_05      1.120465      d
##
## attr(,"class")
## [1] "group"

# 使用boxplot函数建立盒装图
boxplot(n_content.log ~ treatment, # 因变量与自变量之间的关系 y ~ x
        data = n_content.data, # 数据来源
        ylim = c(0.5, 2), # y轴范围
        xlab = "Cultivation type", # y轴标题
        ylab = "log10(N content) (mg)" # x轴标题
        )

# comparison$groups保存有用于作图的信息
comparison$groups

##       n_content.log groups
## tp_01      1.451480      a
## tp_02      1.375007     ab
## tp_04      1.298723     bc
## tp_06      1.270566     bc
## tp_03      1.214698     cd
## tp_05      1.120465      d

# 使用comparison$groups添加文本
# 首先确定x轴位置
comparison$groups$x.pos <- c(1, 2, 4, 6, 3, 5)

# 添加comparison$groups的Tukey结果
text(x = comparison$groups$x.pos, # 字符的x轴位置
     y = comparison$groups$n_content.log + 0.3, # 字符的y轴位置
     labels = comparison$groups$groups # 输出字符
      )

教学内容：

参考资料:

数据下载链接：

R语言教程:

单因素方差检验 One-way ANOVA

方差分析过程

你可能也喜欢

《生物统计学》上机内容（五）

《生物统计学》-投币游戏

《生物统计学》上机内容（四）