第一章 引言

#字体用consolas，Editor theme使用Merbivore；
?boxplot  #求助方法
help("boxplot")  #求助方法
q()  #会询问是否保存工作空间
q(save = "no")  #不会询问是否保存工作空间

#将R代码保存到文本编辑器中，不要保存到工作空间；若计算机需要花费很长时间才能去完成的时候可以通过File——Save Workspace进行保存，通过

getwd()  #获取当前工作目录
rm(list=ls(all=TRUE)) 或者 rm(list=ls()) #移除所有变量

#更改R版本：tools--global Options--General--R version-change--choose a specific version of R语言初学者指南

# ctrl + shift + c 一次性注释很多行，一次性取消很多注释也是如此
#R包安装：install.packages("plotrix")；如果使用install.packages(plotrix)就无法安装成功
#小键盘上的0可以将R中的光标从竖直变为横着；

第二章 R 中的数据输入

Wingcrd<-c(59,55,53.5,55,52.5,57.5,53,55)  #c()函数生成了一个长度为8的向量
Wingcrd[1]  #查看第一个值
Wingcrd[1:5]  #查看前五个值
Wingcrd[-2]  #查看除第二个之外的值
sum(Wingcrd)  #求和
Tarsus<-c(22.3,19.7,20.8,20.3,20.8,21.5,20.6,21.5)
Head<-c(31.2,30.4,30.6,30.3,30.3,30.8,32.5,NA)
Wt<-c(9.5,13.8,14.8,15.2,15.5,15.6,15.6,15.7)
sum(Head)
sum(Head,na.rm = T) #na.rm 移除缺失值

#使用c,cbind,rbind结合变量
BirdData<-c(Wingcrd,Tarsus,Head,Wt)  #c结合变量会生成一个变量
BirdData
Z<-cbind(Wingcrd,Tarsus,Head,Wt)  #cbind结合的变量会以列的形式输出
Z
Z2<-rbind(Wingcrd,Tarsus,Head,Wt)  #rbind结合的变量会以行的形式输出
Z2

Id<-rep(c(1,2,3,4),each=8) 
Id

Id<-rep(1:4,each=8)
Id

a<-seq( from=1 , to = 4 , by=1 )
a
rep(a,each=8)

# ******Q1 ----------------------------------------------------------------------
#Q1: 如何产生 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4

VarName<-c("Wingcrd","Tarsus","Head","Wt")  #以字符形式结合

Id2<-rep(VarName,each=8)
Id2

Z<-cbind(Wingcrd,Tarsus,Head,Wt)  #cbind结合的变量会以列的形式输出,c表示column
Z

dim(Z)  #查看维度
dim(Z)[3]

使用矩阵结合数据

Dmat<-matrix(nrow = 8, ncol = 4)  #生成一个8X4的矩阵
Dmat
Dmat[,1]<-c(59,55,53.5,55,52.5,57.5,53,55)  #对矩阵进行数据填充
Dmat[,2]<-c(22.3,19.7,20.8,20.3,20.8,21.5,20.6,21.5)
Dmat[,3]<-c(31.2,30.4,30.6,30.3,30.3,30.8,32.5,NA)
Dmat[,4]<-c(9.5,13.8,14.8,15.2,15.5,15.6,15.6,15.7)
Dmat  #此时矩阵行和列还没有名字
colnames(Dmat)<-c("Wingcrd","Tarsus","Head","Wt")  #使用colnames函数给Dmat的列加上名字; 注意一点，就是这个函数放在了 <-的左边！
rownames(Dmat)<-c("甲","乙","丙","丁","戊","己","庚","辛")  #使用rownames函数给Dmat的行加上名字
Dmat

使用数据框（data.frame函数）结合数据

Dfrm<-data.frame(Wingcd=Wingcrd,Tarsus=Tarsus,Head=Head,Wt=Wt)  #
Dfrm

#cbind和matrix只能结合相同类型的数据,而data.frame可以结合不同类型的数据

使用list函数结合数据

#R中几乎所有的函数的输出结果都是保存在list中
x1 <- c(1,2,3)   #三个数字的向量
x2 <- c("a","b","c","d")   #三个字符的向量
x3 <- 3
x4 <- matrix(nrow = 2,ncol = 2)   #矩阵
x4[,1] <- c(1,2)  #给矩阵赋值
x4[,2] <- c(3,4)  #给矩阵赋值
x4
y <- list(x1=x1,x2=x2,x3=x3,x4=x4)  #用list将x1,x2,x3,x4这四个变量结合起来,若改为y <- list(x1,x2,x3,x4)，则下一步的y$x2无法正常运行. =左边是自定义的名字，=右边是实质性内容；
y$x2
y
y[[1]] #取list的元素使用[[]]来实现
y[[2]]
y[[2]][1]

AllData <- list(BirdData=BirdData,Id=Id2,Z=Z,VarName=VarName)  # =左边是自定义的名字，=右边是实质性内容；
AllData
AllData$BirdData
AllData[[1]]
tmp<-AllData[[3]]
class(tmp)
AllData[[3]][1,] #因为AllData第三个元素Z是matrix，所以AllData[[3]][1,] 表示取该matrix的第一行
AllData[[3]][,1] #表示取该matrix的第一列

# list 和 data.frame 的差别？
# ******数据载入 -----------------------------------------------------------------
# read.table函数
# read.table读取的是txt文件

#可以使用绝对路径
Squid <- read.table(file = "D:/生信学习/猪年年尾R学习笔记/R语言初学者指南/乌贼数据.txt",header =T,dec=".")  #file =不写也可以；read.table读取的是txt文件！

#也可以先建立.Rproj后缀文件，然后打开该文件进行定位，然后打开Rstudio文件。注意：要在文件所在位置建立Rproj文件，但是要打开的Rstudio文件位置在哪里都可以，也就是说Rstudio文件并不需要和Rproj在同一个文件夹下。
Squid <- read.table(file = "乌贼数据.txt",header =T,dec=".") 
#header =T 有表头；dec=".",小数点字符为.
Squid
head(Squid)  #head函数看前几行
# 当一个目录下面有很多文件需要读取，可以先利用setwd函数设置一下工作路径，然后再read.table
setwd("D:/生信学习/猪年年尾R学习笔记/R语言初学者指南/")  #若上面使用了.Rproj进行了定位，则这里不可以更改路径。
Squid<-read.table(file="乌贼数据.txt",header = T)

# ******Answer of Q1 ---------------------------------------------------------
# 把each换成times就可以
Id<-rep(c(1,2,3,4),each=8) 
Id

Id<-rep(c(1,2,3,4), times=8) 
Id

第三章 访问变量和处理数据子集

#read.table函数应该和names以及str函数始终结合起来；因为这样可以查看变量，并且查看每个变量的属性；保证下一步不会出现错误；
setwd("D:/生信学习/猪年年尾R学习笔记/R语言初学者指南")
Squid<-read.table(file = "乌贼数据.txt",header = T,dec = ".")
Squid
names(Squid) #查看正在处理的变量
str(Squid) #查看数据框中每个变量的属性
head(Squid) #查看前几行
class(Squid)

Squid2<-read.table(file = "乌贼数据.txt",header = T,dec = ",") #我们如果将小数点错误的标为“，” 则GSI会被认为成一个因子(它不再是数值型变量)，继续使用函数作图将会产生错误；
names(Squid2) #查看正在处理的变量
str(Squid2) #查看数据框中每个变量的属性。发现出了问题，GSI成为了因子。
class(Squid2)

访问数据框中的变量

方法一：$符号进行访问

Squid$GSI  #通过$符号访问数据框中的变量；
mean(Squid$GSI)  #计算均值

方法二：先attach，然后输入GSI

attach(Squid)  #使用attach命令将Squid添加到R的搜索路径里面
GSI  #避免每次输入Squid$
boxplot(GSI)
detach(Squid)  #用detach命令将Squid从R的搜索路径里面移除
GSI  #发现找不到

结合布尔运算符访问数据子集

unique(Squid$Sex)  #显示这个变量里面有多少个唯一值

访问所有性别为1的数据

Sel <- Squid$Sex ==1  #生成一个布尔向量用来选择行，若sex为1则该变量的值是TRUE，否则为FALSE
Sel
SquidM <- Squid[Sel,]  #选择Sel为TRUE的行; 注意选的是“行”，所以Sel放 = 前面
#以上两个步骤也可以合成一个步骤：SquidM <- Squid[Squid$Sex ==1,]
SquidM

访问数据Squid的Location为1,2,3的子集(即Location为1，2，3的行，所以后面要加逗号)

Squid123 <- Squid[Squid$Location ==1|Squid$Location ==2|Squid$Location ==3,] #|表示布尔运算符“或”
#以下三种也可以
# Squid123 <- Squid[Squid$Location !=4,] != 表示“不等于”
# Squid123 <- Squid[Squid$Location <=3,]
# Squid123 <- Squid[Squid$Location >=1 & Squid$Location <=3,]
Squid.1 <- Squid[Squid$Sex ==1 & Squid$Location ==1,] #需要满足Sex为1和Location为1两个条件，因为要提取符合条件的行，所以写在 = 前面
Squid.1
Squid.2 <- Squid[Squid$Sex ==1 & (Squid$Location ==1|Squid$Location ==2),] #需要满足Sex为1和Location为1或者2两个条件
Squid.2

数据排序(order)

Squid
Ord1 <- order(Squid$MONTH) #对月份进行排序
Ord1  
Squid[Ord1,]  #月份排序之后的Squid，为何要放在行上？？？因为要以排序后的月份对Squid每一行进行排列！
Squid$MONTH[Ord1] #月份排序之后的月份
Squid$Sample[Ord1]  #月份排序之后的Sample

数据集组合(merge)

setwd("D:/生信学习/猪年年尾R学习笔记/R语言初学者指南")
Sq1 <- read.table(file = "乌贼数据1.txt",header = T) #读取文件1
Sq2 <- read.table(file = "乌贼数据2.txt",header = T) #读取文件2
Sq1
Sq2
SquidMerged <- merge(Sq1, Sq2, by="Sample",all=T)  #使用变量Sample作为相同的标识符合并两个数据集,all=T表示如果有缺失值将用NA来填充；如果all=F则表示如果有缺失值将忽略；
#所谓的缺失值是指整行的缺失(这一点非常重要，意思是Sq1缺失了几个整行或者Sq2缺失了几个整行！)
SquidMerged

数据输出(write.table)

SquidM <- Squid[Squid$Sex ==1,]  #提取Sex为1的数据
write.table(SquidM,file = "MaleSquid.txt",sep = "",quote = F,append = F, na = "NA") #将上一步提取出来的数据输出到MaleSquid.txt里面
getwd()  #当前路径就是该文件所在路径
#通过write.table输出之后就可以发现多了一个文件，名为MaleSquid.txt。

重新编码分类变量(factor)

Squid
Squid$fLocation <- factor(Squid$Location, levels = c(1,2,3,4), labels = c("甲","乙","丙","丁")) #创建一个新的列，名字叫fLocation，它的内容就是将Location的1,2,3,4换为了甲乙丙丁
Squid$fLocation
Squid
Squid$fSex <- factor(Squid$Sex,levels = c(1,2),labels = c("M","F")) #创建一个新的列，名字叫fLocation，它的内容就是将Sex的1,2换为了M和F
Squid$fSex
Squid
boxplot(GSI ~ fSex, data =Squid) #把Squid的GSI和fSex分别拿出来，对应着进行boxplot画图；

第四章 简单的函数

setwd("c:/RBook/")
Veg <- read.table(file="Vegetation2.txt", header =TRUE)
Veg
head(Veg)
names(Veg)
str(Veg)

m <- mean(Veg$R)  #Veg$R表示将Veg中的R这一子集单独提取出来；mean(Veg$R)表示对Veg中R这一子集进行求取平均值；
m
m1 <- mean(Veg$R[Veg$Transect == 1])  #从外往里看，首先这是求平均值，求Veg中R这一子集的平均值，但不是求取所有R的平均值，有一个条件，就是被求取平均值的R要满足Transect为1, []之中为R要满足的条件。注意条件要用[ ]括起来；
m2 <- mean(Veg$R[Veg$Transect == 2])
m3 <- mean(Veg$R[Veg$Transect == 3])
m4 <- mean(Veg$R[Veg$Transect == 4])
m5 <- mean(Veg$R[Veg$Transect == 5])
m6 <- mean(Veg$R[Veg$Transect == 6])
m7 <- mean(Veg$R[Veg$Transect == 7])
m8 <- mean(Veg$R[Veg$Transect == 8])
c(m, m1, m2, m3, m4, m4, m5, m6, m7, m8)

tapply函数

tapply函数可以根据INDEX的不同水平对X进行FUN的运算；tapply(y,x,FUN=mean)

Veg$Transect
Me <- tapply(Veg$R, Veg$Transect, mean)  #等同于 Me <- tapply(X = Veg$R, INDEX = Veg$Transect ,FUN = mean)
Me
Sd <- tapply(Veg$R, Veg$Transect, sd)   #等同于 Sd <- tapply(X = Veg$R, INDEX = Veg$Transect, sd)
Sd
Le <- tapply(Veg$R, Veg$Transect, length) #等同于 Le <- tapply(X = Veg$R, INDEX = Veg$Transect, length)
Le
cbind(Me, Sd, Le)

sapply函数

sapply函数对y的每一个变量使用FUN的函数；sapply(y,FUN=mean)

Me <- mean(Veg$R) #求整个序列的均值
Mi <- min(Veg$R)  #求整个序列的最小值
Ma <- max(Veg$R)  #求整个序列的最大值
Sd <- sd(Veg$R)  #求整个序列的标准差
Le <- length(Veg$R)  #求整个序列的长度
c(Me,Mi,Ma,Sd,Le)
names(Veg)  #通过names函数我们可以发现有很多变量
length(Veg)  #通过length函数我们发现有24个变量，如果我们想对每一个变量都求平均值，就需要使用mean()函数24次
#sapply函数可以避免以上囧状

sapply(Veg[, 5:10], FUN = mean)  #对Veg的5到10列求取mean值
sapply(Veg[, 5:10], function(x) {mean(x)})  #与上面等价

lapply(Veg[, 5:10], FUN = mean)  #lapply和sapply函数一样功能，只不过输出结果形式不一样；sapply函数的输出是一个向量，lapply函数的输出是一个列表；
lapply(Veg[, 5:10], function(x) {mean(x)})

总结：tapply函数计算的是一个变量观察值子集的均值(或其他函数)；
lapply和sapply函数计算的是一个或多个变量全部观察值的均值（或其他函数）；
lapply和sapply中包含数据的变量必须是数据框

summary函数

summary函数可以提供变量基本信息

Z <-cbind(Veg$R, Veg$ROCK, Veg$LITTER, Veg$ML)  #使用cbind将几个变量结合起来
Z
head(Z)
colnames(Z) <- c("R","ROCK","LITTER","ML")  #如果没有这一步输出的结果每一列会没有名字，这一步的作用是给每一列添加名字！
head(Z)
summary(Z) #提供变量信息
names(Veg)

summary(Veg[,c("R","ROCK","LITTER","ML")])  #可以实现同样的功能
summary(Veg[,c(5,6,7,8)])  #可以实现同样的功能

setwd("D:/生信学习/猪年年尾R学习笔记/R语言初学者指南/RBook")
Deer <- read.table(file="Deer.txt")
names(Deer)
str(Deer)

table(Deer$Farm)
table(Deer$Sex, Deer$Year)

第五章 基础绘图工具简介

setwd("c:/RBook/")
Veg <- read.table(file="Vegetation2.txt",header =TRUE)
plot(Veg$BARESOIL, Veg$R, xlab = "BARESOIL", ylab = "R")  #注意自变量和因变量，也可以 plot(x = Veg$BARESOIL, y = Veg$R, xlab = "BARESOIL", ylab = "R")来指定自变量和因变量
plot(x = Veg$BARESOIL, y = Veg$R, xlab = "横坐标", ylab = "纵坐标",
     xlim = c(0,50), ylim =c(0,20),main="标题", pch=16)  #xlim和ylim设定取值范围,pch是指定绘图字符；
?plot  #可以查看很多参数

不同时间使用不同符号来表示

unique(Veg$Time)  #查看总共有几个年份
Veg$Time2 <- Veg$Time  #将这些年份生成一个对于pch有效的新的数值向量
Veg$Time2[Veg$Time == 1958] <-1 #[]里面是条件，此处也可以改为Veg$Time2[Veg$Time2 == 1958] <-1，因为在<-左边Time和Time2是对应起来的，<-右边才是条件判断过之后进行的赋值；
Veg$Time2[Veg$Time == 1962] <-2
Veg$Time2[Veg$Time == 1967] <-3
Veg$Time2[Veg$Time == 1974] <-4
Veg$Time2[Veg$Time == 1981] <-5
Veg$Time2[Veg$Time == 1994] <-6
Veg$Time2[Veg$Time == 2002] <-7
Veg$Time2[Veg$Time == 1989] <-8
plot(x = Veg$BARESOIL, y = Veg$R, xlab = "横坐标", ylab = "纵坐标",
     xlim = c(0,50), ylim =c(0,20),main="标题", pch=Veg$Time2)  #我们将不同年份的数据使用不同的符号进行了标示！

#也可以如下操作
Veg$Time3 <- Veg$Time
Veg$Time3
Veg$Time3[Veg$Time <= 1974] <- 1  ##[]里面是条件，此处也可以改为Veg$Time3[Veg$Time3 <= 1974] <- 1，因为在<-左边Time和Time3是对应起来的，<-右边才是条件判断过之后进行的赋值；
Veg$Time3[Veg$Time > 1974] <- 16
plot(x = Veg$BARESOIL, y = Veg$R, xlab = "横坐标", ylab = "纵坐标",
     xlim = c(0,50), ylim =c(0,20),main="标题", pch=Veg$Time3)

改变绘图符号颜色

#在plot函数里面对col选项使用向量
Veg$Time3 <- Veg$Time
Veg$Time3
Veg$Time3[Veg$Time <= 1974] <- 1  ##[]里面是条件，此处也可以改为Veg$Time3[Veg$Time3 <= 1974] <- 1，因为在<-左边Time和Time3是对应起来的，<-右边才是条件判断过之后进行的赋值；
Veg$Time3[Veg$Time > 1974] <- 16

Veg$col3[Veg$Time <= 1974] <- 1  #生成颜色向量1
Veg$col3[Veg$Time > 1974] <- 2   #生成颜色向量2
plot(x = Veg$BARESOIL, y = Veg$R, xlab = "横坐标", ylab = "纵坐标",
     xlim = c(0,50), ylim =c(0,20),main="标题", pch=Veg$Time3, col=Veg$col3)

改变绘图符号尺寸

# 使用cex选项改变尺寸
plot(x = Veg$BARESOIL, y = Veg$R, xlab = "横坐标", ylab = "纵坐标",
     xlim = c(0,50), ylim =c(0,20),main="标题", pch=Veg$Time3, col=Veg$col3, cex=5)

#也可以对cex使用向量，即对不同年份使用不同尺寸来表示
Veg$cex2 <- Veg$Time
unique(Veg$Time)
Veg$cex2[Veg$Time <= 1974] <- 1.5
Veg$cex2[Veg$Time > 1974] <- 4
plot(x = Veg$BARESOIL, y = Veg$R, xlab = "横坐标", ylab = "纵坐标",
     xlim = c(0,50), ylim =c(0,20),main="标题", pch=Veg$Time3, col=Veg$col3, cex=Veg$cex2)

添加一条平滑线

plot(x = Veg$BARESOIL, y = Veg$R, xlab = "横坐标", ylab = "纵坐标",
     xlim = c(0,50), ylim =c(0,20),main="标题", pch=Veg$Time3, col=Veg$col3, cex=Veg$cex2)
##不进行排序直接拟合会画出多条线
M.Loess <- loess(R ~ BARESOIL, data = Veg)  #R ~ BARESOIL表示R可以作为BARESOIL的函数
Fit <- fitted(M.Loess)  #提取拟合值并将其赋值给Fit
lines(Veg$BARESOIL, Fit)  #第一个参数绘制横坐标，第二个参数绘制纵坐标，画线,但是发现是多条线,因为line命令的第一个参数是按照表里面顺序连接点的

#所以我们要先进行排序，然后才可以进行线的拟合
plot(x = Veg$BARESOIL, y = Veg$R, xlab = "横坐标", ylab = "纵坐标", xlim = c(0,50), ylim =c(0,20),main="标题", pch=Veg$Time3, col=Veg$col3, cex=Veg$cex2)
M.Loess <- loess(R ~ BARESOIL, data = Veg)  #R ~ BARESOIL表示R可以作为BARESOIL的函数
Fit <- fitted(M.Loess)  #提取拟合值并将其赋值给Fit
Ord1 <- order(Veg$BARESOIL)
lines(Veg$BARESOIL[Ord1], Fit[Ord1], lwd = 5, lty = 2, col="black")  #lwd表示线的宽度，lty表示线的类型,col可以用“颜色单词”表示，也可以直接用数字表示！

第六章 循环与函数

补充一个知识点

sort、order函数的分别

sort()函数是对向量进行从小到大的排序
order()函数返回的值表示位置

举个例子

a <- c(6, 8, 2, 5)
sort(a)
order(a) #order()函数返回的值表示位置，返回值是3，4，1，2表示第三个数排在第一位，第四个数排在第二位，第一个数排在第三位，第二个数排在第四位。


setwd("C:/RBook/")
Owls <- read.table(file="Owls.txt",header = TRUE)
names(Owls)
str(Owls)
head(Owls)
dim(Owls)
unique(Owls$Nest)

#想要探索Owls文件中每一个Nest子集里面Arrival Time和Negotiation behaviour的关系
#例子1：探索Owls文件中第一个Nest子集(AutavauxTV)里面Arrival Time和Negotiation behaviour的关系
Owls.ATV<-Owls[Owls$Nest=="AutavauxTV",]  #提取子集！这个非常的好理解，可以打开Owls.txt对应着看一下；[]里面是条件，Owls$Nest表示文件中的第一列，
# 但是通过unique(Owls$Nest)函数也可以发现第一列有好几个不同的行，所以=="AutavauxTV",表示要名字为AutavauxTV的这几行！
plot(x = Owls.ATV$ArrivalTime, y = Owls.ATV$NegPerChick,
     xlab = "Arrival Time", ylab = "Negotiation behaviour",
     main =  "AutavauxTV", cex=1.2, col = "red", pch = 16)

#例子2：探索Owls文件中第二个Nest子集(Bochet)里面Arrival Time和Negotiation behaviour的关系
Owls.Bot <- Owls[Owls$Nest=="Bochet",]  #提取Bochest子集
plot(x = Owls.Bot$ArrivalTime, y = Owls.Bot$NegPerChick,
     xlab = "Arrival Time", ylab = "Negotiation behaviour",
     main =  "Bochet", cex=1.2, col = "black", pch = 16)

#例子3：探索Owls文件中第三个Nest子集(Champmartin)里面Arrival Time和Negotiation behaviour的关系
Owls.Cha <- Owls[Owls$Nest=="Champmartin",]  #提取Bochest子集
plot(x = Owls.Cha$ArrivalTime, y = Owls.Cha$NegPerChick,
     xlab = "Arrival Time", ylab = "Negotiation behaviour",
     main =  "Champmartin", cex=1.2, col = "purple", pch = 16)
#                              .
#                              .
#                              .
#                              
#                 这样下去需要重复25遍才能结束！

# Nest.i <- "子集名称"
# Owls.i <- Owls[Owls$Nest==Nest.i,]
# plot(x = Owls.i$ArrivalTime, y = Owls.i$NegPerChick,
#      xlab = "Arrival Time", ylab = "Negotiation behaviour",
#      main =  Nest.i, cex=1.2, col = "purple", pch = 16)

图像的保存

#将上面的第一个保存下来
setwd("D:/生信学习/猪年年尾R学习笔记/R语言初学者指南/AllGraphs")  #设定好图片保存的路径, 有时候会报错，重新将地址复制过来就好了。
jpeg(filename = "AutavauxTV.jpg")  #打开画布，需要先打开画布，然后再plot才可以！
Owls.ATV<-Owls[Owls$Nest=="AutavauxTV",]  #提取子集
plot(x = Owls.ATV$ArrivalTime, y = Owls.ATV$NegPerChick,
     xlab = "Arrival Time", ylab = "Negotiation behaviour",
     main =  "AutavauxTV", cex=1.2, col = "red", pch = 16)  #作图
dev.off()  #关闭画布，不关闭则无法打开图片

构造循环

AllNests <- unique(Owls$Nest)  #使用unique确定Nest名字唯一性；
for (i in 1:27){
  Nest.i <- AllNests[i]  #使Nest.i 和第i个Nest名字对应起来
  Owls.i <- Owls[Owls$Nest == Nest.i,]  #提取第i个子集
  YourFileName <- paste(Nest.i,".jpg",sep="")  #paste可以将变量连接为字符串；使用paste函数对每一个图片的进行命名，可以分别改变一下里面参数，看一下什么变化。
  jpeg(file=YourFileName)  #打开画布，这一步要在plot之前
  plot(x = Owls.i$ArrivalTime, y = Owls.i$NegPerChick,
       xlab = "Arrival Time", ylab = "Negotiation time",
       main =  Nest.i, pch=16)  #画图
  dev.off()  #关闭画布
}

函数

对每一个变量制作一个表格给出缺失值和零的个数

例子1：确定每个变量中缺失值的数量

setwd("C:/RBook")
Veg <- read.table(file = "Vegetation2.txt", header = T)
names(Veg)
str(Veg)
head(Veg)

NAPerVariable <- function(X1){  #function构造一个函数，函数名字是 NAPerVariable
  D1 <- is.na(X1)  #这里进行一个判断； is.na(X1)生成一个与X1维数相同的布尔矩阵，并将X1矩阵中的缺失值计为TRUE，非缺失值计为FALSE；
  colSums(D1)  #计算每一列元素的和；colSums是R中已有的函数，colsum函数的对象通常为数值矩阵，也可以是布尔矩阵，它会将布尔矩阵中的TRUE计为1，FALSE计为0；并且变量的任何位置出现了空值则colSum函数的输出就为空值，可以通过设定na.rm = TRUE来实现；
}

#X1和D1的解释
# X1是NAPerVariable这个函数的变量，把什么变量给NAPerVariable函数，则X1就是什么
# D1是通过is.na(X1)生成的与X1维数相同的布尔矩阵，并将X1矩阵中的缺失值计为TRUE，非缺失值计为FALSE；

NAPerVariable(Veg[,1:24]) #因为NAPerVariable是刚刚写好的函数，所以直接拿来用，在这里X1就是Veg[,1:24]

例子2：确定每个变量中缺失值的数量

setwd("c:/RBook/")
Parasite <- read.table(file="CodParasite.txt", header = TRUE)
names(Parasite)
str(Parasite)

NAPerVariable<-function(X1) {
  D <- is.na(X1)
  colSums(D)
}

NAPerVariable(Parasite) #在这里X1就是Parasite

例子3：确定每个变量中零值的数量

ZerosPerVariable<-function(X1) {
  D1 =(X1 == 0)  #在这里进行判断，等于0就是TRUE，否则就是FALSE；这个括号有没有都可以的；
  colSums(D1, na.rm = T)  #通过设定na.rm = TRUE来移除缺失值
}

ZerosPerVariable(Parasite) #在这里X1就是Parasite

例子4：同时确定每个变量中缺失值和零值的数量

VariableInfo<-function(X1,Choice1) {  #在这里Choice1是一个有两种选择的变量。
  if (Choice1 =="Zeros"){ D1 = X1 == 0 }
  if (Choice1 =="NAs")  { D1 <- is.na(X1)}
  colSums(D1, na.rm = TRUE)
}
#函数 VariableInfo 有两个参数：X1 和 Choice1。X1是数据框，Choice1是一个有两种选择的变量。

VariableInfo(Parasite,"Zeros")  #在这里函数 VariableInfo的两个参数：X1 和 Choice1分别对应着 Parasite 和 "Zeros"
VariableInfo(Parasite,"NAs")
VariableInfo(Parasite,"error test") #因为Choice1只有两种可能性的选择，"Zeros"或者"NAs"，"error test"不符合两种选择中的任何一个，所以会显示error

例子5：同时确定每个变量中缺失值和零值的数量，并且当输入错误时可以显示特定的提示信息

VariableInfo<-function(X1,Choice1 = "Zeros") {  #这里写的Choice1 = "Zeros"而不仅仅写的Choice1，就相当于设置了默认值为Choice1 = "Zeros"
  if (Choice1 =="Zeros"){ D1 <- X1 == 0 }
  if (Choice1 =="NAs")  { D1 <- is.na(X1)}
  if (Choice1 != "Zeros" & Choice1 != "NAs") {print("What the fuck! You made a typo! You're so stupid!!!")} else {
    colSums(D1, na.rm = TRUE)}
}

VariableInfo(Parasite)  #因为默认值是Choice1 = "Zeros"，所以这里不给choice1的值会按Choice1 = "Zeros"来进行计算。
VariableInfo(Parasite,"Zeros")
VariableInfo(Parasite,"NAs")
VariableInfo(Parasite,"error test")  #若输入错误则会有特定的提示信息！

例子6：对例子5进行简化，使if数量达到最少

VariableInfo<-function(X1, Choice1 = "Zeros") {
  ifelse (Choice1 =="Zeros", D1 <- (X1 == 0) , D1 <- is.na(X1))  # ifelse替代了前面两个if命令
  if (Choice1 != "Zeros" & Choice1 != "NAs") {print("What the fuck! You made a typo! You're so stupid!!!")} else {
    colSums(D1, na.rm = TRUE)}}

VariableInfo(Parasite)
VariableInfo(Parasite,"NAs")
VariableInfo(Parasite,"Zeros")
VariableInfo(Parasite,"error test")  #若输入错误则会有特定的提示信息！


setwd("C:/RBook/")
Benthic <- read.table("RIKZ.txt",header=T)
Species <- Benthic[,2:76]
n <- dim(Species)
n
sum(Species[1, ], na.rm = TRUE)
sum(Species[2, ], na.rm = TRUE)
sum(Species[3, ], na.rm = TRUE)
sum(Species[4, ], na.rm = TRUE)
#        .
#        .
#        .
# 构造循环来避免该命令写45遍
TA <- vector(length = n[1])  #构造了一个长度为45的空向量
TA  #此时里面45个位置都是空的
for (i in 1:n[1]){
  TA[i] <- sum(Species[i,], na.rm = TRUE)
}

TA

# 其实有更简洁的方式
TA <- rowSums(Species, na.rm = T)  # rowSums函数计算每一行的和
TA  

统计每一行的数值大于0的个数

dim(Species) #总共45行75列
sum(Species[1, ] > 0, na.rm = TRUE)  #Species[1, ] > 0是一个判断式，起到了判断的作用，对Species第一行的75个值进行判断，若该值大于0则为TRUE，返回1，否则为FALSE返回0。然后sum进行了求和，得到的值就是TRUE的个数，即第一行75个数值中满足 >0 这一条件的个数。
sum(Species[2, ] > 0, na.rm = TRUE)
sum(Species[3, ] > 0, na.rm = TRUE)
sum(Species[4, ] > 0, na.rm = TRUE)
#        .
#        .
#        .
# 构造循环来避免该命令写45遍

Richness <- vector(length = n[1]) #构造了一个长度为45的空向量
Richness  #此时里面45个位置都是空的
for (i in 1:n[1]){
  Richness[i] <- sum(Species[i, ] > 0, na.rm = TRUE)
}  #for循环对这45个空位置进行一一填充
Richness

# 其实有更简洁的方式
Richness <- rowSums(Species > 0, na.rm = TRUE) # Species > 0对每一行数据进行判断，若判断成立则返回TRUE，值为1，否则返回FALSE值为0，rowSums函数在这里计算的结果就是满足条件的值的个数。
Richness

RS <- rowSums(Species, na.rm = TRUE)  #因为这里不是判断式，所计算的是每一行的和
RS

prop <- Species / RS
prop
H <- -rowSums(prop * log10(prop), na.rm = TRUE)
H

suppressMessages(install.packages("vegan"))
library(vegan)
H <- diversity(Species)
H

Choice <- "Richness"

if (Choice == "Richness") {
  Index <- rowSums(Species >0 , na.rm = TRUE)}
if (Choice == "Total Abundance") {
  Index <- rowSums(Species, na.rm = TRUE) }
if (Choice=="Shannon") {
  RS <- rowSums(Species, na.rm = TRUE)
  prop <- Species / RS
  Index <- -rowSums(prop * log10(prop), na.rm = TRUE)}


Index.function <- function(Spec, Choice){
  if (Choice == "Richness") {
    Index <- rowSums(Spec>0, na.rm = TRUE)}
  if (Choice == "Total Abundance") {
    Index <- rowSums(Spec, na.rm = TRUE) }
  if (Choice=="Shannon") {
    RS <- rowSums(Species, na.rm = TRUE)
    prop <- Species / RS
    Index <- -rowSums(prop * log10(prop), na.rm = TRUE)}
  list(Index = Index, MyChoice = Choice)
}

第七章 图形工具

setwd("C:/RBook/")
BFCases <- read.table(file="BirdFluCases.txt", header = TRUE,sep="\t")
names(BFCases)
dim(BFCases)
str(BFCases)
head(BFCases)
Cases <- rowSums(BFCases[,2:16])  # rowSum对每一行进行求和计算；对BFCases的2到16列的每一行求和！
Cases  #此时还没有名字
names(Cases) <- BFCases[,1]  # names函数将BFCases的第一列内容赋给了Cases，作为Cases的名字；
Cases  #此时有了名字
head(BFCases)

pie函数画饼图

Piechart

par(mfrow = c(2,2), mar = c(3, 3, 2, 1)) #mfrow是指面板能生成的图形个数与排版；mar是指每个图形周围空间大小；(改一下参数试下就明白了！)
pie(Cases , main = "Ordinary pie chart")
pie(Cases , col = gray(seq(0.4,1.0,length=6)), clockwise=T, main = "Grey colours")  #clockwise=T是指显示的方向
pie(Cases , col = rainbow(6), clockwise = TRUE, main="Rainbow colours")

3D Exploded Pie Chart

library(plotrix)
pie3D(Cases , labels = names(Cases), explode = 0.1, main = "3D pie chart", labelcex=0.6)

par函数对画布进行设置

op <- par(mfrow = c(2,2), mar = c(3, 3, 2, 1))  # mfrow是指面板能生成的图形个数与排版；mar是指每个图形周围空间大小；(改一下参数试下就明白了！)
pie(Cases , main = "Ordinary pie chart")
pie(Cases , col = gray(seq(0.4,1.0,length=6)), clockwise=TRUE, main = "Grey colours")
pie(Cases , col = rainbow(6),clockwise = TRUE, main = "Rainbow colours")
pie3D(Cases , labels = names(Cases ), explode = 0.1, main = "3D pie chart", labelcex=0.6)
par(op)  #图行参数设置保存在了第一行代码的op里面，这里再进行一下par(op)表示结束par设置，就像par从来没有使用过一样。

绘制条形图

BFDeaths <- read.table(file="BirdFluDeaths.txt", header = TRUE, sep="\t")
Deaths <- rowSums(BFDeaths[,2:16]) #rowSum对每一行进行求和计算；对BFDeaths的2到16列的每一行求和！
Deaths #还没有名字
names(Deaths) <- BFDeaths[,1]  #names函数将BFCases的第一列内容赋给了Cases，作为Cases的名字；
Deaths #有名字了

Counts <- cbind(Cases, Deaths)  #以列的形式结合向量；
Counts

par(mfrow = c(2,2), mar = c(3, 3, 2, 1))  #par函数进行画布设置
barplot(Cases , main = "Bird flu cases")
barplot(Counts)
barplot(t(Counts), col = gray(c(0.5,1))) #gray(c(0.5,1)))是指两个灰度
barplot(t(Counts), beside = TRUE)  #beside是指两个是上下形式还是左右形式；上面那个默认为F
# t函数将数据进行了转置，如果不转置则是分成了Cases和Deaths两组，经过转置之后以年份分为了6组

在条形图上显示均值和标准差以及外部线框

setwd("C:/RBook/")
Benthic <- read.table(file = "RIKZ2.txt", header = TRUE, sep="\t")
dim(Benthic)
head(Benthic)

Bent.M <- tapply(Benthic$Richness, INDEX=Benthic$Beach, FUN=mean)  #Benthic$Richness表示Benthic文件中Richness这一元素！
head(Bent.M)
Bent.sd <- tapply(Benthic$Richness, INDEX=Benthic$Beach, FUN=sd)
head(Bent.sd)
MSD<- cbind(Bent.M, Bent.sd)
MSD

bp <- barplot(Bent.M, xlab = "Beach", ylab = "Richness", col = rainbow(9), ylim = c(0,20))
bp  # bp表示沿x轴每个条形图的中点。
arrows(bp, Bent.M, bp, Bent.M + Bent.sd, lwd = 1.5, angle=90,length=0.1) 
# arrows给图形加标准差垂直线。 bp, Bent.M表示标准差垂直线的起点；bp, Bent.M + Bent.sd表示标准差垂直线的终点，lwd是线宽，angle=90,length=0.1表示标准差垂直线的顶端线与标准差垂直线的角度和顶端线的长度。
box()  #给图形加一个线框框住

Benth.le <- tapply(Benthic$Richness, INDEX=Benthic$Beach, FUN=length)
Bent.se <- Bent.sd / sqrt(Benth.le)  #计算标准误

stripchart(Benthic$Richness ~ Benthic$Beach, vert = TRUE, pch=1, method = "jitter", jit = 0.05, xlab = "Beach", ylab = "Richness")
points(1:9, Bent.M, pch = 16, cex = 1.5)  #在均值上加点

arrows(1:9, Bent.M, 1:9, Bent.M + Bent.se, lwd = 1.5, angle=90, length=0.1)  #加号表示正的se
arrows(1:9, Bent.M, 1:9, Bent.M - Bent.se, lwd = 1.5, angle=90, length=0.1)  #减号表示负的se

Add lines for sd.

boxplots

setwd("C:/RBook/")
Owls <- read.table(file = "Owls.txt", header= TRUE, sep="\t")
boxplot(Owls$NegPerChick, main = "Negotiation per chick")
par(mfrow = c(2,2), mar = c(3, 3, 2, 1))  #布置画布
head(Owls)
Owls$LogNeg <- log10(Owls$NegPerChick + 1)  #新加了一列名为LogNeg ，把Owls文件中NegPerChick这一子集提取出来，做相应的运算，赋值给左边.
head(Owls)
boxplot(NegPerChick~SexParent, data = Owls)  #把NegPerChick的数据按照SexParent进行画图
boxplot(NegPerChick~FoodTreatment, data = Owls)  #把NegPerChick的数据按照FoodTreatment进行画图
boxplot(NegPerChick~SexParent * FoodTreatment,data = Owls)  #把NegPerChick的数据按照SexParent 和 FoodTreatment进行画图
boxplot(NegPerChick~SexParent * FoodTreatment, names = c("F/Dep","M/Dep","F/Sat","M/Dep"),data = Owls)#把NegPerChick的数据按照SexParent 和 FoodTreatment进行画图，并进行命名；

在条形图上加标签

par(mar = c(2,2,3,3))
boxplot(NegPerChick~Nest, data = Owls, axes = FALSE, ylim = c (-3.5, 9))   #axes表示轴线
axis(2, at = c(0, 2, 4, 6, 8))  #1表示行，2表示列，在图片左边生成一个轴，axis的第一个参数是位置；at参数指的是要标记的刻度；
text(x = 1:27, y = -2, labels = levels(Owls$Nest),cex=0.75, srt=60)     # y=-2表示该标签在y轴的相对位置，放置标签，cex表示字体大小，srt表示倾斜角度！
levels(Owls$Nest)
    

#Boxplot for RIKZ data
setwd("C:/RBook/")
Benthic <- read.table(file = "RIKZ2.txt", header= TRUE, sep="\t")
head(Benthic)
names(Benthic)
Bentic.n <- tapply(Benthic$Richness, Benthic$Beach, FUN =length) #对文件Benthic，以Beach对Richness进行分组，统计每组数据个数
Bentic.n
bp <- boxplot(Richness ~ Beach, data = Benthic, col = "grey",  
              xlab = "Beach", ylab = "Richness")  #boxplot函数可以将信息存储在一个list里面，也就是此处的bp包含好几个变量！
bp
str(bp)

bpmid <- bp$stats[2, ] + (bp$stats[4,] - bp$stats[2,]) / 2
bpmid
text(1:9, bpmid, Bentic.n, col = "white", font = 2)  #x=1:9是向量的长度，y=bpmid表示9元素向量的位置，labels =Bentic.n表示9元素向量的内容
?text

心得！
对于一个函数没有必要去背，比如这里text(1:9, bpmid, Bentic.n, col = "white", font = 2) 是什么意思就很难背下来每一个参数对应哪一个，只要?text一下就可以很容易的对应上

text(x, y = NULL, labels = seq_along(x$x), adj = NULL,
      pos = NULL, offset = 0.5, vfont = NULL,
      cex = 1, col = NULL, font = NULL, ...)  #只要一对应就好了，知道1:9, bpmid, Bentic.n,分别对应着x，y，labels，然后再在帮助文档看一下什么意思就可以了。

点图

很重要的一点，这个点图上下的摆动是随机分布的(这种随机是指限制在组内的随机)，但是左右的摆动是按照他的刻度来的！

?dotchart
setwd("C:/RBook/")
Deer <-read.table("Deer.txt", header = TRUE, sep="\t")
par(mfrow=c(1,2))  #布置画布
dotchart(Deer$LCT, xlab="Length (cm)", ylab = "Observation number")

dotchart(Deer$LCT, groups = factor(Deer$clas1_4))  #groups选项允许数据根据分类变量进行分组！因为Sex变量里面有缺失值，所以这一条命令执行不成功。
Deer$clas1_4
Isna <- is.na(Deer$clas1_4)  #is.na是检测缺失值的函数，若返回值为TRUE则为缺失值，FALSE则为非缺失值；
!Isna  #符号！将Isna进行了一个颠倒。颠倒了之后下一步的作图只会去对TRUE的进行绘图。

dotchart(Deer$LCT[!Isna], groups = factor(Deer$clas1_4[!Isna]), 
         xlab = "Length (cm)", ylab = "Observation number grouped by sex")  


setwd("C:/RBook/")
Benthic <- read.table(file = "RIKZ2.txt", header= TRUE,sep="\t")
head(Benthic)
sum(Benthic$Beach)  #因为Benthic$Beach是数值型向量，所以可以求和；

Benthic$fBeach <- factor(Benthic$Beach)  #factor将数值型向量变为了因子。可以通过sum函数发现Benthic$Beach可以求和，Benthic$fBeach却不可以；
sum(Benthic$fBeach)  #因为Benthic$Beach是因子，所以不能求和；
par(mfrow=c(1,2))
dotchart(Benthic$Richness,groups=Benthic$fBeach,  #如果省略掉Benthic$fBeach <- factor(Benthic$Beach) 这一步，则groups=Benthic$fBeach就应该写成groups=factor(Benthic$Beach)
         xlab="Richness", ylab = "Beach")

Bent.M<-tapply(Benthic$Richness,Benthic$Beach,FUN = mean)
Bent.M

dotchart(Benthic$Richness,groups=Benthic$fBeach,
         gdata = Bent.M, gpch=19, xlab="Richness", ylab="Beach")  # gdata用来表示一个值，这里指的平均值

legend("bottomright",c("values","mean"),pch=c(1,19),bg="red")  #legend函数添加图例，bg表示背景色的颜色

plot函数

setwd("C:/RBook/")
Benthic <- read.table(file = "RIKZ2.txt", header= TRUE,sep="\t")
Benthic$fBeach <- factor(Benthic$Beach)
plot(Benthic$Richness,Benthic$fBeach)
par(mfrow = c(2,2), mar = c(4,4,2,2) +.5)

plot(y = Benthic$Richness, x = Benthic$NAP,
     xlab = "Mean high tide (m)", ylab = "Species richness",
     main = "Benthic data")

M0 <- lm(Richness ~ NAP, data = Benthic)  #lm利用线性回归建立模型，结果存储在列表M0里面
abline(M0)  #abline函数添加拟合线

plot函数中type，axes函数选项

plot(y = Benthic$Richness, x = Benthic$NAP,
     type = "n", axes = FALSE,  #type表示图形的类型,n就表示不画图；axes表示轴线
     xlab = "Mean high tide", ylab = "Species richness")
?plot  #可以看到，如果将type = "n"换成type = "p"就会出图；

     xlab = "Mean high tide", ylab = "Species richness")
points(y = Benthic$Richness, x = Benthic$NAP)

plot(y = Benthic$Richness, x = Benthic$NAP,
     type = "n", axes = FALSE,
     xlab = "Mean high tide", ylab = "Species richness",
     xlim = c(-1.75,2), ylim = c(0,20))
points(y = Benthic$Richness, x = Benthic$NAP)

axis命令绘制轴线

axis(2, at = c(0, 10, 20), tcl = 1)  #axis第一个参数2表示在左边添加轴，at表示添加的刻度；tcl表示z轴上刻度的长短；
axis(1, at = c(-1.75,0,2), labels = c("Sea","Water line","Dunes") )

点/文本/线的添加

setwd("C:/RBook/")
Birds <- read.table(file="loyn.txt", header= TRUE, sep="\t")
Birds$LOGAREA<- log10(Birds$AREA)
Birds$fGRAZE <- factor(Birds$GRAZE)

M0 <- lm(ABUND~ LOGAREA + fGRAZE, data = Birds)
summary(M0)

plot(x = Birds$LOGAREA, y = Birds$ABUND, xlab="Log transformed AREA", ylab="Bird abundance")

改变字体和字体大小

windowsFonts()  #查看目前已经安装好的字体类型
title("Fitted model", cex.main = 2, family = "sans", font.main = 1)  #family指字体类型，font.main指字体大小

#计算五个拟合值
LAR<-seq(-1, 3, by = 0.1)

ABUND1 <- 15.7 +  7.2 * LAR
ABUND2 <- 16.1 +  7.2 * LAR
ABUND3 <- 15.5 +  7.2 * LAR
ABUND4 <- 14.1 +  7.2 * LAR
ABUND5 <- 3.8 +  7.2 * LAR

#把拟合值作为线添加到图形上
lines(LAR, ABUND1, lty = 1, lwd = 1, col =1)
lines(LAR, ABUND2, lty = 2, lwd = 2, col =2)
lines(LAR, ABUND3, lty = 3, lwd = 3, col =3)
lines(LAR, ABUND4, lty = 4, lwd = 4, col =4)
lines(LAR, ABUND5, lty = 5, lwd = 5, col =5)

legend.txt <- c("Graze 1","Graze 2","Graze 3","Graze 4","Graze 5")

图例添加

legend("topleft",  #topleft表示位置
       legend = legend.txt,  #要标注的内容
       col = c(1,2,3,4,5),  #拟合线的颜色
       lty = c(1,1,1,1,5),  #图例线的形式，lty=1表示实心；
       lwd = c(1,2,3,4,5),  #图例线条的粗细
       bty = "o",  #添加围绕图例的盒子
       cex = 0.8  #图例文本大小
)

添加特殊符号

setwd("C:/RBook/")
Whales <- read.table(file="TeethNitrogen.txt", header= TRUE)
head(Whales)

N.Moby1 <- Whales$X15N[Whales$Tooth == "Moby"]  #因为[,]的逗号左边是行，右边是列；但是要求被取值的内容不能是一维度，否则哪里来的行和列？因为此处Whales$X15N是一维数据，所以[]中不能加逗号；
N.Moby1

#N.Moby2 <- Whales[Whales$Tooth == "Moby",]  #因为[,]的逗号左边是行，右边是列；但是要求被取值的内容不能是一维度，否则哪里来的行和列？
#N.Moby2

强烈注意：

关于[]取子集什么时候加逗号，什么时候不加逗号

因为逗号是用来分隔行和列的，逗号前的为行，逗号后的为列；所有加逗号取子集的前提是被取值的内容不是一维数据，有行有列是二维啊；所以，一维数据取子集不用加逗号，因为一维没有行列；

Age.Moby <- Whales$Age[Whales$Tooth == "Moby"]
plot(x = Age.Moby, y = N.Moby, xlab = "Age", ylab = expression(paste(delta^{15}, "N")))

pairs 函数生成多面板的散点图

setwd("C:/RBook/")
Benthic <- read.table(file = "RIKZ2.txt", header= TRUE)
head(Benthic)
class(Benthic)
pairs(Benthic[, 2:9])  #等同于plot(Benthic[,2:9])，因为Benthic是一个数据框，plot函数能够识别它并能调用函数plot.data.frame
#生成的图片对角线显示的是变量的名称，表示对角线上方或下方面板的x轴变量名称，同时还表示对角线左方或右方面板y轴的变量名称
#即每一个图片都可以画一个横轴，也可以画一个纵轴，会跟对角线有交点，交点就是x变量和y变量
#这种多组图有一半的图是多余的，因为每幅图都会出现两次；分别出现在对角线上下方，但坐标轴是颠倒的。

pairs(Benthic[, 2:9], diag.panel = panel.hist, upper.panel = panel.smooth, lower.panel = panel.cor)



# 皮尔逊相关系数应用于pairs函数 -------------------------------------------------
panel.hist <- function(x, ...)
{
  usr <- par("usr"); on.exit(par(usr))
  par(usr = c(usr[1:2], 0, 1.5) )
  h <- hist(x, plot = FALSE)
  breaks <- h$breaks; nB <- length(breaks)
  y <- h$counts; y <- y/max(y)
  rect(breaks[-nB], 0, breaks[-1], y, col="cyan", ...)
}


panel.cor <- function(x, y, digits=2, prefix="", cex.cor)
{
  usr <- par("usr"); on.exit(par(usr))
  par(usr = c(0, 1, 0, 1))
  r <- abs(cor(x, y))
  txt <- format(c(r, 0.123456789), digits=digits)[1]
  txt <- paste(prefix, txt, sep="")
  if(missing(cex.cor)) cex <- 0.8/strwidth(txt)
  text(0.5, 0.5, txt, cex = cex * r)
}

pairs(Benthic[, 2:9], diag.panel = panel.hist,
      upper.panel = panel.smooth, lower.panel = panel.cor)

# pairs函数只能显示双向关系，显示多项关系需要coplot



# coplot 显示多项关系 -----------------------------------------------------------

# * coplot显示三向关系 ----------------------------------------------------------
setwd("C:/RBook/")
Benthic <- read.table(file = "RIKZ2.txt", header= TRUE, sep="\t")

#对于不同Beach绘制Richness对于NAP的散点图；y~x的顺序确定自变量和因变量；as.factor(Beach)把海滩变量作为条件生成图形；
coplot(Richness ~ NAP | as.factor(Beach), pch=19, data = Benthic)  

#对于不同grainsize绘制Richness对于NAP的散点图；
coplot(Richness ~ NAP | grainsize, pch=19, data = Benthic)  

# * coplot 每个面板里加上一条回归线 ---------------------------------------------------
# panel.lm定义在每个面板里如何显示数据
panel.lm = function(x, y, ...) {
  tmp<-lm(y~x,na.action=na.omit)  #线性回归函数lm用来把数据暂时存储在变量tmp中，该分析中任何NAs都将被忽略；
  abline(tmp, lwd = 2)  #abline函数绘制线
  points(x,y, ...)}  #points绘制点

coplot(Richness ~ NAP | as.factor(Beach), pch=19, panel = panel.lm, data=Benthic) #对于不同Beach绘制Richness对于NAP的散点图；
coplot(Richness ~ NAP | as.factor(Beach), pch=19, span =1, panel = panel.smooth, data=Benthic)


# * coplot 显示四向关系 ---------------------------------------------------------
# panel.lm定义在每个面板里如何显示数据
panel.lm = function(x, y, ...) {
  tmp<-lm(y~x,na.action=na.omit)  #线性回归函数lm用来把数据暂时存储在变量tmp中，该分析中任何NAs都将被忽略；
  abline(tmp, lwd = 2)  #abline函数绘制线
  points(x,y, ...)}  #points绘制点

#PH对应于SDI，海拔和森林覆盖之间的关系
setwd("C:/RBook/")
pHEire<-read.table(file="SDI2003.txt", header=TRUE, sep="\t")
pHEire$LOGAlt <- log10(pHEire$Altitude)
pHEire$fForested <- factor(pHEire$Forested)
coplot(pH ~ SDI | LOGAlt * fForested, panel=panel.lm,data=pHEire)
coplot(pH ~ SDI | LOGAlt * fForested, panel=panel.lm,data=pHEire, number=2)  #LOGAlt区间的位置和数目可以通过number来调控；


# * 增加coplot的修饰 -----------------------------------------------------------
panel.lm2=function(x,y,col.line="black",lwd=par("lwd"),
                   lty=par("lwt"), ...){
  tmp<-lm(y~x,na.action=na.omit)
  points(x,y, ...)
  abline(tmp,col=col.line,lwd=lwd,lty=lty)}


CI <- co.intervals(pHEire$LOGAlt,3)
GV <- list(CI,c(2,1))

pHEire$Temperature

#根据不同的温度，将点变成不同颜色
pHEire$Temp2 <- cut(pHEire$Temperature, breaks = 2)  #breaks = 2，所以cut函数将温度数据分成两部分；
pHEire$Temp2
pHEire$Temp2.num <- as.numeric(pHEire$Temp2)  #将pHEire$Temp2从因子转为数值，因为因子不能作为色彩和灰度值；
pHEire$Temp2.num
coplot(pH ~ SDI | LOGAlt * fForested, panel=panel.lm, data=pHEire, given.values = GV, cex=1.5,pch=19, col=gray(pHEire$Temp2.num/3))


# 图形排列 --------------------------------------------------------------------
#mfrow命令可以在一个屏幕上绘制多个图形；而layout可以生成复杂的图形排列；
MyLayOut <- matrix(c(2,0,1,3), nrow = 2, ncol=2, byrow = TRUE)  #matrix函数生成一个矩阵，2行，2列，按行排列；
MyLayOut
nf <- layout(mat = MyLayOut,widths = c(3, 1),  # widths 指定列宽；
             heights = c(1, 3), respect = TRUE)  # height 指定行高，respect = TRUE确保垂直方向上的一个单位与水平方向上的一个单位相同；

layout.show(nf)

xrange<-c(min(Benthic$NAP),max(Benthic$NAP))
yrange<-c(min(Benthic$Richness),max(Benthic$Richness))

#生成第一个图形
par(mar=c(4,4,2,2))
plot(Benthic$NAP,Benthic$Richness,xlim=xrange,ylim=yrange, xlab = "NAP", ylab="Richness")

#生成第二个图形
par(mar=c(0,3,1,1))
boxplot(Benthic$NAP,horizontal=TRUE,axes=F, frame.plot=F,ylim=xrange,space=0)

#生成第三个图形
par(mar=c(3,0,1,1))
boxplot(Benthic$Richness,axes=F,ylim=yrange,space=0,horiz=TRUE)



# —— ----------------------------------------------------------------------

# 格包简介 --------------------------------------------------------------------

setwd("C:/RBook")
Env <- read.table(file="RIKZENV.txt",header = TRUE)
names(Env)
str(Env)
sum(is.na(Env))  #2975 NAs

## where are the NAs?

# Stations in each area
table(Env$Station,Env$Area)

tapply(Env$SAL,list(Env$Station,Env$Year,Env$Month),length)
tapply(Env$SAL,list(Env$Station,Env$Month),length)
tapply(Env$SAL,list(Env$Area,Env$Month),length)

# count values
ndata<-function(x) sum(!is.na(x))
tapply(Env$SAL,list(Env$Station,Env$Month),ndata)
tapply(Env$SAL,list(Env$Area,Env$Month),ndata)

# count NAs
nnas<-function(x) sum(is.na(x))
tapply(Env$SAL,list(Env$Station,Env$Month),nnas)
tapply(Env$SAL,list(Env$Area,Env$Month),nnas)



# 多面板散点图：xyplot -----------------------------------------------------------
library(lattice)
Env$MyTime<-Env$Year+Env$dDay3/365

#根据station不同来绘制SAL对MyTime的图形；
xyplot(SAL ~ MyTime | factor(Station), type="l",  # y轴 ~ x轴，管道命令后为条件变量(此处为station)，条件变量通常为因子；
       strip = function(bg, ...)  #...表示将信息传递给其他特定的函数
         strip.default(bg = 'white', ...),  #对每一个带使用白色背景
       col.line=1,  #图形线条颜色
       data = Env)

xyplot(SAL ~ MyTime | factor(Station), data = Env)
xyplot(SAL ~ MyTime | factor(Station), type="l", strip = TRUE, col.line=1, data = Env)
xyplot(SAL ~ MyTime | factor(Station), type="l", strip = FALSE, col.line=1, data = Env) #查看strip功能，strip用来绘制带状；


# 多面板盒型图：bwplot -----------------------------------------------------------
bwplot(SAL ~ factor(Month) | Area,
       strip = strip.custom(bg = 'white'),
       cex = .5, layout = c(2, 5),
       data = Env, xlab = "Month", ylab = "Salinity",
       par.settings = list(
         box.rectangle = list(col = 1),
         box.umbrella  = list(col = 1),
         plot.symbol   = list(cex = .5, col = 1)))

bwplot(SAL ~ factor(Month) | Area, layout = c(2, 5), data = Env) #若采用黑白作图，则上面的代码可以简化为这个；



# 多面板点图：dotplot -----------------------------------------------------------
dotplot(factor(Month) ~ SAL  | Station,
        subset = Area=="OS", jitter.x = TRUE,  #subset是用来选定全部数据的一个子集；jitter.x = TRUE表示当多个观察值在同一月份相同的时候，在水平方向增加一些变化；
        data = Env, strip = strip.custom(bg = 'white'),
        col = 1, cex = 0.5, ylab = "Month",
        xlab = "Salinity")

dotplot(SAL ~ factor(Month) | Station,    
        subset = Area=="OS", jitter.x = T,
        data = Env, strip = strip.custom(bg = 'white'),
        col = 1, cex = 0.5, xlab = "Month",
        ylab = "Salinity")    #与上一个相比将x轴y轴对调一下；



# 多面板直方图: histogram -------------------------------------------------------
histogram( ~ SAL | Station, data = Env,
           subset = Area=="OS", layout = c(1,4),  #layout里面为1列2行的布局，与先行后列不同
           nint = 30, xlab = "Salinity", strip = FALSE,  #nint是条形的数目
           strip.left = TRUE, ylab="Frequencies")  #通过strip = FALSE, strip.left = TRUE将带状移到了面板左边；


#Example 1
xyplot(SAL ~ Month | Year, data = Env,  
       type = c("p"), subset = (Station =="GROO"),  #type有三个函数，分别是p，g，smooth；相当于函数执行了panel.xyplot，panel.grid，panel.loess
       xlim = c(0, 12), ylim = c(0, 30),pch = 19,
       panel = function (...){
         panel.xyplot(...)  #画出数据点
         panel.grid(..., h = -1, v = -1)  #panel.grid加网格线，h和v为正值表示水平与垂直网格线的数目；若为负值则网格与轴坐标对其；
         panel.loess(...)}, span = 0.9)  #panel.loess加拟合线，span表示拟合线的平滑程度；

#Example 2
#使用不同的颜色和大一点的点来表示潜在异常值，方法1和方法2是殊途同归
setwd("C:/RBook")
library(lattice)
Env <- read.table(file="RIKZENV.txt",header = TRUE)
#方法1
dotplot(factor(Month) ~ SAL | Station, pch = 16,
        subset = (Area=="OS"), data = Env,
        ylab = "Month", xlab = "Salinity",
        panel = function(x,y,...) {
          q1 <- summary(x,na.rm=TRUE)[2]
          q3 <- summary(x,na.rm=TRUE)[5]
          R <- q3 - q1
          L <- median(x,na.rm=TRUE) - 3 * (q3 - q1)
          MyCex <- rep(0.4,length(y))
          MyCol <- rep(1,length(y))
          MyCex[x < L] <- 1.5
          MyCol[x < L] <- 2
          panel.dotplot(x,y, cex = MyCex, col = MyCol, ...)})

# 方法2
dotplot(factor(Month) ~ SAL | Station, pch = 16,
        subset = (Area=="OS"), data = Env,
        ylab = "Month", xlab = "Salinity",
        panel = function(x, y, ...) {
          Q <- quantile(x, c(0.25, 0.5, 0.75) ,  #
                        na.rm = TRUE)
          R <- Q[3] - Q[1]
          L <- Q[2] - 3 * (Q[3] - Q[1])  #L为盐度数据的截止水平，即中位数减去第三四分位数和第一四分位数差的三倍
          MyCex <- rep(0.4, length(y))
          MyCol <- rep(1, length(y))
          MyCex[x < L] <- 1.5  #当盐度低于标准修改点的尺寸
          MyCol[x < L] <- 2    #当盐度低于标准修改点的颜色
          panel.dotplot(x, y, cex = MyCex, 
                        col = MyCol, ...)})



# 多面板 PCA 分析双标图 -----------------------------------------------------------
setwd("C:/RBook")
Sparrows<-read.table(file="Sparrows.txt", header=TRUE)
names(Sparrows)
str(Sparrows)
library(lattice)
xyplot(Wingcrd ~ Tarsus | Species * Sex,  #管道命令后面是两个条件变量
       xlab = "Axis 1", ylab = "Axis 2", data = Sparrows,
       xlim = c(-1.1, 1.1), ylim = c(-1.1, 1.1),
       panel = function(subscripts, ...){
         zi <- Sparrows[subscripts, 3:8]
         di <- princomp(zi, cor = TRUE)
         Load <- di$loadings[, 1:2]
         Scor <- di$scores[, 1:2]
         panel.abline(a = 0, b = 0, lty = 2, col = 1)
         panel.abline(h = 0, v = 0, lty = 2, col = 1)
         for (i in 1:6){
           llines(c(0, Load[i, 1]), c(0, Load[i, 2]), 
                  col = 1, lwd = 2)
           ltext(Load[i, 1], Load[i, 2], 
                 rownames(Load)[i], cex = 0.7)}
         sc.max <- max(abs(Scor))
         Scor <- Scor / sc.max
         panel.points(Scor[, 1], Scor[, 2], pch = 1,
                      cex = 0.5, col = 1)
       })



# 三维点图/表面图/等高线图 -----------------------------------------------------------
setwd("C:/RBook")
library(lattice)
Env <- read.table(file="RIKZENV.txt",header = TRUE)   
cloud(CHLFa ~ T * SAL | Station, data = Env,
      screen = list(z = 105, x = -70),
      ylab = "Sal", xlab = "T", zlab = "Chl. a",
      ylim = c(26,33),
      subset = (Area=="OS"),
      scales = list(arrows = FALSE))




# 常见问题 --------------------------------------------------------------------


# * 改变面板顺序 ----------------------------------------------------------------
setwd("C:/RBook")
Hawaii <- read.table("waterbirdislandseries.txt", header = TRUE)
library(lattice)
names(Hawaii)
Birds <- as.vector(as.matrix(Hawaii[,2:9]))  #as.vector命令不能应用于数据框；所以as.matrix命令首先将数据转化为一个矩阵，然后再由as.vector将其转化为一个长向量！
Time <- rep(Hawaii$Year, 8)  #rep函数的作用是生成一个单独的长向量；
MyNames <- c("Stilt_Oahu","Stilt_Maui","Stilt_Kauai_Niihau","Coot_Oahu","Coot_Maui","Coot_Kauai_Niihau","Moorhen_Oahu","Moorhen_Kauai")
ID <- rep(MyNames, each = 48)

xyplot(Birds ~ Time|ID, ylab = "Bird abundance", layout = c(3, 3), type = "l", col = 1)

#通过改变levels顺序就可以改变这些面板的顺序；
ID2 <- factor(ID, levels=c("Stilt_Oahu","Stilt_Kauai_Niihau","Stilt_Maui","Coot_Oahu","Coot_Kauai_Niihau","Coot_Maui","Moorhen_Oahu","Moorhen_Kauai"))
xyplot(Birds ~ Time|ID2, ylab = "Bird abundance",layout = c(3, 3), type = "l", col = 1)


# * 改变坐标轴界限刻度 -------------------------------------------------------------
xyplot(Birds ~ Time|ID2, ylab = "Bird abundance",
       layout = c(3, 3), type = "l", col = 1,
       scales = list(x = list(relation = "same"),
                     y = list(relation = "free"),
                     tck=-1))  #tck是刻度，负朝里，正朝外；数字代表长短

# * 一个面板中绘制多条线 ------------------------------------------------------------
#将同一种鸟的所有时间序列陈列在一个单独的面板中
Species <-rep(c("Stilt","Stilt","Stilt","Coot","Coot","Coot","Moorhen","Moorhen"),each = 48)
xyplot(Birds ~ Time|Species, ylab = "Bird abundance",
       layout = c(2, 2), type = "l", col = 1,
       scales = list(x = list(relation = "same"),
                     y = list(relation = "free")),
       groups = ID, lwd=c(1,2,3))




# 更新图形 --------------------------------------------------------------------
MyPlot <- xyplot(Birds ~ Time|Species, ylab = "Bird abundance",
                 layout = c(2, 2), type = "l", col = 1,
                 scales = list(x = list(relation = "same"),
                               y = list(relation = "free")),
                 groups = ID, lwd=c(1,2,3))

print(MyPlot)
update(MyPlot, layout=c(3,3)) #对象的很多属性可以通过update函数来修改！可以保持原始图不变。