2018年8月22日水曜日

生態系の熱帯化:藻場が狭まり、サンゴ群集が拡大する

渾身の論文がPNAS(米国科学アカデミー紀要)出版になりました(プレスリリースも出ました):
Naoki H. Kumagai, Jorge García Molinos, Hiroya Yamano, Shintaro Takao, Masahiko Fujii, Yasuhiro Yamanaka (in press) Ocean currents and herbivory drive macroalgae-to-coral community shift under climate warming. Proceedings of the National Academy of Sciences of the United States of America.
(https://doi.org/10.1073/pnas.1716826115)

研究内容や説明についてはプレスリリースに譲りますが、以下は書き切れなかった、もう少し個人的な部分や苦労についてメモしておきます。

研究開始から約5年、初投稿から1年9ヶ月も掛かった大作です(Nは1ヶ月でエディターキック、Sの査読まで進んだものの手が届かず、PNASで2回の改訂ののち受理)。これまでに培った、移動分散、メタ個体群の考え方、データベースの構築、多くの観察データを集めて用いる統計モデリング技術、ベイズ推定、GISなど、全てを込めました。

この論文では、主に1950〜2015年の、数多くの地道な観察記録(439文献、22,253調査記録)をしらみつぶしに探索し、観察年・位置情報を付加した上で、45種の生物種ごとの日本国内の分布変化を網羅しました。私はおもに海藻と魚類の分布変化記録の収集と整備を担当しました(サンゴも一部だけ)。興味のある方はSupporting InformationのFig S2を見ていただけると45種の個々の分布変化が分かります。
お陰で日本のあらゆる海岸線の地名には詳しくなり、市町村合併や海岸線の変化(埋め立て)の歴史も垣間見ました。

この論文ではさらに、藻場を構成する海藻や、より南方から分布を拡げてくる造礁サンゴ、海藻を食害する魚類の分布変化速度を指標として、藻場とサンゴ群集の分布変化を解析しています。気候変動影響のもとで、海流の輸送といった大スケールの物理的要因および魚類による海藻の食害といった生物間相互作用が組み合わさることで、観察された分布変化をモデルによってうまく再現できることを示しました。これらの解析・モデルを最終的な形態にもっていく過程では、共著者の皆さまに大きく助けていただきました。この研究を通じて、私自身のモデリング技術、GIS技術とその理解も大きく向上したので、実によい勉強機会になったと思います。

PNASなどのshort formatスタイルは限られた字数の中に、必要なことを凝縮しつつ、しかも読みやすくしなければならないという、実際やってみると非常に難しい作業でした。字数が短い=簡単と考える人も居るかもしれませんが、記述的な新発見の論文でもない限り、short formatの論文の執筆はフルペーパーを書くよりはるかに難しい(あるいは慣れていないと難しい)作業だと感じました。これについても論文を改訂する過程で共著者の皆さまには大きく助けていただきました。自分自身の文章力も向上できたと思います。

この論文で述べたように、今後の日本の温帯藻場は温暖化の直接的影響だけでなく、魚類やウニなどの摂食圧の影響がいっそう深刻になると予想されます。またサンゴの分布拡大も気候の変化速度よりも遅いため、海藻藻場もサンゴ群集も何かしら人の手を加えて保全する手立てを考えていく必要があると思います。この論文がそういう議論の必要性を感じる1つの切っ掛けになればと思います。

なお、実際の海中では同じ海(例えば1回のダイビングの中で共存するくらいの範囲)に藻場とサンゴ群集が居合わせることはよくあるけれど、いざこれを写真に収めようとすると適した場所はなかなかありません。査読2回目の改訂の際に偶然、この論文で取り上げた全生物要素がバランスよく生息するシンボル的な海(宇和島市津島町田之浜)を訪れることができました。論文のFig.1A中央、プレスリリースの写真1中央の写真は、南方性コンブ類、温帯性ホンダワラ類、南方性ホンダワラ類、南方性サンゴの全要素が凝縮された奇跡の一枚です。
(参照:田之浜の調査でお世話になった、愛媛ダイビングセンターさんのブログでも取り上げていただきました)

2018年5月29日火曜日

Rいろは:グラフ作成ggplot2編

###### ggplot2とは? #################
# 標準仕様でもそれらしい科学的グラフを作ってくれるRパッケージ
# (plot関数の場合、相当なアレンジをしないと使い物にならない)
# 色分けや複合グラフ、重ね合わせなどに強い、応用の利きやすいメリットがある
# さらにcowplotパッケージを適用して、より学術用に適したグラフ作成を紹介します

# 項や足し算形式でグラフ要素を追加する形態になっている。例えば、散布図、棒グラフのようなグラフ形式や、回帰直線など。
# 項の指定方法などが(plot関数と違って)英単語を元にしているものが多く連想しやすい

# 使い方の早見表も参考になります

install.packages("ggplot2", dep=T) # 未インストールの場合
install.packages("cowplot", dep=T) # 未インストールの場合
install.packages("reshape", dep=T) # 未インストールの場合
require(ggplot2) # 呼び出し

###### まず定番のdata"iris"を用いて、originalのggplot2で描いてみます #########

data(iris) # iris呼び出し
head(iris) # irisの構成をチェック
  #Sepal.Length Sepal.Width Petal.Length Petal.Width Species
#1          5.1         3.5          1.4         0.2  setosa
#2          4.9         3.0          1.4         0.2  setosa
#3          4.7         3.2          1.3         0.2  setosa
#4          4.6         3.1          1.5         0.2  setosa
#5          5.0         3.6          1.4         0.2  setosa
#6          5.4         3.9          1.7         0.4  setosa



##### 散布図 ####################
ggplot(data=iris, aes(x=Sepal.Length, y=Sepal.Width, fill=Species)) + geom_point(pch=21, size=2)

# ggplot(data=iris, # データ元の指定
#  aes(x=Sepal.Length, y=Sepal.Width, # x, yに用いる変数の指定(aes()内で変数を指定)
#      fill=Species)) # Speciesによって異なる色を使用
# + geom_point( # + geom_xxx: xxxの部分でグラフ種類を指定
#  pch=21,  # pch=21: 塗りと枠に異なる色を適用可能(21~25) # 全種に適用する場合、aes外に記述
#  size=2) # プロットサイズ(少し大きくしました)

# 以下でも同様のグラフになります(fill=Speciesを後半に記述)
ggplot(data=iris, aes(x=Sepal.Length, y=Sepal.Width)) + geom_point(aes(fill=Species), pch=21, size=2)

# このままでは背景がうるさく、学術用に不向き
require(cowplot) # cowplotパッケージを呼び出し、同じようにグラフを作成します
theme_set(theme_cowplot()) # このコードを通します(すると以降、cowplotスタイルに変更される)

ggplot(data=iris, aes(x=Sepal.Length, y=Sepal.Width, fill=Species)) + geom_point(pch=21, size=2)
# 背景が白地のスッキリしたものに変わりました。文字のサイズもちょうどいいです。



##### 箱ひげ図 ####################
ggplot(data=iris, aes(x=Species, y=Sepal.Width, fill=Species)) + geom_boxplot()



##### 棒グラフ+SD ####################
# サンプルデータの準備
MEAN1 <- with(iris, tapply(Sepal.Width, Species, mean)) # 各種の平均("対象データ", "グループ分け", "関数名")
SD1 <- with(iris, tapply(Sepal.Width, Species, sd)) # 各グループの標準偏差
iris2 <- data.frame(Species=names(MEAN1), me=MEAN1, sd=SD1)
# グラフの作成
ggplot(data=iris2, aes(x=Species, y=me, fill=Species)) + geom_bar(stat="identity") + geom_linerange(aes(ymax=me + sd, ymin=me - sd))



##### 2要因の棒グラフ + SD ####################
# サンプルデータの準備
iris3 <- rbind( # 1.2倍したサンプルデータを追加
data.frame(Species=names(MEAN1), me=MEAN1, sd=SD1, treat="C"),
data.frame(Species=names(MEAN1), me=1.2*MEAN1, sd=1.2*SD1, treat="T") )
# グラフの作成
ggplot(data=iris3, aes(x=Species, y=me, fill=treat)) + geom_bar(stat="identity", position="dodge") + geom_linerange(aes(ymax=me + sd, ymin=me - sd), position=position_dodge(width=0.9))
# position, dodge: 接したグラフを作成するためのkeyword



##### ヒストグラム ####################
ggplot(iris, aes(x=Sepal.Width, fill=Species)) + geom_histogram()



##### 時系列データのグラフ ####################
data(airquality) # data"airquality"を使用
head(airquality) # データの中身をチェック
  #Ozone Solar.R Wind Temp Month Day
#1    41     190  7.4   67     5   1
#2    36     118  8.0   72     5   2
#3    12     149 12.6   74     5   3
#4    18     313 11.5   62     5   4
#5    NA      NA 14.3   56     5   5
#6    28      NA 14.9   66     5   6

# 時間形式に変換
airquality$Time <- as.POSIXct(strptime(with(airquality, paste(2018, Month, Day, sep="-")), "%Y-%m-%d")) # (2020.03.03少しコードを修正)
ggplot(airquality, aes(x=Time, y=Ozone)) + geom_point(pch=21, cex=2, fill="purple") + geom_line()



##### 時系列データのグラフ(複数系列) ####################
require(reshape) # 要インストール
airquality2 <- melt( # データ構成を変換
airquality[,c("Ozone","Temp","Time")], id.var="Time")
head(airquality2)
        #Time variable value
#1 2018-05-01    Ozone    41
#2 2018-05-02    Ozone    36
#3 2018-05-03    Ozone    12
#4 2018-05-04    Ozone    18
#5 2018-05-05    Ozone    NA
#6 2018-05-06    Ozone    28

# 2パネル構成
base <- ggplot(airquality2, aes(x=Time, y=value)) # ごちゃごちゃしてきたので分割します
poi <- geom_point(pch=21, cex=2, aes(fill=variable)) 
lin <- geom_line(aes(colour=variable)) 
base + poi + lin + facet_wrap(~ variable, scale="free_y", ncol=1) # variableでパネルを分割




###### 散布図+回帰直線 ##################
summary(iris) # 値の幅をチェック
  #Sepal.Length    Sepal.Width     Petal.Length    Petal.Width          Species
 #Min.   :4.300   Min.   :2.000   Min.   :1.000   Min.   :0.100   setosa    :50
 #1st Qu.:5.100   1st Qu.:2.800   1st Qu.:1.600   1st Qu.:0.300   versicolor:50
 #Median :5.800   Median :3.000   Median :4.350   Median :1.300   virginica :50
 #Mean   :5.843   Mean   :3.057   Mean   :3.758   Mean   :1.199                
 #3rd Qu.:6.400   3rd Qu.:3.300   3rd Qu.:5.100   3rd Qu.:1.800                
 #Max.   :7.900   Max.   :4.400   Max.   :6.900   Max.   :2.500                

md1 <- lm(Sepal.Length ~ Sepal.Width, iris[iris$Species=="setosa",]) # 1種に限定
coef(md1)
#(Intercept) Sepal.Width
  #6.5262226  -0.2233611

# 散布図 + 回帰直線
x.interval <- seq(2, 5, 0.1) # x(Sepal.Width)の最小〜最大を0.1刻みで細かく分割
est.Y <- predict(md1, newdata=data.frame(Sepal.Width=x.interval))
fitted <- data.frame(Sepal.Width=x.interval, Sepal.Length=est.Y)
# fitted <- data.frame(Sepal.Width=x.interval, Sepal.Length=coef(md1)[1] + coef(md1)[2]*x.interval) # 上と同じ
base2 <- ggplot(iris[iris$Species=="setosa",], aes(x=Sepal.Width, y=Sepal.Length)) 
poi2 <- geom_point(fill="turquoise", pch=21, size=2) 
lin2 <- geom_line(data=fitted, colour="turquoise", lwd=1) # ここは"data="を明記する必要あり
base2 + poi2 + lin2



##### 散布図 + 回帰直線 + 信頼区間 ##################
est.Y2 <- predict(md1, newdata=data.frame(Sepal.Width=x.interval), interval="confidence")
head(est.Y2)
       #fit      lwr      upr
#1 4.019981 3.753101 4.286860
#2 4.089030 3.839589 4.338470
#3 4.158079 3.925981 4.390177
#4 4.227128 4.012251 4.442004
#5 4.296177 4.098369 4.493984
#6 4.365226 4.184291 4.546160

fitted2 <- data.frame(Sepal.Width=x.interval, Sepal.Length=est.Y2[,1], upperCI=est.Y2[,3], lowerCI=est.Y2[,2])
band <- geom_ribbon(data=fitted2, aes(ymin=lowerCI, ymax=upperCI), alpha=0.5, fill="turquoise", linetype="blank") # alpha: 半透明(0~1)
base2 + poi2 + band + lin2 # 重ねたい順で



##### 散布図 + 回帰直線(複数グループ) ##################
x.interval3 <- seq(1,5,0.1)
new.dat <- data.frame( # 種数分、繰り返す
Species=rep(unique(iris$Species),each=length(x.interval3)), 
Sepal.Width=rep(x.interval3, length(unique(iris$Species))))
md3 <- lm(Sepal.Length ~ Sepal.Width*Species, iris) # 全種を対象、Speciesによる違い(ANCOVA)
est.Y3 <- predict(md3, newdata=new.dat)
fitted3 <- data.frame(Species=new.dat$Species, Sepal.Width=new.dat$Sepal.Width, Sepal.Length=est.Y3) # Species列を追加
base3 <- ggplot(iris, aes(x=Sepal.Width, y=Sepal.Length)) 
poi3 <- geom_point(aes(fill=Species), pch=21, size=2) 
lin3 <- geom_line(data=fitted3, aes(colour=Species), lwd=1) # ここは"data="を明記する必要あり
base3 + poi3 + lin3




##### 複数タイプのグラフの組み合わせ(cowplot利用) ##################
g1 <- ggplot(data=iris, aes(x=Sepal.Length, y=Sepal.Width)) + geom_point(aes(fill=Species), pch=21, size=2) 
g2 <- ggplot(data=iris, aes(x=Species, y=Petal.Width)) + geom_boxplot(aes(fill=Species)) # 箱ひげ図
+ theme(axis.text.x = element_text(angle=90, vjust=0.5)) # label方向

g12 <- plot_grid(g1, g2, labels="auto", align="h") 
# labels="auto":小文字でラベリング(大文字は"AUTO")
# 直接入力の例:labels=c("a","b")
# align="h":図の端を揃える("h", "v", "hv")

#### グラフの調整(例:凡例を統一、サイズ調整)
g1n <- g1 + theme(legend.position="none") # 一方の凡例消去
bd <- panel_border(colour=1, size=1) # 外枠追加
g1n2 <- plot_grid(g1n + bd, g2+ bd, labels="auto", 
   rel_widths=c(1, 1.4), # 横幅を1:1.4に変更
   scale=0.9, # パネル間の間隔を広げる
   vjust=0.2) # ラベル位置(a, b)の調整
 
 
 
#### ファイルの保存 ##############################
save_plot("directry/filename.pdf", g1n2) # cowplotの関数、お任せで適度なサイズになる。指定も可能(2020.03.03: ggsaveからsave_plotへ変更)

#### 水平線・垂線 ##############################
abl <- geom_hline(yintercept=0, linetype="dotted") # 水平 ("dotted": 点線)
abl2 <- geom_vline(xintercept=0, linetype="dotted") # 垂直


#### 軸とラベルの調整 ############################
x.ax <- scale_x_continuous(breaks=c(1:5), limits=c(1, 5), expand=c(0.001, 0.001), name="Sepal width (mm)")
# ex.
base3 + poi3 + lin3 + x.ax


#### 色の指定 ##############################
# グレースケール
set.fillG <- scale_fill_grey(start=0.5, end=1) # グレー 〜 白の場合(塗りつぶし色)
# "fill"の部分を"colour"にすると線の色の指定("shape": 形状、

set.fillC <- scale_fill_manual(values=c("turquoise", "salmon", "blue")) # 塗りつぶし色
set.colC <- scale_colour_manual(values=c("turquoise", "salmon", "blue")) # 線の色
# ex.
base3 + poi3 + lin3 + set.fillC + set.colC

#### cf. 色コードの選択 ##########################
"見やすい"色の組み合わせ例の参照サイト(ColorBrewer)
fill="#2ca25f" のようにしてカラーコードで色の指定が可能



#### 90° 回転(グラフを横向きにする)#######################
 + coord_flip()


#### 凡例の調整 ##############################
+ theme(legend.position="none") # 凡例消し
+ theme(legend.position="bottom") # ラベルの右揃え(cf. vjust)

# 長すぎるラベルを改行する( \n が使える)
scale_x_continuous(name = "Length\n (cm)")

# 上付き、ギリシャ文字など
expression(m^2) 
expression(Delta) # Δ



#### cowplotの調整 ##############################
+ panel_border(colour=1, size=1) # defaultの線・文字色がなぜかグレーなので黒に変更
+ theme_cowplot(font_size = 12) # 文字が大きく感じる場合はサイズ指定

# 日本語の利用(ただしアウトライン化されてしまう模様 )
ggsave(..., family="Japan1GothicBBB",...)