生物统计附试验设计

梅步俊 梁永厚 吴志红

目录

  • 1 绪论
    • 1.1 生物统计学发展史
    • 1.2 学习生物统计学的必要性
    • 1.3 常用术语和基本概念
    • 1.4 R软件的介绍
    • 1.5 进一步阅读的文献
    • 1.6 习题
  • 2 资料的描述性统计分析
    • 2.1 位置测度
    • 2.2 离散性测度
    • 2.3 R软件的应用
    • 2.4 习题
  • 3 随机变量与概率分布
    • 3.1 随机变量
    • 3.2 概率分布
    • 3.3 二维随机变量
    • 3.4 正态分布
    • 3.5 一些重要的概率分布
    • 3.6 R软件的应用
    • 3.7 习题
  • 4 参数估计方法
    • 4.1 估计量的评价准则
    • 4.2 点估计
    • 4.3 区间估计
    • 4.4 R软件的应用
    • 4.5 习题
  • 5 统计假设测验
    • 5.1 假设检验的基本问题
    • 5.2 统计检验的基本步骤
    • 5.3 抽样分布
    • 5.4 样本平均数与总体平均数差异显著性检验
    • 5.5 两个样本平均数的差异显著性检验
    • 5.6 百分数资料差异显著性检验
    • 5.7 总体参数的区间估计
    • 5.8 非参数检验
    • 5.9 R软件的应用
    • 5.10 习题
  • 6 方差分析
    • 6.1 单因素方差分析
    • 6.2 多重比较
    • 6.3 多因素方差分析
    • 6.4 方差分析需要满足的条件
    • 6.5 习题
  • 7 协方差分析
    • 7.1 协方差分析概述
    • 7.2 协方差分析的基本原理
    • 7.3 协方差分析的计算过程
    • 7.4 R软件的应用
    • 7.5 习题
  • 8 相关与回归分析
    • 8.1 变量之间的相互关系
    • 8.2 直线相关
    • 8.3 回归分析的性质
    • 8.4 一元正态线性回归统计模型
    • 8.5 多元线性回归统计模型
    • 8.6 自变量的选择与逐步回归
    • 8.7 曲线回归
    • 8.8 应用直线回归与相关的注意事项
    • 8.9 R软件的应用
    • 8.10 习题
  • 9 实验设计
    • 9.1 实验设计概述
    • 9.2 生物实验计划
    • 9.3 完全随机设计
    • 9.4 随机单位组设计
    • 9.5 拉丁方设计
    • 9.6 交叉设计
    • 9.7 正交设计
    • 9.8 R软件的应用
    • 9.9 习题
  • 10 附件
    • 10.1 复习题
R软件的介绍

第四节  R软件的介绍 

一、unix终端调用R

R [options] [<infile] [>outfile]或者R CMD BATCH infile

解释:没有在[]中指定任何内容,R将从终端打开,等待输入R命令。输出将被打印在屏幕上。要从R中退出,请键入q(),然后您就会出现问是否保存图像。如果选择保存,则创建的对象将保存在中文件“.RData”,到目前为止输入的所有命令都将保存在“.Rhistory”文件中。注意:这两个文件是隐藏文件。

如果给出“infile”,则R将在“infile”中执行R命令,并将输出写入“outfile”或者如果它是空的,则在屏幕上打印。您必须在[选项]中指定是否保存对象:--save,--no-save。

示范:

> a <- matrix(1:8,2,4)

> a

[,1] [,2] [,3] [,4]

[1,]  1  3   5  7

[2,]  2  4   6  8

> b <- matrix(1:8*0.1,2,4)

> b

[,1] [,2] [,3] [,4]

[1,]  0.1 0.3 0.50.7

[2,]  0.2 0.4 0.60.8

> a + b

[,1] [,2] [,3] [,4]

[1,] 1.1 3.3 5.5 7.7

[2,] 2.2 4.4 6.6 8.8

或者,可以在R控制台中键入以下命令以运行中的所有命令文件“demo-startup.R”:

source("demo-startup.R",echo=TRUE)

文件“demo-startup.Rout”中的上述内容将被打印在屏幕上。

R CMD BATCH file.R与。类似

R - save <file.R> file.Rout,但在“file.Rout”中写入屏幕上显示的所有消息,包括错误消息。

使用nohup在后台运行程序,例如:nohup R BATCH infile&

二、Windows调用 R

使用Windows命令行:将包含R.exe的目录添加到环境变量路径。您可以像使用Unix终端一样运行R。但nohup设施可能不存在。

点击桌面上的R启动,您需要首次更改工作目录。保存图像后(文件“.RData”将被创建),第二次单击文件“.RData”将从这个目录打开R。

三、获得帮助

从R控制台输入?关键字或help.search(关键字)

从互联网上:http://cran.r-project.org

四、对象和操作

1.数字和矢量

下面我使用示例来演示如何使用数字向量:

> a <- 1

> a

[1] 1

> b <- 2:10

> b

[1] 2 3 4 5 6 7 8 9 10

> #concatenate two vectors

> c <- c(a,b)

> c

[1] 1 2 3 4 5 6 7 8 9 10

> #vector arithmetics

> 1/c

[1] 1.0000000 0.5000000 0.3333333 0.2500000 0.20000000.1666667 0.1428571

[8] 0.1250000 0.1111111 0.1000000

> c^2

[1] 1 4 9 16 25 36 49 64 81 100

> c^2 + 1

[1] 2 5 10 17 26 37 50 65 82 101

> #apply a function to each element

> log(c)

[1] 0.0000000 0.6931472 1.0986123 1.3862944 1.60943791.7917595 1.9459101

[8] 2.0794415 2.1972246 2.3025851

> sapply(c,log)

[1] 0.0000000 0.6931472 1.0986123 1.3862944 1.60943791.7917595 1.9459101

[8] 2.0794415 2.1972246 2.3025851

> #operation on two vectors

> d <- (1:10)*10

> d

[1] 10 20 30 40 50 60 70 80 90 100

> c + d

[1] 11 22 33 44 55 66 77 88 99 110

> c * d

[1] 10 40 90 160 250 360 490 640 810 1000

> d ^ c

[1] 1.000000e+01 4.000000e+022.700000e+04 2.560000e+06 3.125000e+08

[6] 4.665600e+10 8.235430e+121.677722e+15 3.874205e+17 1.000000e+20

> #more concrete example: computing variance of ’c’

> sum((c - mean(c))^2)/(length(c)-1)

[1] 9.166667

> #of course, there is build-in function for computingvariance:

> var(c)

[1] 9.166667

> #subsetting vector

> c

[1] 1 2 3 4 5 6 7 8 9 10

> c[2]

[1] 2

> c[c(2,3)]

[1] 2 3

> c[c(3,2)]

[1] 3 2

> c[c > 5]

[1] 6 7 8 9 10

> #let’s see what is "c > 5"

> c > 5

[1] FALSE FALSE FALSE FALSE FALSE TRUE TRUE TRUE TRUETRUE

> c[c > 5 & c < 10]

[1] 6 7 8 9

> c[as.logical((c > 8) + (c < 3))]

[1] 1 2 9 10

> log(c)

[1] 0.0000000 0.6931472 1.0986123 1.3862944 1.60943791.7917595 1.9459101

[8] 2.0794415 2.1972246 2.3025851

> c[log(c) < 2]

[1] 1 2 3 4 5 6 7

> #modifying subset of vector

> c[log(c) < 2] <- 3

> c

[1] 3 3 3 3 3 3 3 8 9 10

> #extending and cutting vector

> length(c) <- 20

> c

[1] 3 3 3 3 3 3 3 8 9 10 NA NA NA NA NA NA NA NA NA NA

> c[25] <- 1

> c

[1] 3 3 3 3 3 3 3 8 9 10 NA NA NA NA NA NA NA NA NA NA NANA NA NA 1

> #getting back original vector

> length(c) <- 10

> c

[1] 3 3 3 3 3 3 3 8 9 10

> #introduce a function ‘‘seq’’

> seq(0,10,by=1)

[1] 0 1 2 3 4 5 6 7 8 9 10

> seq(0,10,length=20)

[1] 0.0000000 0.5263158 1.0526316 1.5789474 2.10526322.6315789

[7] 3.1578947 3.6842105 4.2105263 4.7368421 5.26315795.7894737

[13] 6.3157895 6.8421053 7.3684211 7.8947368 8.42105268.9473684

[19] 9.4736842 10.0000000

> 1:10

[1] 1 2 3 4 5 6 7 8 9 10

> #seq is more reliable than ":"

> n <- 0

> 1:n

[1] 1 0

> #seq(1,n,by=1)

> #Error in seq.default(1, n, by = 1) : wrong sign in’by’ argument

> #Execution halted

> #function ‘‘rep’’

> c<- 1:5

> c

[1] 1 2 3 4 5

> rep(c,5)

[1] 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5

> rep(c,each=5)

[1] 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 4 4 4 4 4 5 5 5 5 5

2.缺失值

> a <- 0/0

> a

[1] NaN

> is.nan(a)

[1] TRUE

> b <- log(0)

> b

[1] -Inf

> is.finite(b)

[1] FALSE

> c <- c(0:4,NA)

> c

[1] 0 1 2 3 4 NA

> is.na(c)

[1] FALSE FALSE FALSE FALSE FALSE TRUE

3.字符向量

字符串使用双引号(")或单引号(‘)输入,但是使用双引号打印(或有时不带引号)。他们使用C风格的转义序列,使用\作为转义字符,所以\\被输入和打印为\\,而在双引号内“作为\输入\。其他有用的转义序列是\ n,新行,\ t,制表符和\ b,退格。

> A <- c("a","b","c")

> A

[1] "a" "b" "c"

> paste("a","b",sep="")

[1] "ab"

> paste(A,c("d","e"))

[1] "a d" "b e" "c d"

> paste(A,10)

[1] "a 10" "b 10" "c 10"

> paste(A,10,sep="")

[1] "a10" "b10" "c10"

> paste(A,1:10,sep="")

[1] "a1" "b2" "c3""a4" "b5" "c6" "a7" "b8""c9" "a10"

4.矩阵

> A <- matrix(0,4,5)

> A

[,1] [,2] [,3] [,4] [,5]

[1,]  0  0  0  0  0

[2,]  0  0  0  0  0

[3,]  0  0  0  0  0

[4,]  0  0  0  0  0

> A <- matrix(1:20,4,5)

> A

[,1] [,2] [,3]  [,4]  [,5]

[1,]  1  5  9   13   17

[2,]  2  6  10  14  18

[3,]  3  7  11  15  19

[4,]  4  8  12  16  20

> #subsectioning and modifying subsection

> A[c(1,4),c(2,3)]

[,1]  [,2]

[1,]  5  9

[2,]  8  12

> A[c(1,4),c(2,3)] <- 1

> A

[,1] [,2] [,3]  [,4] [,5]

[1,]  1  1  1   13  17

[2,]  2  6  10  14  18

[3,]  3  7  11  15  19

[4,]  4  1  1   16  20

> A[4,]

[1] 4 1 1 16 20

> A[4,,drop = FALSE]

[,1] [,2] [,3] [,4] [,5]

[1,]  4   1  1  16  20

> #combining two matrices

> #create another matrix using another way

> A2 <- array(1:20,dim=c(4,5))

> A2

[,1] [,2] [,3] [,4] [,5]

[1,]  1  5   9  13 17

[2,]  2  6  10  14 18

[3,]  3  7  11  15 19

[4,]  4  8  12  16 20

> cbind(A,A2)

[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]

[1,]  1  1   1  13 17  1  5  9  13  17

[2,]  2  6  10  14 18  2  6  10  14  18

[3,]  3  7  11  15 19  3  7  11  15  19

[4,]  4  1   1  16 20  4  8  12  16  20

> rbind(A,A2)

[,1] [,2] [,3] [,4] [,5]

[1,]  1  1   1  13 17

[2,]  2  6  10  14 18

[3,]  3  7  11  15 19

[4,]  4  1   1  16 20

[5,]  1  5   9  13 17

[6,]  2  6  10  14 18

[7,]  3  7  11  15 19

[8,]  4  8  12  16 20

> #operating matrice

> #transpose matrix

> t(A)

[,1] [,2] [,3] [,4]

[1,]  1  2   3  4

[2,]  1  6   7  1

[3,]  1  10  11  1

[4,]  13  14  1516

[5,]  17  18  1920

> A

[,1] [,2] [,3] [,4] [,5]

[1,]  1  1   1  13 17

[2,]  2  6  10  14 18

[3,]  3  7  11  15 19

[4,]  4  1   1  16 20

> A + 1

[,1] [,2] [,3] [,4] [,5]

[1,]  2  2  2  14  18

[2,]  3  7  11  15  19

[3,]  4  8  12  16  20

[4,]  5  2  2   17  21

> x <- 1:4

> A*x

[,1] [,2] [,3] [,4] [,5]

[1,]  1  1   1  13  17

[2,]  4  12  20  28  36

[3,]  9  21  33  45  57

[4,]  16  4   4  64  80

> #the logical here is coercing the matrix"A" into a vector by joining the column

> #and repeat the shorter vector,x, as many times asmaking it have the same

> #length as the vector coerced from "A"

> #see another example

> x <- 1:3

> A*x

[,1] [,2] [,3]  [,4]  [,5]

[1,]  1  2   3   13  34

[2,]  4  18  10  28  54

[3,]  9  7   22  45  19

[4,]  4  2   3   16  40

> A^2

[,1] [,2] [,3]  [,4]  [,5]

[1,]  1  1   1  169  289

[2,]  4  36  100 196  324

[3,]  9  49  121 225  361

[4,]  16  1  1  256  400

> A <-matrix(sample(1:20),4,5)

> A

[,1] [,2] [,3] [,4] [,5]

[1,]  7  1   9  3  20

[2,]  11 10  12 15  4

[3,]  6  16  5  8  19

[4,]  2  17  14 13  18

> B <-matrix(sample(1:20),5,4)

> B

[,1] [,2]  [,3][,4]

[1,]  11  13  3  16

[2,]  2   4   7  10

[3,]  9   12  1  14

[4,]  5   15  17  8

[5,]  19  6   20  18

> C <- A %*% B

> C

[,1]  [,2] [,3][,4]

[1,]  555 368 488632

[2,]  400 576 450636

[3,]  544 436 651732

[4,]  589 565 720826

> solve(C)

[,1]          [,2]        [,3]        [,4]

[1,] 0.0113721434 -0.011305781 -0.03593289 0.03184764

[2,] 0.0037070714 -0.006309519 -0.03288573 0.03116506

[3,] -0.0008241436 -0.012304505 -0.02598824 0.03313549

[4,] -0.0099265186 0.023103180 0.07077051 -0.07169985

> #solving linear equation

> x <- 1:4

> d <- C %*% x

> solve(C,d)

[,1]

[1,] 1

[2,] 2

[3,] 3

[4,] 4

> #altenative way (but not recommended)

> solve(C) %*% d

[,1]

[1,] 1

[2,] 2

[3,] 3

[4,] 4

> #SVD (C = UDV’) and determinant

> svd.C <- svd(C)

> svd.C

$d

[1] 2332.552515 204.076932 98.799790 7.614026

$u

[,1]       [,2]       [,3]     [,4]

[1,] -0.4430091 0.41432949 0.7886463 0.1005537

[2,] -0.4432302 -0.84065092 0.2206557 -0.2194631

[3,] -0.5143024 0.34835579 -0.3848769 -0.6826500

[4,] -0.5854766 0.01689212 -0.4256969 0.6897202

$v

[,1]       [,2]       [,3]      [,4]

[1,] -0.4492025 0.45643331 0.66652400 0.3816170

[2,] -0.4172931 -0.83456029 0.09103668 0.3479769

[3,] -0.5024522 0.30793204 -0.73787822 0.3290219

[4,] -0.6096108 -0.01885767 0.05471573 -0.7905854

> #calculating determinant of C

> prod(svd.C$d)

[1] 358092916

5.列表

R列表是由称为其组件的对象的有序集合组成的对象。没有特别需要组件具有相同的模式或类型,例如,列表可以由一个数字向量、一个逻辑值、一个矩阵、一个复数向量、一个字符数组、一个函数等。

> a <- 1:10

> b <- matrix(1:10,2,5)

> c <- c("name1","name2")

> alst <- list(a=a,b=b,c=c)

> alst

$a

[1] 1 2 3 4 5 6 7 8 9 10

$b

[,1] [,2] [,3] [,4] [,5]

[1,]  1  3   5  7  9

[2,]  2  4   6  8  10

$c

[1] "name1" "name2"

> #refering to component of a list

> alst$a

[1] 1 2 3 4 5 6 7 8 9 10

> alst[[2]]

[,1] [,2] [,3] [,4] [,5]

[1,]  1  3   5  7  9

[2,]  2  4   6  8  10

> blst <- list(d=2:10*10)

> #concatenating list

> ablst <- c(alst,blst)

> ablst

$a

[1] 1 2 3 4 5 6 7 8 9 10

$b

[,1] [,2] [,3] [,4] [,5]

[1,]  1  3   5  7  9

[2,]  2  4   6  8  10

$c

[1] "name1" "name2"

$d

[1] 20 30 40 50 60 70 80 90 100

列表通常用于返回大型程序的结果,例如线性回归拟合的结果。

6.数据框架

数据框架可能出于多种目的可被视为具有可能具有不同模式和属性的列的矩阵。它可以以矩阵形式显示,其行和列使用矩阵索引约定提取。

> name <-c("john","peter","jennifer")

> gender <-factor(c("m","m","f"))

> hw1 <- c(60,60,80)

> hw2 <- c(40,50,30)

> grades <- data.frame(name,gender,hw1,hw2)

> grades

name gender hw1 hw2

1 john    m    60 40

2 peter    m    60  50

3 jennifer  f    80  30

> #subsectioning a data frame

> grades[1,2]

[1] m

Levels: f m

> grades[,"name"]

[1] john peter jennifer

Levels: jennifer john peter

> grades$name

[1] john peter jennifer

Levels: jennifer john peter

> grades[grades$gender=="m",]

name gender hw1 hw2

1 john  m    60  40

2 peter  m    60  50

> grades[,"hw1"]

[1] 60 60 80

> #divide the subjects by "gender", andcalculating means in each group

> tapply(grades[,"hw1"],grades[,"gender"],mean)

f m

80 60

7.从外部文件读取数据

文件“grades”显示如下:

name gender hw1 hw2

1 john    m  60  40

2 peter   m   60  50

3 jennifer f    80  30

创建grades数据框架:

> grades <- read.table("grades")

> grades

name gender hw1 hw2

1 john  m     60  40

2 peter  m     60  50

3 jennifer f      80 30

其他从不同格式的外部文件中读取数据的函数:read.csv,read.delim。它们是read.table的特定形式。有关更多详细信息,请键入?read.table。

8.编写你自己的函数

在文件“demo-fun.R”中,定义了一个函数:

#looking for the maximum value of a numeric vector x

find.max <- function(x)

{

n <- length(x)

x.m <- x[1]

ix.m <- 1

if(n > 1)

{

for( i in seq(2,n,by=1) )

{

if(x[i] > x.m)

{

x.m <- x[i]

ix.m <- i

}

}

}

#return the maximum value and the index

list(max=x.m,index.max=ix.m)

}

我们现在可以使用这个函数:

> #sourcing functions in file "demo-fun.R"

> source("demo-fun.R")

> x <- runif(12)

> x

[1] 0.06779191 0.09266746 0.63309784 0.859863120.81900862 0.80315468

[7] 0.71691262 0.35424083 0.59253821 0.23433190 0.608917870.35025762

> #calling "find.max"

> find.max(x)

$max

[1] 0.8598631

$index.max

[1] 4

9.用R制作图形

为了演示如何使用R生成图并在文件中保存图,我使用R绘制图来说明以下两个函数:

此处。R命令如下所示:

demofun1 <- function(x)

{

(1 + 2*x^2 + 3*x^3 + exp(x) ) / x^2

}

demofun2 <- function(x)

{

(1 + 2*(10-x)^2 + 3*(10-x)^3 + exp(10-x) ) / (10-x)^2

}

#open a file to draw in it

postscript("fig-latexdemo.eps",paper="special",

height=4.8, width=10, horizontal=FALSE)

#specify plotting parameters

par(mfrow=c(1,2), mar = c(4,4,3,1))

x <- seq(0,10,by=0.1)

#make "Plot 1"

plot(x, demofun1(x), type="p", pch = 1,

ylab="y", main="Plot 1")

#add another line to "Plot 1"

points(x, demofun2(x), type="l", lty = 1)

#make "plot 2"

plot(x, demofun1(x), type="b", pch = 3, lty=1 ,

ylab="y", main="Plot 2")

#add another line to "Plot 2"

points(x, demofun2(x), type="b", pch = 4, lty =2)

dev.off()

生成的文件包含结果图:

图1.1:该图演示了使用不同线和点的两个非线性函数

五、安装新的软件包并加载新的软件包

1.从源文件安装新的软件包

下载一个源文件,比如说apacakge.tar.gz。

运行R CMD INSTALL apackage.tar.gz -l /my/rlibrary

编译源文件并将其安装到目录/ my / library。

您必须安装GCC编译器,而不需要付出很多努力就可用于Windows系统。该包不需要从CRAN处获得。

2.从CRAN安装新软件包

如果软件包可从CRAN获得,则可以在不编译的情况下安装新软件包。预编译的软件包可以下载。为此,请输入install.packages("apackage",lib="/my/rlibrary")。

3.加载一个包

安装新的软件包apackage后,输入library("apackage",lib.loc="/my/library"),然后你可以使用这个包提供的功能。

4.调用C函数(仅适用于unix系统)

R运行很慢。我们最好使用C代码来进行密集计算。R向量的指针(矩阵将被强制为向量)传递给C代码,以便C程序可以使用它们。这是通过功能“.C”实现的。

下面是一个简单的演示。文件“sum.c”显示如下:

void newsum(int la[1], double a[], double s[1])

{

int i;

s[0] = 0;

for(i=0;i<la[0];i++)

s[0] += a[i];

}

然后编译上面的C程序:R CMD SHLIB sum.c

将生成两个文件,sum.o和sum.so。现在我们可以将“sum.so”加载到R环境中并使用.C调用C函数:

dyn.load("sum.so")

> a <- c(1,2,3,4)

> .C("newsum",length(a),a,sum=0)

[[1]]

[1] 4

[[2]]

[1] 1 2 3 4

$sum

[1] 10

通常,人们编写一个“包装器(wrapper)”函数,以便在R中定期调用C函数,例如:

newsum <- function(a)

{

.C("newsum",length(a),a,sum=0)$sum

}

六、创建您自己的软件包并提交给CRAN

它适用于高级用户。但是,这并不困难。按照如下程序:

1.创建一个目录“apackage”,它有3个目录:R(包含所有R源),src(包含所有C或Fortran文件),man(每个R函数的文档文件)和一个“DESCRIPTION”文件。编写文档文件和“描述”需要有格式。

2.(建议)运行R CMD check apackage以检查包装并进行建议更改。CRAN接受的软件包在此步骤中不得有任何警告和错误。

3.运行R CMD build apackage。

要提交给CRAN,只需将其(使用匿名ftp)发布到ftp://cran.r-project.org/incoming/。

有人会运行你的源文件,如果它通过R CMD check没有警告和错误,它会作为CRAN提供的软件包在线发布。Mac和Windows二进制文件将被创建。