上節(jié)課數(shù)量生態(tài)學(xué)筆記||緒論我們簡單了解《數(shù)量生態(tài)學(xué)》的基本內(nèi)容渡紫，特別介紹了書中用到的數(shù)據(jù)集Doubs、甲螨數(shù)據(jù)集考赛。關(guān)于R并未做過多的介紹惕澎，因為這是一本生態(tài)學(xué)的書。但是關(guān)于學(xué)習(xí)方法我推薦一種卡片學(xué)習(xí)法：

將知識點整理在一個小的可用隨身攜帶的卡片上欲虚，可以隨時翻閱集灌，可以建立鏈接。也就是學(xué)習(xí)就像拼圖，在整理欣喧、記憶腌零、鏈接中形成自己的知識樹。

上節(jié)課的卡片其實就是開始的那張導(dǎo)圖唆阿，我鼓勵制作自己的學(xué)習(xí)卡片（手寫）益涧。可以包括生態(tài)學(xué)知識點驯鳖，R函數(shù)闲询，并不是為了記憶而是為了聯(lián)系。沒用存量談不上體系浅辙。

今天我們來學(xué)習(xí)書中的第二章：探索性數(shù)據(jù)分析扭弧。

在我們提到數(shù)據(jù)分析的時候還腦海里閃現(xiàn)的往往是簡潔的報表、漂亮的數(shù)據(jù)圖记舆，再不濟也會聯(lián)想到假設(shè)檢驗（Hypothesis Testing）和建模（modeliling）鸽捻。然而，數(shù)據(jù)采集完之后泽腮，數(shù)據(jù)表整理好之后御蒲，我要做的第一步并不是去做復(fù)雜的統(tǒng)計分析。而是要做一些統(tǒng)計來了解數(shù)據(jù)的概況诊赊，對數(shù)據(jù)有一個大致的認識厚满。用這認識去指導(dǎo)我們后面的分析實踐，這就是數(shù)據(jù)探索碧磅。如SPSS中就有一個菜單叫做描述統(tǒng)計碘箍。本章基本上也是屬于這個范疇。

數(shù)據(jù)概況

首先我們載入數(shù)據(jù)：

rm(list=ls())
setwd("D:\\Users\\Administrator\\Desktop\\RStudio\\數(shù)量生態(tài)學(xué)\\Run")

#導(dǎo)入物種多度數(shù)據(jù)
spe<-read.csv("../DATA/DoubsSpe.csv",row.names = 1)
str(spe)
#導(dǎo)入環(huán)境數(shù)據(jù)
env<-read.csv("../DATA/DoubsEnv.csv",row.names = 1)
#導(dǎo)入空間坐標(biāo)數(shù)據(jù)
spa<-read.csv("../DATA/DoubsSpa.csv",row.names = 1)

我們對物種群落數(shù)據(jù)做一個簡單的描述統(tǒng)計续崖，同時也是看看我們的和數(shù)據(jù)格式是否正常敲街。

# 基礎(chǔ)函數(shù)
# ********
spe #在控制臺顯示整個數(shù)據(jù)框的內(nèi)容，但對于大樣本的數(shù)據(jù)框
#并不建議直接顯示
spe[1:5,1:10]      #只展示前5行和前10列
head(spe)           #只展示前幾行
nrow(spe)           #提取數(shù)據(jù)框總行數(shù)
ncol(spe)           #提取數(shù)據(jù)框總列數(shù)
dim(spe)            #提取數(shù)據(jù)框的維度（顯示數(shù)據(jù)框多少行严望，多少列）
colnames(spe)       #提取列名多艇，在這里是物種名
rownames(spe)       #提取行名，在這里一行代表一個樣方
summary(spe)         #以列為單位像吻，對列變量進行描述性統(tǒng)計
#比較多度的中位值和平均值峻黍。大部分是對稱分布嗎？

如果多度分布是對稱的拨匆，中位數(shù)應(yīng)該和平均值差別不大姆涩。大家看這里的數(shù)據(jù)，顯然多數(shù)數(shù)據(jù)并不是對稱的惭每。

# 多度數(shù)據(jù)總體分布情況
# *******************
# 整個多度數(shù)據(jù)值的范圍
range(spe)
# 計算每種多度值的數(shù)量
ab <- table(unlist(spe))
ab
# 所有種混和在一起的多度分布柱狀圖
barplot(ab, las=1, xlab="多度等級", ylab="頻度", col=gray(5:0/5))
# 多度數(shù)據(jù)中0值的數(shù)量
sum(spe==0)
# 多度數(shù)據(jù)中0值所占比例
sum(spe==0) / (nrow(spe)*ncol(spe))
#請觀察多度頻率分布柱狀圖骨饿，如何解讀為什么0值（缺失）在數(shù)據(jù)框內(nèi)頻
#率這么高?

其實造成缺失的因素有很多亏栈，但是有兩種需要我們的注意：

• 真實的環(huán)境適合這個物種生存，只是我們采樣的時候沒采到（比如人家冬眠了宏赘，出去玩了绒北，尚未遷徙到這里）。
• 真實的環(huán)境不適合這個物種生存察署，在這里生存就會被淘汰闷游。

所以對于零值我們要小心處理，關(guān)鍵還是理解數(shù)據(jù)的生物學(xué)意義贴汪。

樣方的分布

數(shù)據(jù)探索也是也個數(shù)據(jù)和模型相互磨合的過程脐往，不僅看用來描述我們實驗本身的數(shù)據(jù)，也可以用來描述實驗設(shè)計扳埂。

 樣方位置地圖
# **************
# 生成空的繪圖窗口（橫縱坐標(biāo)軸比例1:1业簿，帶標(biāo)題）
# 從spa數(shù)據(jù)框獲取地理坐標(biāo)x和y
plot(spa, asp=1, type="n", main="樣方位置",
     xlab="x坐標(biāo) (km)", ylab="y坐標(biāo) (km)")
# 加一條連接各個樣方點的藍色線（代表Doubs河）
lines(spa, col="light blue")
# 添加每個樣方的編號
text(spa, row.names(spa), cex=0.8, col="red")
# 添加文本
text(70, 10, "上游", cex=1.2, col="red")
text(20, 120, "下游", cex=1.2, col="red")

30g個樣方沿著Doubs河的空間分布。繪制這幅圖用到的plot函數(shù)是R的基礎(chǔ)繪圖函數(shù)聂喇∠皆矗可以?plot()查看其幫助文檔，asp是用來調(diào)整繪圖版面的長寬比列的希太。

下面我們把物種數(shù)據(jù)映射到采樣點之上，看看物種是怎樣隨著河流變化的酝蜒。

#某些魚類的分布地圖
# ******************
# 將繪圖窗口分割為4個繪圖區(qū)域誊辉，每行兩個
par(mfrow=c(2,2))
plot(spa, asp=1, col="brown", cex=spe$TRU, main="褐鱒",
     xlab="x坐標(biāo) (km)", ylab="y坐標(biāo) (km)")
lines(spa, col="light blue")
plot(spa, asp=1, col="brown", cex=spe$OMB, main="茴魚",
     xlab="x坐標(biāo) (km)", ylab="y坐標(biāo) (km)")
lines(spa, col="light blue")
plot(spa, asp=1, col="brown", cex=spe$BAR, main="鲃魚",
     xlab="x坐標(biāo) (km)", ylab="y坐標(biāo) (km)")
lines(spa, col="light blue")
plot(spa, asp=1, col="brown", cex=spe$BCO, main="歐鳊",
     xlab="x坐標(biāo) (km)", ylab="y坐標(biāo) (km)")
lines(spa, col="light blue")
#觀察所生成的4張圖，你就會明白為什么Verneaux 選擇這4種魚類作為不同區(qū)域的生態(tài)指示種亡脑，看了后面將要展示的環(huán)境因子空間分布情況會更清楚堕澄。

從這個圖上我們清楚地看到，褐鱒霉咨、茴魚蛙紫、鲃魚、歐鳊的多度是沿著Doubs河從上游到下游分布的途戒，也就明白為什么Verneaux 選擇這4種魚類作為不同區(qū)域的生態(tài)指示種坑傅。注意之前提到的零值問題，這里是同一條河流不會有遷移的障礙喷斋，這幾種魚的生活史也較接近不存在有冬眠的不一致的情況唁毒。

另一個引起我們注意的就是plot()函數(shù)的參數(shù)cex的用法，它的作用是定義數(shù)據(jù)點標(biāo)識的大小星爪。提問浆西，為什么這個標(biāo)識是圓形的而不是其他的呢？可以調(diào)嗎顽腾？

每個物種在多少樣方中出現(xiàn)近零？，我們可以看看物種的相對頻度。

# 比較物種頻度
# **************
# 計算每個物種出現(xiàn)的樣方數(shù)
#按列進行計數(shù)久信，因此函數(shù)apply()第二個參數(shù)MARGIN應(yīng)該設(shè)定為2
spe.pres <- apply(spe > 0, 2, sum)
# 按照升序的方式重新排列結(jié)果
sort(spe.pres)
# 計算頻度百分比
spe.relf <- 100*spe.pres/nrow(spe)
round(sort(spe.relf), 1)    # 設(shè)置排列結(jié)果為1位小數(shù)
#繪柱狀圖
par(mfrow=c(1,2))       # 將繪圖窗口垂直一分為二
hist(spe.pres, main="物種出現(xiàn)數(shù)", right=FALSE, las=1,
     xlab="出現(xiàn)數(shù)", ylab="物種數(shù)量",
     breaks=seq(0,30,by=5), col="bisque")
hist(spe.relf, main="物種相對頻度", right=FALSE, las=1,
     xlab="出現(xiàn)率(%)", ylab="物種數(shù)量",
     breaks=seq(0, 100, by=10), col="bisque")

我問鼓寺，這兩個圖的縱軸都是“物種數(shù)量”為什么最大值還不一樣呢噪馏？頻度圖讓我們了解每個物種存在于多少個樣方內(nèi)。接下來我們看看每個樣方內(nèi)存在多少物種（物種的豐度）。思考頻度與豐度的差異悬槽。主義apply函數(shù)的應(yīng)用，apply函數(shù)家族在R中應(yīng)用很普遍刘陶。

# 樣方比較：物種豐富度
# ********************
# 計算每個樣方內(nèi)物種數(shù)
# 以行匯總动壤，apply()函數(shù)第二個參數(shù)MARGIN應(yīng)該設(shè)定為1
sit.pres <- apply(spe > 0, 1, sum)
#按照升序的方式重新排列結(jié)果
sort(sit.pres)
par(mfrow=c(1,2))   #將繪圖窗口垂直一分為二
# 繪制樣方沿著河流的分布位置和所含物種豐富度
plot(sit.pres,type="s", las=1, col="gray",
     main="物種豐富度-上下游的梯度",
     xlab="樣方沿著河流的位置", ylab="物種豐富度")
text(sit.pres, row.names(spe), cex=.8, col="red")
# 使用地理坐標(biāo)繪制氣泡地圖
plot(spa, asp=1, main="物種豐富度地圖", pch=21, col="white",
     bg="brown", cex=5*sit.pres/max(sit.pres), xlab="x坐標(biāo) (km)",
     ylab="y坐標(biāo) (km)")
lines(spa, col="light blue")
#你能否辨析沿著河流哪里是物種豐富度的熱點地區(qū)？

我們可以清楚地看出沿河流物種的整體分布酥诽。

最后我們用vegan包中的diversity()函數(shù)計算生物多樣性指數(shù)鞍泉。

#計算生物多樣性指數(shù)
# *****************
# 載入所需要的vegan程序包
library(vegan) # 如果未載入，需要執(zhí)行這一步
#訪問diversity() 幫助界面
?diversity
N0 <- rowSums(spe > 0)               #物種豐富度
H <- diversity(spe)                    # Shannon熵指數(shù)
N1 <- exp(H)                        # Shannon 多樣性指數(shù)
N2 <- diversity(spe, "inv")              # Simpson多樣性指數(shù)
J <- H/log(N0)                          # Pielou 均勻度
E1 <- N1/N0                            # Shannon均勻度 (Hill比率)
E2 <- N2/N0                            # Simpson均勻度 (Hill比率)
div <- data.frame(N0, H, N1, N2, E1, E2, J)
div

大家看到在計算的結(jié)果中肮帐，第8個樣方幾個指數(shù)出現(xiàn)了inf看看原始數(shù)據(jù)spe的第8個樣方有什么規(guī)律咖驮，說出為什么會出現(xiàn)INF。

環(huán)境數(shù)據(jù)

現(xiàn)在我們已經(jīng)對物種數(shù)據(jù)有了基本的了解训枢，下面我們看一下環(huán)境數(shù)據(jù)托修。

# 部分環(huán)境變量的氣泡地圖
# *******************************************
par(mfrow=c(2,2))
plot(spa, asp=1, main="海拔", pch=21, col="white", bg="red",
     cex=5*env$alt/max(env$alt), xlab="x", ylab="y")
lines(spa, col="light blue")
plot(spa, asp=1, main="流量", pch=21, col="white", bg="blue",
     cex=5*env$deb/max(env$deb), xlab="x", ylab="y")
lines(spa, col="light blue")
plot(spa, asp=1, main="氧含量", pch=21, col="white", bg="green3",
     cex=5*env$oxy/max(env$oxy), xlab="x", ylab="y")
lines(spa, col="light blue")
plot(spa, asp=1, main="硝酸鹽濃度", pch=21, col="white", bg="brown",
     cex=5*env$nit/max(env$nit), xlab="x", ylab="y")
lines(spa, col="light blue")
#哪幅圖最能展示上下游的梯度？如何解釋其他環(huán)境變量的空間分布格局恒界？

我們可以說是海拔反映了環(huán)境變量的梯度睦刃。流量和海拔都好解釋，想一下含氧量和硝酸鹽濃度只靠這幾張圖能解釋嗎十酣？在（150,200）處的硝酸鹽濃度很高而氧含量很低涩拙，為什么？

環(huán)境變量量沿河流分布情況耸采。

#線條圖
# *****
par(mfrow=c(2,2))
plot(env$das, env$alt, type="l", xlab="離源頭距離 (km)",
     ylab="海拔 (m)", col="red", main="海拔")
plot(env$das, env$deb, type="l", xlab="離源頭距離 (km)",
     ylab="流量 (m3/s)", col="blue", main="流量")
plot(env$das, env$oxy, type="l", xlab="離源頭距離 (km)",
     ylab="氧含量 (mg/L)", col="green3", main="氧含量")
plot(env$das, env$nit, type="l", xlab="離源頭距離 (km)",
     ylab="硝酸鹽濃度 (mg/L)", col="brown", main="硝酸鹽濃度")

如果要了解任意任意兩個環(huán)境變量之間的關(guān)系兴泥，我們可以使用強大的矩陣散點圖繪制函數(shù)pairs().

# 所有變量對之間的二維散點圖
# **************************
#載入自編的函數(shù)R腳本
source("panelutils.R")  # panelutils.R腳本文件必須與當(dāng)前R工作空間在同一文件
#夾下
# 帶頻度分布的柱狀圖和光滑擬合曲線的雙變量散點圖
op <- par(mfrow=c(1,1), pty="s")
pairs(env, panel=panel.smooth, diag.panel=panel.hist,
      main="雙變量散點圖（帶頻度分布圖和平滑曲線）")
par(op)
#從柱狀圖能否看出哪些變量符合正態(tài)分布？
#需要注意的是虾宇，對于回歸分析和典范排序搓彻，并沒有要求解釋變量符合正態(tài)
#分布。是否有很多散點圖顯示出變量之間的線性關(guān)系或至少是單調(diào)關(guān)系文留？

簡單的轉(zhuǎn)化可以改善某些變量的數(shù)據(jù)分布好唯，另外變量之間的剛量不同很多分析要將其標(biāo)準(zhǔn)化。

# 某個環(huán)境變量簡單轉(zhuǎn)化
# ********************
range(env$pen)
# 坡度變量對數(shù)轉(zhuǎn)化(y = ln(x))
# 比較轉(zhuǎn)化前后數(shù)值的柱狀圖和箱線圖
par(mfrow=c(2,2))
hist(env$pen, col="bisque", right=FALSE, main="坡度頻度分布圖",ylab="頻度",xlab="坡度")
hist(log(env$pen), col="light green", right=F, main="對數(shù)化坡度頻度分布圖",ylab="頻度",xlab="對數(shù)化坡度")
boxplot(env$pen, col="bisque", main="坡度箱線圖", ylab="坡度")
boxplot(log(env$pen), col="light green", main="對數(shù)化坡度箱線圖",
        ylab="對數(shù)化坡度")

# 所有環(huán)境變量的標(biāo)準(zhǔn)化
# *********************
# 中心化和標(biāo)準(zhǔn)化=標(biāo)準(zhǔn)化變量 (z-scores)
env.z <- decostand(env, "standardize")
apply(env.z, 2, mean)   # 平均值 = 0
apply(env.z, 2, sd)     # 標(biāo)準(zhǔn)差 = 1
# 使用scale()函數(shù)也可以運行相同的標(biāo)準(zhǔn)化(輸出的是矩陣)
env.z <- as.data.frame(scale(env))

探索性數(shù)據(jù)分析|百科
卡片學(xué)習(xí)法
學(xué)習(xí)卡片大法
Detailed Table