《R語言初學者指南》學習筆記
2020-06-18 本文已影响0人
余顺太
第一章 引言
#字体用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函数来修改！可以保持原始图不变。
