2020年12月20日日曜日

tidyverse版: 複数生物種の個体数集計データを縦に伸ばすor縮める(集計データ← →データフレーム変換)

 (以前の記事をtidyverseで作り直してみました。まだ旧来のRコードの方が馴染みがあるので、tidyverseに慣れている方からすると不格好かもしれません。)

多数の種が含まれる生物群集データを解析する時、エクセルなどにデータを整理していると、各種をインデックスにするか要因のひとつとして整理するか、ケースによって変わってくるかと思います。とくに各種を応答変数とするか説明変数に使用するかあたりで必要になってくるかと。完全なデータフレーム型にしておいた方が統計解析はやりやすいけれど、膨大な行数となってエクセルの限界にあっさりと到達したりもするし…(cf. エクセルのvlookup関数)。


こんな時に、一部集計データとデータフレーム型とを自在に一発変換できたら効率がよいと思い、やり方をまとめてみました。正直、入力の手間が一番省けるのは部分的には集計されたデータだったりしますよね(下記のdata1のような形)。

データの整形にはdplyrパッケージの関数などを使用します。tidyverseパッケージとして関連パッケージも一括でインストールが可能です。

require(tidyverse) # 要インストール、tidyverseのパッケージ群が一括で呼び出される。

# 例えば、こんな種毎に集計されているデータがある時
# (tidyverse標準のtibble形式で作成)
data1 <- tibble(
year = c(rep("y08", 3), rep("y09", 3)), 
site = c(1:3, 1:3), 
depth = seq(1, 5, length=6), 
sp1 = c(6:1), sp2 = c(1:6), sp3 = c(0:5))

# 出力するとこんな感じです。<chr>は文字型、<int>は整数型、<dbl>は連続変数型を示す。
> data1
# A tibble: 6 x 6
  year   site depth   sp1   sp2   sp3
  <chr> <int> <dbl> <int> <int> <int>
1 y08       1   1       6     1     0
2 y08       2   1.8     5     2     1
3 y08       3   2.6     4     3     2
4 y09       1   3.4     3     4     3
5 y09       2   4.2     2     5     4
6 y09       3   5       1     6     5



########################################################

# data1をデータフレーム型に縦に伸ばす場合
data2 <- data1 %>% pivot_longer(c(-year, -site, -depth), names_to="species", values_to="abundance")
# pivot_longerは集計された種名の表を縦長に伸ばす。種名を表す新たな種をspecies, 個体数をabundanceとする。
# 形式を保持したい変数は、c(-year, -site, -depth)のように記す(c()は必要)。ひとつだけならば-yearのように書けば良いが、複数の場合は-c(year, site, depth)ではないことに注意。一般的なRの記法と異なっており、厄介な点。
# "%>%"はパイプと呼ばれる。パイプの前のデータを後の処理に渡す意味。いちいちdata1$yearのように書かなくていい、with関数やattach関数と似ている。可読性は良くなっているかもだけれど、デバッグがしづらい気もする。

> data2 # 出力してみます
# A tibble: 18 x 5
   year   site depth species abundance
   <chr> <int> <dbl> <chr>       <int>
 1 y08       1   1   sp1             6
 2 y08       1   1   sp2             1
 3 y08       1   1   sp3             0
 4 y08       2   1.8 sp1             5
 5 y08       2   1.8 sp2             2
 6 y08       2   1.8 sp3             1
 7 y08       3   2.6 sp1             4
 8 y08       3   2.6 sp2             3
 9 y08       3   2.6 sp3             2
10 y09       1   3.4 sp1             3
11 y09       1   3.4 sp2             4
12 y09       1   3.4 sp3             3
13 y09       2   4.2 sp1             2
14 y09       2   4.2 sp2             5
15 y09       2   4.2 sp3             4
16 y09       3   5   sp1             1
17 y09       3   5   sp2             6
18 y09       3   5   sp3             5


########################################################
# もう一段階、伸ばしてみましょう。一個体あたり一行というデータセットへ。多項ロジットなど、カテゴリカルな解析に便利そうです。

data3 <- select(data2, -abundance) %>% slice(unlist(map2(1:nrow(data2), data2$abundance, rep)))
# data2のabundanceの数に応じてrepで行数を増やし、sliceで行を選択。map2はおおよそmapplyに相当、2変数を取るsapply。slice(unlist(map2(,,,,)))はもう少しスッキリ書けないものだろうか。

> data3 
# A tibble: 57 x 4
   year   site depth species
   <chr> <int> <dbl> <chr>  
 1 y08       1   1   sp1    
 2 y08       1   1   sp1    
 3 y08       1   1   sp1    
 4 y08       1   1   sp1    
 5 y08       1   1   sp1    
 6 y08       1   1   sp1    
 7 y08       1   1   sp2    
 8 y08       2   1.8 sp1    
 9 y08       2   1.8 sp1    
10 y08       2   1.8 sp1    
# … with 47 more rows




########################################################

# data3を一つ前の状態に戻してみます(観測地点あたりで集計、といったところ)
data4 <- data3 %>% count(year, site, depth, species, name="abundance") 
# year, site, depth, species毎にデータ頻度をカウントし、カウント結果をabundance変数として追加

> data4

# A tibble: 17 x 5
   year   site depth species abundance
   <chr> <int> <dbl> <chr>       <int>
 1 y08       1   1   sp1             6
 2 y08       1   1   sp2             1
 3 y08       2   1.8 sp1             5
 4 y08       2   1.8 sp2             2
 5 y08       2   1.8 sp3             1
 6 y08       3   2.6 sp1             4
 7 y08       3   2.6 sp2             3
 8 y08       3   2.6 sp3             2
 9 y09       1   3.4 sp1             3
10 y09       1   3.4 sp2             4
11 y09       1   3.4 sp3             3
12 y09       2   4.2 sp1             2
13 y09       2   4.2 sp2             5
14 y09       2   4.2 sp3             4
15 y09       3   5   sp1             1
16 y09       3   5   sp2             6
17 y09       3   5   sp3             5

# sp3 のゼロ個体データだけは消えてしまったが、当然と言えば当然




########################################################

# data4の種数部分を集計表状に横に伸ばす(最初の状態=data1の形にする場合)

data5 <- data4 %>% pivot_wider(names_from=species, values_from=abundance, values_fill=0)
 

# pivot_widerは縦長データから横長データへ変換する、見出しとなるspeciesと対応する値abundanceを展開する。生じるNAはvalues_fill=0で置き換えた。

> data5 # 出力してみます。
# A tibble: 6 x 6
  year   site depth   sp1   sp2   sp3
  <chr> <int> <dbl> <int> <int> <int>
1 y08       1   1       6     1     0
2 y08       2   1.8     5     2     1
3 y08       3   2.6     4     3     2
4 y09       1   3.4     3     4     3
5 y09       2   4.2     2     5     4
6 y09       3   5       1     6     5



########################################################
# 出現する種のリストにひどく変動がある場合など、1つのセルの中に種名を羅列したくなるケースもあると思います。そういうデータからの変換例も追加しておきます。

data6 <- tibble(site=paste0("s", 1:10), month=rep(c(1:2), each=5), 
species=c("スズメ", "","スズメ、ヒヨドリ、シジュウカラ", "ムクドリ、スズメ", "ヒヨドリ、スズメ", "スズメ、キジバト", "スズメ、ムクドリ", "スズメ、ヒヨドリ", "", "ムクドリ"))

# 例えばこんなかんじのデータセット。空欄さえあります…
# A tibble: 10 x 3
   site  month species                         
   <chr> <int> <chr>                           
 1 s1        1 "スズメ"                        
 2 s2        1 ""                              
 3 s3        1 "スズメ、ヒヨドリ、シジュウカラ"
 4 s4        1 "ムクドリ、スズメ"              
 5 s5        1 "ヒヨドリ、スズメ"              
 6 s6        2 "スズメ、キジバト"              
 7 s7        2 "スズメ、ムクドリ"              
 8 s8        2 "スズメ、ヒヨドリ"              
 9 s9        2 ""                              
10 s10       2 "ムクドリ"                      

# 注:このような日本語データのでの入力の場合、普通に読み込もうとすると文字化けしやすいので、readxlパッケージのread_excel関数などで読み込むのがトラブルが少なくていいと思っています。csvでは機種依存文字やmac/windows間での文字コード問題が混乱のもとなので。

# こんなデータでも次のようにすれば、通常の集計データに変換できます。data1の形態へ変換してみましょう。

sptabs <- data6 %>% # data6をベースに集計作業
pull(species) %>% # 種名部分をベクトルとして抽出
map(str_split, pattern="、") %>% # "、"で文字列を分割
map_dfr(table) %>% # 要素毎に集計テーブルを作成
map_dfr(replace_na, 0) %>% # naを0で置き換え
select(-starts_with("...")) # 出現無しを削除
data7 <- select(data6, -species) %>% # 元データの見出しを抽出
mutate(sptabs) # 種集計テーブルと結合

data7 # 出力してみます(出来上がり)

# A tibble: 10 x 7
   site  month スズメ  シジュウカラ ヒヨドリ ムクドリ キジバト
   <chr> <int> <table> <table>      <table>  <table>  <table> 
 1 s1        1          1       0            0        0        0       
 2 s2        1          0       0            0        0        0       
 3 s3        1          1       1            1        0        0       
 4 s4        1          1       0            0        1        0       
 5 s5        1          1       0            1        0        0       
 6 s6        2          1       0            0        0        1       
 7 s7        2          1       0            0        1        0       
 8 s8        2          1       0            1        0        0       
 9 s9        2          0       0            0        0        0       
10 s10       2          0       0            0        1        0       

2020年9月3日木曜日

Zoom / Teamsでスライドショーにならない/スライドが進まない問題の解決:macで外部ディスプレイ使用時?

コロナ禍のもと、オンラインツールを用いてプレゼンをする機会は多いと思います。Zoom / Teamsでスライドショーにならない/スライドが進まないことがあり、探しても解決策が見つからなかったため(書いてあっても単に私が理解していなかっただけかもですが)ここで紹介しておきます。いずれも外部ディスプレイの使用時のトラブル例です。macでと書いていますが、Windowsで同様のことが起こるのかどうか単に未確認です。


Zoom、Teamsの両方に共通の解決方法が分かってきたので、その手順を先に書いておきます:

0. とっさの場合には、外部ディスプレイを引っこ抜いてmac本体単体で共有をやり直すのが手っ取り早いです。以下は、外部ディスプレイを保持したまま正常に共有でスライドショーにするための対処法です。

1. まずPowerPointでスライドショーを開始します。するとメイン画面は発表者ツール、外部ディスプレイが全画面のプレゼンになっているはずです。

2. PowerPointの発表者ツールの左上に、ツールのボタンが並んでいますが、「ディスプレイの入れ替え」をクリックします。

すると、メイン画面が全画面のプレゼン、外部ディスプレイが発表者ツールに入れ替わるはずです。



3. この状態でZoom / Teamsへ移動し、通常通りにファイル共有ボタンをクリックします。この状態ではPowerPointのファイルを選択する形ではなく、ウィンドウ1または2の選択する形で全画面プレゼンになっているメイン画面を選択します。なお、PowerPointですでにプレゼンモードになっている状態でZoom / Teamsへ移動するには、command + タブ(このマークのボタン ->| )を押すと移動先のアプリを選択できるようになります。

4. これで相手側にも正常にZoom / Teams上でスライドショーが表示されているはずです。

5. プレゼンが終了したら、先に共有を終了してから、次にスライドショーを終了することをお勧めします。先にスライドショーを終了すると、デスクトップをさらけ出すことになります。

6. 一度この手順を踏んでおくと、その後、同じ環境では設定が保存されているようです。別のファイルを開こうとした際にも、初めから外部ディスプレイに発表者ツールが現れました。ただし、いざリアルのプロジェクターを用いたプレゼンに戻る際には注意が必要そうです。


参考までに、自分自身ではどういう現象が起きたのかを書いておきます。

●macのZoomファイル共有でPowerPointプレゼンのスライドショーが進まない場合

手元では共有したPowerPointをスライドショーに切り替え、手元ではスライドショーが進んでいるのに、相手側からは"最初のページから進んでいない"と言われました。後で実験してみると、相手側ではスライドショーモードにすら切り替わっていないことが分かりました。スライドショーを切って、編集モードの状態でスライドを進めると、相手側でも進みます。


●macのTeamsファイル共有でPowerPointプレゼン全画面できない場合 

手元では共有したPowerPointをスライドショーに切り替え、手元ではスライドショーが進んでいるのに、相手側からは"編集モードのままだ"と言われました。ディスプレイ設定をミラーリング、拡張のいずれもダメ、PowerPoint側を全画面/プレゼンツールへ切り替えなど試しましたが、自分の手元でしか全画面プレゼンにならずに困りました。

2020年8月30日日曜日

{cmdstanr}: Stan高速コンパイル、{rstan}代替としても有用

前回、最近mac版のRパッケージ、とくにコンパイル必要パッケージのインストールが困難になっているという記事を書いた。R3.6ではひどい状況だったが、R4.0.2の現在、CRANの手順に従えば、おおよそ問題はなくなっている。残る大きな問題は、rstanの動作不良だ。手元のR4.0.2では、rstanのインストール自体はできても、モデルのコンパイルができなくなっている。Macだけでなく、最近Windowsでもrstanの問題が出ている模様。

対策方法を調べていったところ、cmdstanrというコマンドラインベース(コンソール?macならTerminal)のRパッケージに辿り着いた。rstanのようにいちいちRcppでコンパイルしないので動作が速いし、またバージョンアップの際の依存環境との相性問題も生じにくいようだ。mac, linux, windowsのいずれでも使用可能なようだ。コマンドラインベースとは言うものの、Rから全ての作業を完結できることが分かったので、使い勝手も悪くない。

cmdstanrの公式ページが分かりやすいので、とくに説明も要らなさそうだが、実行サンプルを載せておく。

# インストール: まず以下をR上で実行(公式ページの記述そのまま)

install.packages("cmdstanr", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))

# または、開発版のインストール(上記ではコンパイル失敗するモデルがあったが、そういう場合でも開発版では成功している→2020.09.17追記)

devtools::install_github("stan-dev/cmdstanr")

# 次に、以下でcmdstanrを呼び出し、

require(cmdstanr)

# R上で、下記のインストール関数を実行(coresの数値は、PCのコアをいくつ使用するかの指定、そこまで時間は掛からない)

install_cmdstan(cores = 2)

cmdstan_path() # インストールされるディレクトリを示す。

# 2回目以降の使用で謎エラーが出てモデルがコンパイルできない場合、いったんcmdstanrを削除しインストールし直すことで使えている。いい対処法ではないだろうけれど、今のところ他にうまい手が見つからない。


####################################

# 計算実行:まず、モデル式を記述(公式サンプルではファイルからの読み込み)、

stan_program <- "

data {

  int<lower=0> N;

  int<lower=0,upper=1> y[N];

}

parameters {

  real<lower=0,upper=1> theta;

}

model {

  y ~ bernoulli(theta);

}

"


modfile <- write_stan_file(stan_program)

mod <- cmdstan_model(modfile) # 一時ファイルへ書き出し、この方が運用しやすいと思う

#( threads_per_chainを用いる場合(並列化)は以下のoptionが必要)

mod <- cmdstan_model(wt_code, cpp_options = list(stan_threads = TRUE)) 


# あとはrstanと同様の計算実行(公式の計算例、そのまま)

data_list <- list(N = 10, y = c(0,1,0,0,0,0,0,0,0,1))


fit <- mod$sample(

  data = data_list,

  seed = 123, # 乱数の種

  chains = 4,

  parallel_chains = 2, # 並列化するchainの数

  refresh = 500 # コンソール上に経過を表示するiterationの更新間隔

)

# samplingの詳細設定の解説も充実
## 指定可能な代表的なもの
# iter_warmup: burn-in数
# iter_sampling: butn-in後のsampling数
# thin: sampleを保存する間隔
# init: 初期値、init=function() list(theta=runif(1))のように
# threads_per_chain: 各chainに割り当てるCPUのthreads数


# 計算結果はtibble形式になっている
fit$summary() # 計算結果の要約、tibble形式
fit$sampler_diagnostics() # 収束診断
fit$draws() # sample値(次元:iteration, chain, parameter)、
# "draws_array"形式だが、matrix(fit$draws()[,,1])あるいはmatrix(fit$draws("theta"))のようにparameter名で指定し、parameter毎に取り出し可能


# しかし、以下によってrstanの形式に変換することも可能
stanfit <- rstan::read_stan_csv(fit$output_files())

# rstanの形式にできた方が都合が良いだろうが、rstanのインストールが全くできない場合はtibbleの方がありがたい場合もあるだろう。

# 推定結果を保存する場合、以下のようにするのが推奨されている。
fit$save_object(file = "fit.RDS")
# 再読み込みの場合
fit2 <- readRDS("fit.RDS")

2020年4月3日金曜日

当面要注意、mac版Rのアップデート

(2021.06.03追記)Apple SiliconでのRはRosetta2上でRが動作して高速といった書き込みを見掛けます。開発版Rの公開は始まっているようですが、安定版がいつ頃出るのか気になるところです。

(2020.08.29追記)クリーンインストール状態のR4.0.2をセットアップ試みました。CRANに、ソースからのパッケージコンパイルをする場合の方法が追記されていて(https://cran.r-project.org/bin/macosx/)、それ通りでrstan以外は適正にパッケージがインストールできるようになった模様。Xcode, command line tools, XQuartz, GNU Fortran 8.2をインストールせよとあります。注意点は、GNU Fortran 8.2はリンク先でfor Marvericsとなっていますが、Catalinaにインストールしても今のところ問題は起きていません(cf. https://thecoatlessprofessor.com/programming/cpp/r-compiler-tools-for-rcpp-on-macos/)。
Catalinaでrstanがうまく動作しない問題は深刻なようで、範囲も広くて把握できていないので、とりあえず放置しています。より問題の少なそうなrstan代替のパッケージ紹介もしておきました。

まとめると、上記のようにcranのガイドに従って手順を踏みさえすれば、rstan以外のインストール問題は解決しているようです。それでも以前のように、Rインストーラひとつで終わりとはいかなくなっているので、ハードルが上がってしまっており残念なことです。


*************(以下は、過去の参考履歴として残しておきます)****************************
(2020.04.14追記)CRANからRcpp1.0.4.6が通常版として配布され、ソースからのコンパイルが失敗する問題はある程度解決しているかもしれません。しかし依然、バイナリー版を見つけられない問題は残っていますし、ソースのコンパイルが通常の方法のみで行けない問題も残っています(tidyverseの必須パッケージの1つxml2など)。

現在(2020.03.31に確認)mac版のRをアップグレード後に、パッケージのインストールがうまくいかなくなるトラブルが生じるかもしれません。どれくらいの範囲、環境まで影響するかわからないのですが、CatalinaだけでなくEl Capitanでも発生しています。沢山の追加パッケージをインストールする必要のある方、とくにC++、Fortranコンパイル系のパッケージを利用している方は、しばらく様子を見た方がいいかもしれません。実態がわからず、復旧に相当な労力を要しました。通常ならばmac版Rのパッケージの再インストールなど、install.packages関数のリストのコピペ一発で済むはずでしたが、今回は日単位の時間を浪費する羽目に合いました。最近Linux版Rのセットアップをしたところでしたが、それよりもはるかに困難でした。

生じた問題は大きく分けて2つあります。

(1)install.packages関数がバイナリーパッケージを見つけられなかった。
バイナリーよりもソースが新しいversionである場合に、ソースパッケージが選択されることはこれまでも時折ありましたが、今回はCRANのアドレスは正しいのに、必ずソースパッケージを探しに行ってしまう現象が生じています。試しにネットワーク環境を変えてみても状況は変わらずで、また一方で同じネットワーク環境でも未アップデートのR(R3.5.1@El Capitan)では問題なくバイナリーを取得できました。エラーメッセージは以下のような感じに出てきました。

install.packages("パッケージ名") を実行すると、以下のようなエラーメッセージが出て、バイナリーをスキップしてソースファイルのインストールが試みられました:

警告:   リポジトリー https://cran.ism.ac.jp/bin/macosx/el-capitan/contrib/3.6 に対する索引にアクセスできません :
   URL 'https://cran.ism.ac.jp/bin/macosx/el-capitan/contrib/3.6/PACKAGES' を開けません 
 ソースパッケージ ‘パッケージ名’ をインストール中です 
 URL 'https://cran.ism.ac.jp/src/contrib/パッケージ名.tar.gz' を試しています 
(この後、ソースのコンパイルが始まる。つまりバイナリーの取得には失敗するけれど、ソースは取得できているので、CRAN自体には接続できている)

(2)大半のRパッケージでソースのコンパイルに失敗した。
強制的にソースを取得させられるのに、それがことごとく失敗しました。展開されたコード内に、'uuid_t'; 'uid_t'; といった大量のエラーが生じてコンパイルできない現象が生じました。

これらの現象は、El Capitan、クリーンインストール後のCatalina両方で発生しました(つまり私は原因がわからずに、いい機会だからとOSのクリーンインストールまでしてしまった!)。初めは、El CapitanのR3.6.1を3.6.3にアップしたときにこの症状が発生したので、いったんRを完全消去後に3.6.1を入れ直しました。しかし症状は改善せず。コンパイルに失敗しているので開発環境をクリアにする必要を感じて、この機会にCatalinaにアップしたら、まさかの再発に苦しめられたという顛末です。


****** 解決策 ******************************************************
(1)については、さしあたってCRANにある対象パッケージのバイナリーファイルのアドレスをコピーしてきて実行すれば、とりあえずインストールできました。多くのユーザにとって、これが最短の解決法だと思います。パッケージ名を入れ替えるだけなら簡単ですが、バージョンがわからないので、"cran パッケージ名"でググるなどして調べる必要があります(ターミナルで、R CMD INSTALLで行けないかと思ったけれど、これはソースしか取得できない模様)。私はこっちの解決法に気づくのが遅れて、(2)に取り組んでしまいましたが...。
# dplyrの場合の例、URLを用いたバイナリーからのインストールコード:
install.packages("https://cran.r-project.org/bin/macosx/el-capitan/contrib/3.6/dplyr_0.8.5.tgz”) 

(2)はより深刻です。心折れかけたところでu_riboさんに相談したところ、関連するトピックを見つけて教えてもらいました(https://github.com/RcppCore/Rcpp/issues/1060)。
これを読み解くと、問題として、Mac OSX SDKとRcppのバージョンが合っていない、Rcppに依存しているパッケージの大多数が影響を受けているようです。トピ主のパワーユーザがdplyrとhttpuvで試した限りでは、Mac OSX SDK 10.11の使用(ただし古いOSでしか使用不能とのこと)、あるいはMac OSX SDK 10.5とRcpp 1.0.4.4の組み合わせならば、とりあえずコンパイルに成功したとのこと。Rcpp 1.0.4.4はテスト版のようで、インストールするなら覚悟が必要な代物です。
こちらの問題の対処はRcpp 1.0.4.4だけでなく、複数の対処が必要です。重要な手順はいくつか判明していますが、効果があったかわからないものも含んでいます。それらを検証できる腕もないので、他の方の環境で、同じやり方で解決できるのかどうかは不明です。何の保証もなくお勧めできない最終手段だと思います。

(a)Rのインストールからやり直した方がよさそうなので、ターミナルから以下のコードでRを消去しました。"存在しない"とエラーが返ってくるものもありました。
rm -rf /Applications/R.app
sudo rm -rf /Library/Frameworks/R.framework
sudo rm /usr/bin/{R,Rscript}
sudo rm /usr/local/bin/R
sudo rm ~/library/R/

(b)Rのパッケージファイル(.pkg)のインストール(https://cran.ism.ac.jp/bin/macosx/)、現時点の最新版のR3.6.3を選択しました。Catalina用だけは分けられているので注意です。
なお、こちら(https://ryanhomer.github.io/posts/build-openmp-macos-catalina-complete)のおすすめに従い、カスタムインストールでTcl/TkとTexinfoを外しました。

(c)開発環境Command Line Toolsのインストール。ターミナルで以下のコードを実行、Gバイト単位のファイルを取得するので時間が掛かります。
sudo xcode-select --install

(d)ターミナルで以下のディレクトリを開き、SDKのversionを確認します:
open /Library/Developer/CommandLineTools/SDKs/ 
# 私のCatalina環境の場合、MacOSX10.15.sdkが既にインストールされていました。より下のversionの場合、別の対処が必要かもしれません(https://github.com/RcppCore/Rcpp/issues/1060)。

(e)(入れていなかったら)homebrewのインストール。ターミナルで以下のコードを打つ
ruby -e "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/install)"

(f)Tcl/TkとTexinfo、その他のインストール(homebrewユーザへは、bとは別途、こうした方がよいと勧められていました)
brew install tcl-tk texinfo
brew install llvm libomp
brew install gcc

(g)~/.R/Makevarsファイルの書き換え。こちらの4.の項目を実行。必要だったかどうか不明(https://thecoatlessprofessor.com/programming/cpp/r-compiler-tools-for-rcpp-on-macos/

h)clang7.0、gfortran6.1の.pkgファイルのインストール(これは随所で勧められていました)
https://cran.r-project.org/bin/macosx/tools/clang-7.0.0.pkg
https://cran.r-project.org/bin/macosx/tools/gfortran-6.1.pkg

(i)Rcpp1.0.4.6のインストール Rcpp1.0.4.4のインストール
以下のRコードで入りましたが、まだテスト版とのことなので、入れるなら覚悟して入れてください: 問題のあったversionのRcppはアップデートされ、通常のインストール方法でRcpp1.0.4.6がインストールされるようになっています。
install.packages("Rcpp", repos="https://rcppcore.github.io/drat")

(j)tidyverseの必須パッケージの1つxml2のコンパイルを試みる際に以下のようなエラーが頻発するので、これも実行。こちらを参考(https://github.com/r-lib/xml2/issues/223)。
ちなみに、ANTICONF ERRORはLinux版のRでのパッケージコンパイルでよく見かけますが、ここを見ると対処方法のヒントが得られやすいです。macでこれをやらなければならないかと思うと気が重いですが。

# ターミナルで、以下を実行
brew install pkg-config
# Rで、以下を実行
install.packages("xml2", configure.vars = c("INCLUDE_DIR=/Library/Developer/CommandLineTools/SDKs/MacOSX10.15.sdk/usr/include/libxml2"))

* installing *source* package ‘xml2’ ...
** package ‘xml2’ successfully unpacked and MD5 sums checked
Found pkg-config cflags and libs!
Using PKG_CFLAGS=-I/usr/include/libxml2
Using PKG_LIBS=-L/usr/lib -lxml2 -lz -lpthread -licucore -lm
------------------------- ANTICONF ERROR ---------------------------
Configuration failed because libxml-2.0 was not found. Try installing:
 * deb: libxml2-dev (Debian, Ubuntu, etc)
 * rpm: libxml2-devel (Fedora, CentOS, RHEL)
 * csw: libxml2_dev (Solaris)
If libxml-2.0 is already installed, check that 'pkg-config' is in your
PATH and PKG_CONFIG_PATH contains a libxml-2.0.pc file. If pkg-config
is unavailable you can set INCLUDE_DIR and LIB_DIR manually via:
R CMD INSTALL --configure-vars='INCLUDE_DIR=... LIB_DIR=...'
--------------------------------------------------------------------

(k)さらに、頻発する "make: gfortran: No such file or directory"というエラーに対処するため、~/.R/Makevarsにパスを追記をする。こちらを参考( https://github.com/merliseclyde/BAS/issues/1 )。これによって大多数のエラーが解消されました。
open ~/.R/Makevars
# これによって開かれたファイルに以下のパスを書き込んで保存。置き換えでなく追記です。
F77 = /usr/local/gfortran/bin/gfortran
FC = /usr/local/gfortran/bin/gfortran

(l)ここまで95%くらいのパッケージのインストールに成功しましたが、ここまで来て、パッケージのURLから直接取得する方法にたどり着いたので、これ以上の対処はせず。


(2)まったくひどい対処方法だと思うので、現状では個々のRパッケージをダウンロードしてきてインストールする(1)の方法が無難だと思います。

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",...) 

2017年3月30日木曜日

Rで機種依存文字の混じったデータファイルを読み込む方法

Rへ読み込む際に非常に厄介な機種依存文字(=環境依存文字)、データを作成する際には絶対に避けてほしい代物ですが、丸に数字など、日本のあらゆる文書で愛されており、根絶は非常に困難です。

見渡せる程度の小さなファイルならば手作業で特定し、検索・置換で除去するのですが、巨大なファイルの場合、いったいどんな文字が悪さをしているのか見つけることさえできない場合がありますね。

今回、到底手作業で対処できないサイズのファイルに取り組んでいた際にいい方法を思いついたので挙げておきます(小ネタです)。

1) まず、機種依存文字のファイルは.xlsxファイルとして保存しておきます(ここでは"データファイル名.xlsx"とします)。

2) 通常、read.csv関数などでデータファイルを読み込むところを、openxlsxパッケージのread.xlsx関数を使用します(要インストール、install.packages("openxlsx") などで)。

3) 以下で読み込むことができます、fileEncodingの指定も特に不要の様子。
data <- read.xlsx("データファイル名.xlsx")

cf. もちろんファイル出力にも対応しています(outputというデータフレームを出力する場合の例)write.xlsx(output, "出力ファイル名.xlsx")

Excelファイルを直接読むのはまだ抵抗はあるのですが、気にしなければ汎用性のある方法なので便利です。あとはこのパッケージがちゃんと存続してくれることを祈るばかりですが。なお、他にもExcelファイルを読めるパッケージは複数存在しているようですが、他のは基本的にJavaに依存しているようなので却下しました。

(2019.09.17 追記)
列名に漢字が混じる場合にフリガナが追加されるエラーが出ていて対処方法が分からず。代替関数として、readxlパッケージのread_excel関数の方が良いかもしれません。こちらならこのエラーは起こらず、またJava不使用です。tidyverseの一部のようなので、将来性もあるかもしれません。一方、読み込むとデータフレームではなくtibleという新形式になります。tible非対応のパッケージもまだ多いので、通常のデータフレームに変換するには、as.data.frame(tible_no_data)としてやればよいです。

cf. 機種依存文字が混じっているファイルをread.csvで読もうとすると、"In scan(file, what, nmax, sep, dec, quote, skip, nlines, na.strings,  : 入力コネクション 'ファイル名.csv' に不正な入力がありました "のようなエラーメッセージが出る。
Rじゃなくとも、メールの送受信でも文字化けの温床となります。データ解析に限らず、トラブルの原因にしかなりません。