在 R 當中有非常多的套件，可以提供非常美觀或者複雜的視覺化效果，但如果今天只想要簡單呈現數據間的關係，就可以使用內建的繪圖系統，而不需要額外呼叫其他的套件。 在此我們使用 base R 來進行散點圖的繪製與迴歸分析。 首先準備要分析的數據資料。 R 內建有許多的資料集，可以使用 data() 來查看 R 的內建資料集與各資料集的簡短說明。 如果要查看所有套件的內建資料集，可以使用 data(package = .packages(all.available = TRUE)) 指令。 這邊使用的是 "trees" 資料集。 data(trees) head(trees)

image:https://i.imgur.com/mVZVshx.png :::warning 更新：現在氣象局的 CODiS 氣候資料服務系統已經更改網頁介面，這篇介紹的方法已經沒辦法再使用，只能用 python 自動抓取，或者手動下載資料再用 R 整理，目前尚無解。 ::: 風花圖，又叫做風玫瑰圖（wind rose），是一種呈現風向、風速與頻度分布的資訊繪圖，從圖面上可以了解某地在特定時間的風向、風速概要，在設計機場跑道方向或者設計建築物座向時都會用到。 中央氣象局氣候年報所使用的風花圖圖例。

手上的數據如果有許多變數，又想快速了解數據之間、變數之間的關係，那麼使用主成分分析可以達到這個效果。 主成分分析（principle component analysis, PCA），使用線性變換中的正交變換（orthogonal transformation）方法，對變數測值做線性變換，投影出線性不相關的變數，稱為主成分（principle components）。我們希望能透過主成分分析，將多個變數用比較少的主成分來解釋（降維），最大化變異數的同時又盡可能最小化流失的訊息。 首先建立數據的變異數—共變異數矩陣（variance-covariance matrix），之後進行特徵分解（eigendecomposition）或奇異值分解（singular value decomposition, SVD），最後得到特徵向量（主成分向量的方向）與特徵值（主成分向量的大小），之後就可以選擇數個（通常是 2 至 3 個）主成分將數據視覺化。 現在我們可以使用 R 的內建或者特定套件功能來進行主成分分析以及後續的視覺化。 用內建函數做 PCA R 有兩個內件函數可以用來做主成分分析，分別是 princomp() 與 prcomp()，兩者有些不同，而目前官方是建議以 prcomp() 函數來做主成分分析。

Microsoft Office 裡的萬用字元有問號 ? 與星號 *。 要取代這些萬用字元，在 Linux 系統或者 Notepad++ 中可以在前面加上反斜線 \；而==在 MS Office 中則是要加上 tilde，即波浪符號 ~==。 例如要取代 ? 時，搜尋框要輸入 ~?。 <span style="font-size:30px">🐕‍🦺</span><font color="dcdcdc">2023.07.16</font>

因為 臺灣物種名錄 網站搬家 + 大改版，原本的物種名錄資料格式也有很大的變動，例如名錄檔案分為「學名」與「物種」兩個檔案，前者記錄了更多的資訊如種下階層（subsp., var., fo.）、原始發表文獻、地位等等。對於整理調查資料來說，應該使用物種名錄檔案即可；但因為學名名錄檔才有紀錄是否為接受名（usage_status），所以我們也只好使用這個檔案來製作本機端的資料庫。如果新版網站延續之前 taibnet 的慣例，那名錄資料會每月更新一次，可以定期下載來更新本地端的比對用物種資料庫。 首先下載檔案： wget https://taicol.tw/static/upload/TaiCOL_name_20230818.zip unzip TaiCOL_name_20230818.zip 查看每一行的名稱，共有 49 行，可以使用 awk 指令查看： awk -F, 'NR==1 {for (i=1; i<=NF; i++) {printf("$%d %s\n", i, $i)}}' TaiCOL_name_20230818.csv $1 name_id

在用 R 來畫地形圖——以新竹縣市為例當中，使用 tmap 套件來視覺化新竹縣市的分層設色地形圖。如果想要在旁邊加入一張小地圖，呈現大新竹地區在台灣的位置應跟怎麼做呢？ 之前有提過，tmap 套件源於 grid 套件。因此可以用 grid 中的 viewport() 功能（類似 ggplot2 中的 annotation_custom() 功能），把已經畫好的小地圖，依照指定的視區（就是插圖的大小與位置）畫到主圖上。 先畫一個主圖： colss <- RColorBrewer::brewer.pal(n = 10, name = "RdYlGn") colss1 <- c("#00341b", rev(colss), "blueviolet") HC_map_base <- tm_shape(HC) + tm_borders(lwd = 2) + tm_graticules(n.x = 4, alpha = 0.1, labels.size = 1) +

R 可以完成許多資料的視覺化，地理資料的分析與視覺化也可以由 R 來做。如果熟悉 R 的操作，就可以直覺地調整繪圖參數與算繪（rendering），不用在視覺化軟體上拉來拉去，平白消耗電腦效能。 在 R 當中有許多套件可以做簡單的地理資訊處理（geoprocessing），例 raster、sf（sf 的使用簡介）；視覺化套件主要有 tmap（tmap 的使用簡介）以及搭配 ggplot2 使用的 ggspatial（ggspatial 的使用簡介）。 在這邊我們以 raster 與 sf 分別讀入網格與向量資料，並以 tmap 為視覺化的主要套件來製圖。 這篇用的程式碼參考與修改自此。 :::info :::spoiler {state="open"} :::

:::info :::spoiler {state="open"} 題次 ::: 承接上一篇〈生物多樣性的估計量整理（上）：緒論與 α 多樣性〉，這邊整理的是常用的 β 多樣性估計量。 β 多樣性主要探討在同時間的不同群集之間，或者某群集在不同時間的差異有多大，因此許多量化群集間的差異程度（differentiation）、相似程度（similarity）、重複程度（overlap）或轉換程度（turnover）方法被用以估計 β 多樣性；除此之外，自 Whittaker (1960) 將生物多樣性拆分成 α、β 與 γ 三層次之後，通過分解生物多樣性，也可以利用已知的 α 與 γ 多樣性來分解、轉換出 β 多樣性估計量。 群集間的多樣性：β diversity estimators β 多樣性估計量可以表示各個地區群集間的差異。與 α 多樣性類似，β 多樣性的評估大致上也分為只考慮物種數（richness）與同時考慮物種數與各物種個體充裕度（abundance）差異的估計量。

:::info :::spoiler {state="open"} 題次 ::: 在生物多樣性的估計量整理（上）：緒論與 α 多樣性一文中，提到期望物種數的概念，是為了比較不同的群集物種豐度而發展出來的。因為不同群集、或者不同調查趟次所得的總個體數不同，又或者努力量不同，就算估計出了有效物種數，也沒有辦法直接比較生物多樣性的高低，因此需要標準化之後，方可進行比較。 美國海洋無脊椎動物學家 Howard L. Sanders 在 1968 年一篇比較不同海底無脊椎動物群集豐度的文章中，提出了稀釋化（rarefaction）的作法[^1]，透過比較各樣點採獲的物種個體數與調查到物種百分比之間的關係，呈現其物種多樣性的情況。 其實在 Sanders 之前，已經有許多學者，企圖以計量方法呈現調查努力量（通常以調查到的個體數來代表）與調查到的物種數之間的關係。下表取自 Bobrowsky and Ball 1989，稍微調整代數符號：

image:https://i.imgur.com/lJLGp3h.png 在漁業資源調查中，體長是最容易得到的測量值之一。透過統計體長分布隨時序的變化，我們可以估算魚群（或其他經濟性捕撈動物）的成長速度變化、幼體補充時間等等。奧地利生物學家 Ludwig von Bertalanffy 於 1930 年代發展出了魚類的成長方程式 von Bertalanffy growth function（VBGF）如下： $$ 𝐿𝑡=𝐿_\infty [1−𝑒^{−𝐾(𝑡−𝑡_0)+(\frac{CK}{2𝜋})sin2𝜋(𝑡−𝑡_𝑠)}] $$ 其中 $L_t$ 是年齡為 $t$ 之體長；$L_\infty$ 為極限體長（asymptotic length），也是 VBGF 要逼近的目標；$t_0$ 為 $L_t=0$ 的理論年齡，一般在推論時，會假設 $t_0=0$；$K$ 為成長參數（growth coefficient）；$e$ 是自然指數；$C$ 為季節性振蕩幅度參數（seasonal oscillation），與年溫差成正比，溫差 1 度相當於變化 0.1 單位，其值介於 0（無振盪）到 1（最大振盪）之間；$t_s$ 為成長曲線起伏變化的起點，其值視生物棲息地的冬季低迷點（winter point，WP）而定，把一年 12 個月換算為 0 到 1 之間，生物一年中成長最慢的時段即為其 WP 所在。在北半球一般為 2 月中旬，即 WP=0.2。 因為 VBGF 估算起來相當複雜，因此聯合國糧食及農業組織（Food and Agriculture Organization）與世界魚類研究中心（International Center for Living Aquatic Resources Management, ICLARM）聯合研發了一套電腦體長頻度分析系統 ELEFAN（Electronic Length Frequency Analysis）提供全球漁業資源估算者來使用，可以搭配一樣是 FAO 推出的軟體 FishStat 系列（現在在他們網頁上的是 FishStatJ）。類似的系統還有 MULTIFAN，搭配 MULTIFAN-CL 軟體使用。 當初為了估算計畫需要的貝類 VBGF，本來想用 MULTIFAN-CL 來做統計，沒想到註冊了之後遲遲沒收到回信，感覺可能沒人維護了，因此想找找在 R 裡有沒有替代的套件。結果還真的找到了一個 TropFishR，由 Taylor 與 Mildenberger 2017 年發表於 Fisheries Management and Ecology 期刊。

重灌系統、裝好 R 之後，總是只能憑記憶安裝套件，而且常常會在需要的時候才發現少裝了某套件，就要停下手邊的工作來安裝，多少會打亂工作節奏。因此，想把平常使用的套件記錄下來，之後如果又重新裝了 R 或者整個系統，就可以一次裝好。 在開始裝之前，我會先預設 CRAN mirror，在 R 資料夾的 etc 目錄底下，修改 Profile.site 檔案： *[CRAN]:Comprehensive R Archive Network # set a CRAN mirror #Automatic redirection to servers worldwide, currently sponsored by Rstudio local({r = getOption("repos") r["CRAN"] = "https://cran.csie.ntu.edu.tw/" options(repos=r)})

生物多樣性是地球上生物體的種類與其差異（variety and variability）。國高中的生物課本中，我們學到生物多樣性包含遺傳多樣性、物種多樣性與生態系多樣性。 在生態調查報告中，需要呈現調查範圍內的生物多樣性，因為無法帶領每個閱讀報告的人到現場走一遭，需要透過一些估計量（estimators）方便讀者了解該地生物多樣性的程度；若調查方法一致、努力量相近，研究者也得以透過特定的生物多樣性估計量，比較該地歷史生物多樣性的變化，或者比較特定時間相近地點間的生物多樣性狀況。 當然，生態調查報告最為關注物種多樣性，因此這邊依據清大趙蓮菊老師的〈仰觀宇宙之大, 俯察品類之盛: 如何量化生物多樣性〉與中興大學馮豐隆老師介紹的〈生態歧異度及其求算方法之分類〉，並加入趙老師最新的研究成果與其他常用的生物多樣性估計量，做為自己的閱讀筆記。 :::info :::spoiler {state="open"} 題次 :::

高中生物課本告訴我們，若不考慮環境限制因子，一個生物族群的規模隨著時間會呈現指數式成長；而現實情況中因為資源有限，所以實際上的族群量會有個上限，那個上線就是環境負載量。 大學生態學課本裡，會用數學式子來表示。對於沒有任何限制，且大小為 $N$ 的族群，在 $t$ 時間區段的增長速度 $\frac{\mathrm{d}N}{\mathrm{d}t}=rN$，可以發現平均速度（per capita rate of increase）$r$ 為定值，跟族群數量無關，所以隨著時間過去，族群數量會直直衝，在族群大小—時間圖上形成 J 形曲線，表示為 $N_t=N_0\times(1+r)^t$。。 不過現實上，因為有環境負載量，所以實際上的量測的族群成長速率會是$\frac{\mathrm{d}N}{\mathrm{d}t}=r\left(\frac{K-N}{N}\right)N$，比指數成長曲線的速率多了$\left(\frac{K-N}{N}\right)$一項，其中 $K$ 就是環境承載量（carrying capacity）。當族群數量與環境承載量相當時，族群成長速率趨近於 0。考慮環境負載量的族群成長曲線會近似於邏輯式曲線（logistic curve），通式為：$N_t=\frac{N_0Ke^{rt}}{N_0(e^{rt}-1)+K}=\frac{K}{1+\left(\frac{K-P_0}{N_0}\right)e^{-rt}}$；隨時間過去，族群量會越來越近似於環境負載量：$\lim\limits_{t \to \inf}N_t=K$，最後在族群大小—時間圖上形成 S 形曲線。 取自BioNinjia 如何配適不同種類的族群成長曲線 如果手上有不同時間調查的族群量資料，要怎麼配適（fit）不同的成長曲線呢？ 如果要配適指數曲線，那很容易，在 R 當中只需要先使用 log() 函數，將時間與族群數量都取自然對數，那麼原本的指數成長函式，就會成為一個 $y=b+ax$ 形式的二元一次函式，此時使用 lm() 函數，就可以做簡單線性迴歸了。

在用 R 配適邏輯式族群生長曲線文章中，最後輸出了族群成長的邏輯式函數配適結果；如果可以在圖中加上函數的方程式，那就可以傳達更多資訊了。於是在這裡嘗試了幾個作法，其中需要用到 ggpmisc 套件，也可以在 ggplot 中手動輸出方程式。 使用ggpmisc把方程式加在圖上 library(ggplot2) library(ggpmisc) #原先的圖 p <- ggplot(data = data, aes(x = Year, y = Population)) + geom_line(data = predict, aes(x = x, y = y), size = 1) + geom_point(color = "blue",size = 4) + labs(x = "Year", y = "Estimated population size") +

閱讀 metabarcoding（高通量分子條碼）的研究時，因為會定序出多量短片段的序列需要分析，常會看到 OTU 或者 ASV 這類名詞。 OTU 是操作分類單元（operational taxonomic units），而 ASV 是擴增子序列變異（amplicon sequence variant）。OTU 的概念很早以前就出現[^2]，而 ASV 的概念一直到 2013 年才被提出[^1]；雖然中間或後續也有許多名詞如 ZOTU(zero-radius OTU)、ESV (Exact Sequence Variant)、sOTU (sub-OTU)、ISU (individual sequence variant)，不過概念均與 ASV 相當類似，因此下文就主要比較 OTU 與 ASV 的差別。 兩者的不同 OTU 與 ASV 的差別在演算的方式，前者用的是聚類（clustering），後者用的是序列校正或者去噪（denoising）。 OTU 的產生，主要是據序列的相似性，以一定的相似度把序列歸集到一起。一般用的相似度是 97% 或 99%，產生出來的 OTU 就叫 OTU97 或 OTU99。 以 NGS 產出的資料來說，這個方法的缺點主要是敏感度不足，一個 OTU 中可能包含有不同的實際物種；再者，如果序列因為定序過程產生錯誤，會造成不可信賴的聚類結果。目前比較多人用以產出 OTU 的工具，是 USEARCH 中的 UPARSE 演算法。USEARCH 本身為付費軟體，不過有開放 32 位元版免費使用。

之前用R做過氣候圖，這次換個視覺化的方法，繪製全年的每日氣溫分布圖。從圖中可以快速了解一年當中每天的高、低溫與平均溫度。 一開始，當然要從氣象署的觀測資料查詢系統抓出每天的氣溫資料。以下用已撤銷的香山測站為例，參考阿好伯的做法，抓出設站開始到撤銷之前完整年的資料： library(dplyr) library(jsonlite) library(rvest) library(magrittr) library(lubridate) library(ggplot2)

:::info :::spoiler {state="open"} 題次 ::: 這篇稍微介紹較少人使用的親緣關係多樣性（phylogenetic diversity，原譯系統演化多樣性，或譯譜系多樣性）及功能多樣性（functional diversity），最後介紹趙蓮菊老師發展出的屬性多樣性指數（attribute diversity index）如何統一估計這一系列文章提到的所有多樣性。 親緣關係多樣性 親緣關係多樣性是一種帶有物種差異思考的多樣性估計量，一開始是 1970 年代由 Pielou 所提出[^4]。在計算生物多樣性時，若納入物種間的差異，做出的估計稱為親緣多樣性（系統演化多樣性、譜系多樣性）。

NCBI 中有 sequence identifiers 作為序列的唯一識別編號，可以透過 Entrez Molecular Sequence Database System 進行序列編號的大量操作。NCBI 有推出 ncbi-entrez-direct 小工具，方便我們使用 Entrez 的服務。以下以 NCBI Accession Number "OL800728" 為例，進行資料的查找： $ esearch -db nucleotide -query "OL800728"| esummary 輸出以下結果： <?xml version="1.0" encoding="UTF-8" ?> <!DOCTYPE DocumentSummarySet PUBLIC "-//NLM//DTD esummary nuccore 20201205//EN" "https://eutils.ncbi.nlm.nih.gov/eutils/dtd/20201205/esummary_nuccore.dtd"> <DocumentSummarySet status="OK"> <DocumentSummary> <Id>2259973435</Id>

Julia 於 2012 年面市，當初的開發者希望可以開發出一款在效率與易用性之間取得平衡的一款程式。Julia 允許並行（concurrent）、平行（parallel）與分散（distributed）計算，是一款非常強大、高效和靈活的程式語言，漸漸被人們使用於各種科學計算、生醫領域資料統計與數值分析等任務。 這邊紀錄安裝 julia 的步驟，希望未來能新增幾篇使用紀錄與心得。首先到 julia 官方網站確認目前最新穩定版本號，接著就可以安裝了。 以 1.8.5 版為例，可以從網站下載編譯好的軟體直接用，而我選擇下載原始碼自行編譯，過程如下： # 首先安裝相依性軟體 sudo apt update sudo apt install build-essential libatomic1 python gfortran perl wget m4 cmake pkg-config curl # 接著切到 /tmp 資料夾準備編譯

最近的工作需要用到台灣物種名錄（Catalogue of Life in Taiwan，就是以前的臺灣生物多樣性國家資訊網 Taiwan Biodiversity National Information Network, TaiBNET，為了和臺灣生物多樣性國家入口網站 Taiwan Biodiversity Information Facility, TaiBIF 區隔而改中文名）整合的物種分類階層、原生種、保育類、紅皮書等資料來整理每次生態調查到的資料。一樣的事情需要做超過 3 次的話，還是考慮搭配程式來做比較省時，因此就試著用 bash 裡的關聯陣列（associative array）功能來做。 在 bash 4.2 版之後，就有了 associative array 功能。現在大部分的 ubuntu 中都在 4.2 版之後，若要確認，使用 $ bash --version 確定版本後，就可以準備開始製作 associative array。 定義方式有以下 2 種： $ declare -A assArray1 $ assArray1[dog]=woof