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

教学内容

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

参考资料

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

R语言教程

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

##          used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 431956 23.1     912334 48.8   643651 34.4
## Vcells 795200  6.1    8388608 64.0  1868282 14.3

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

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

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

## [1] "F:/2018年生物统计课程/数据/R语言代码"

单样本t检验 One-sample t-test

练习题1：某种元件的寿命X平均为250小时，且服从正态分布。某批次16只元件的寿命如下, 159, 280, 101, 212, 224, 379, 179, 264, 222, 362, 168, 250, 149, 260, 485, 170。试问是否有理由认为该批次元件是正常批次？

# 加载数据
Input = ("
id time
1 159
2 280
3 101
4 212
5 224
6 379
7 179
8 264
9 222
10 362
11 168
12 250
13 149
14 260
15 485
16 170
")

time.data <- read.table(textConnection(Input), header = TRUE)

# 查看数据格式
head(time.data) # 展示数据前六行

##   id time
## 1  1  159
## 2  2  280
## 3  3  101
## 4  4  212
## 5  5  224
## 6  6  379

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

## 'data.frame':    16 obs. of  2 variables:
##  $ id  : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ time: int  159 280 101 212 224 379 179 264 222 362 ...

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

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   101.0   169.5   223.0   241.5   268.0   485.0

# 使用qqnorm函数建立qq图
qqnorm(time.data$time)
qqline(time.data$time, col="red")

# 单样本t检验
# One-sample t-test
t.test(time.data$time, # 需要检验的数据
       alternative = c("two.sided"), # 双尾
       mu = 250, # 总体的平均值
       conf.int = 0.95 # 置信区间
       )

## 
##  One Sample t-test
## 
## data:  time.data$time
## t = -0.34439, df = 15, p-value = 0.7353
## alternative hypothesis: true mean is not equal to 250
## 95 percent confidence interval:
##  188.8927 294.1073
## sample estimates:
## mean of x 
##     241.5

# 使用boxplot函数建立盒装图
boxplot(time.data$time, ylab = "Time (hours)")
abline(h = 250, lty = 2, col ="red")

配对样本t检验 Paired t-test

练习题2：mCCP (meta-Chlorophenylpiperazine)被认为会影响人们食欲和人体对食物的吸收。因此，有些机构希望通过mCCP治疗肥胖症状。在一项有关mCCP减肥效果的研究中，9名适度肥胖的女性以双盲和安慰剂对照实验法给予mCCP。一些受试者先服mCCP两周，随后不服用任何药物(“冲洗时间”)，再继续服用安慰剂两周。另外一些受试者则相反。记录每一名受试者在每一种处理条件下减轻的体重(kg)。现在判断mCCP是否有效果？

# 加载数据
Input = ("
treatment id lost.weight
Control 1 1.1
Control 2 1.3
Control 3 1
Control 4 1.7
Control 5 1.4
Control 6 0.1
Control 7 0.5
Control 8 1.6
Control 9 -0.5
mCCP 1 0
mCCP 2 -0.3
mCCP 3 0.6
mCCP 4 0.3
mCCP 5 -0.7
mCCP 6 -0.2
mCCP 7 0.6
mCCP 8 0.9
mCCP 9 -2
")

mccp.data <- read.table(textConnection(Input), header = TRUE)

# 查看数据格式
head(mccp.data) # 展示数据前六行

##   treatment id lost.weight
## 1   Control  1         1.1
## 2   Control  2         1.3
## 3   Control  3         1.0
## 4   Control  4         1.7
## 5   Control  5         1.4
## 6   Control  6         0.1

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

## 'data.frame':    18 obs. of  3 variables:
##  $ treatment  : chr  "Control" "Control" "Control" "Control" ...
##  $ id         : int  1 2 3 4 5 6 7 8 9 1 ...
##  $ lost.weight: num  1.1 1.3 1 1.7 1.4 0.1 0.5 1.6 -0.5 0 ...

# 调整格式,将mccp.data$treatment改为分类变量(factor)
mccp.data$treatment <- as.factor(mccp.data$treatment)
levels(mccp.data$treatment)

## [1] "Control" "mCCP"

# mccp.data$lost.weight的数据汇总
summary(mccp.data$lost.weight)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -2.0000 -0.1500  0.5500  0.4111  1.0750  1.7000

# 使用qqnorm函数建立qq图
# 将预设的图分成两部分，第一部分用于装control.weight，第二部分用于装mccp.weight
par(mfcol=c(1,2))

# 使用which分别提取数据中的control.weight和mccp.weight部分
control.weight <- mccp.data$lost.weight[which(mccp.data$treatment == "Control")]
mccp.weight <- mccp.data$lost.weight[which(mccp.data$treatment == "mCCP")]

# 分别构建control.weight和mccp.weight的QQ图
qqnorm(control.weight, main = "QQ test of control.weight")
qqline(control.weight, col="red")

qqnorm(mccp.weight, main = "QQ test of mccp.weight")
qqline(mccp.weight, col="red")

# 配对样本t检验
# Paired t-test
t.test(control.weight, mccp.weight, #需要检验的数据，数据并排排列
       alternative = c("two.sided"), # 双尾
       paired = TRUE, # 确定t检验为配对处理
       conf.level = 0.95)

## 
##  Paired t-test
## 
## data:  control.weight and mccp.weight
## t = 4.1703, df = 8, p-value = 0.003121
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.4470396 1.5529604
## sample estimates:
## mean of the differences 
##                       1

# 使用boxplot函数建立盒装图
boxplot(lost.weight ~ treatment, data = mccp.data, ylab = "Lost weight (kg)", ylim = c(-2, 2))

独立样本t检验 Independent t-test

练习题3：表格展示两个小麦品种从播种到抽穗所需要的天数，试问两者所需的天数差异是否显著？

# 加载数据
Input = ("
treatment day
subsp01 101
subsp01 100
subsp01 99
subsp01 99
subsp01 98
subsp01 100
subsp01 98
subsp01 99
subsp01 99
subsp01 99
subsp02 100
subsp02 98
subsp02 100
subsp02 99
subsp02 98
subsp02 99
subsp02 98
subsp02 98
subsp02 99
subsp02 100
")

day.data <- read.table(textConnection(Input), header = TRUE)

# 查看数据格式
head(day.data) # 展示数据前六行

##   treatment day
## 1   subsp01 101
## 2   subsp01 100
## 3   subsp01  99
## 4   subsp01  99
## 5   subsp01  98
## 6   subsp01 100

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

## 'data.frame':    20 obs. of  2 variables:
##  $ treatment: chr  "subsp01" "subsp01" "subsp01" "subsp01" ...
##  $ day      : int  101 100 99 99 98 100 98 99 99 99 ...

# 调整格式,将day.data$treatment改为分类变量(factor)
day.data$treatment <- as.factor(day.data$treatment)
levels(day.data$treatment)

## [1] "subsp01" "subsp02"

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

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   98.00   98.00   99.00   99.05  100.00  101.00

# 使用qqnorm函数建立qq图
# 将预设的图分成两部分，第一部分用于装subsp01.day，第二部分用于装msubsp02.day
par(mfcol=c(1,2))

# 使用which分别提取数据中的subsp01.day和subsp02.day部分
subsp01.day <- day.data$day[which(day.data$treatment == "subsp01")]
subsp02.day <- day.data$day[which(day.data$treatment == "subsp02")]

# 分别构建subsp01.day和subsp02.day的QQ图
qqnorm(subsp01.day, main = "QQ test of subsp01.day")
qqline(subsp02.day, col="red")

qqnorm(subsp01.day, main = "QQ test of subsp02.day")
qqline(subsp02.day, col="red")

# 方差齐性检验
library(car)

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

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

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

leveneTest(day ~ treatment, #数据以y~x的形式输入
           data = day.data)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.1385 0.7142
##       18

# 独立样本t检验
# Independent t-test
# 如果方差齐性不满足，请使用var.equal = F代码
t.test(day ~ treatment, # 需要检验的数据，数据并排排列
       data = day.data, # 数据来源
       alternative = c("two.sided"), # 双尾 
       var.equal = T, #方差齐性满足
       conf.level = 0.95
       )

## 
##  Two Sample t-test
## 
## data:  day by treatment
## t = 0.74741, df = 18, p-value = 0.4645
## alternative hypothesis: true difference in means between group subsp01 and group subsp02 is not equal to 0
## 95 percent confidence interval:
##  -0.5432817  1.1432817
## sample estimates:
## mean in group subsp01 mean in group subsp02 
##                  99.2                  98.9

# 使用boxplot函数建立盒装图 
boxplot(day ~ treatment, data = day.data, ylab = "Day", ylim = c(98, 102))

订阅评论

2 评论

最旧

教学内容

参考资料

R语言教程

单样本t检验 One-sample t-test

配对样本t检验 Paired t-test

独立样本t检验 Independent t-test

你可能也喜欢

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

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

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