《生物统计学》上机内容(三)

教学内容:

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

参考资料:

数据下载链接:

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 # 输出字符
      )

0 评论
内联反馈
查看所有评论