2024年1月2日火曜日
情報量規準による多重比較のRコード
AICなどの情報量規準を用いた多重比較の方法（Dayton 1998）を実装するRコードの紹介です。以前、共著論文で使用する機会があったのですが、久しぶりに掘り返して、だいぶスッキリしたコードへと改定できたので公開することにしました。手作業だと面倒くさかったのでだいぶ手間が省けるはずです。
非常によく見かける解析方法に、ANOVAの後で多重比較（post-hoc tests）というコンボがあリます。いわゆる多重比較の方法は検定（頻度論）に基づいており、Type I errorの回避に特化しているのが気になっていました。一方で情報量規準に基づくモデル選択による解析を行うケースも多いですが、旧来の多重比較（頻度論）も併用するとダブルスタンダードになってしまいます。実際に査読で突っ込まれたこともあります。また多重比較は正規分布を仮定できるデータにしか適用できないことも制約になります。この点についてはHolmの方法など原理的に分布型に制約されない方法もありますが、組み合わせ数が多くなると検出力が低下する問題もあります。
こういった問題をすっきりと解決できる、情報量規準に基づいたモデル選択による多重比較の方法があります。以前、北大・久保先生のスライドで知った手法ですが（ https://kuboweb.github.io/-kubo/ce/2004/index.html ）、下記のDayton（1998）が同様の手法（PCIC）を発表しており、根拠論文として参照可能です。特に手法名が提唱されているとは知らずに使っている例も多いかもしれませんが。また、名前が知られていないからこそ、統計ソフトにもRパッケージでも実装されてこなかったのかもしれません。
Dayton CM (1998) Information criteria for the paired-comparisons problem. The American Statistician, 52(2):144-151. DOI: 10.1080/00031305.1998.10480554（ResearchGateからも取得可能）
下記の10数行のコードで実装できました。lm, glm, glm.nbなどのモデル式を下記のpcic関数に食わせれば実行できるはずです。
require(partitions)
# 組み合わせを表現する関数setpartsを使うのに必要なパッケージ、要インストール
# pcic関数の定義、ここから：
pcic <- function(Model, Data=Model$model, Expl=names(Model$model)[2]) {
# pcic関数の定義、ここから：
pcic <- function(Model, Data=Model$model, Expl=names(Model$model)[2]) {
nLevels <- length(levels(Data[,Expl])) # グループのリスト
groupPTs <- apply(setparts(nLevels), 2, function(x) as.numeric(factor(x, levels=unique(x))))
# グルーピングの全パターン、setpartsの標準ではb,a,aなどが出るので正順に並べ替え
groupPTs <- apply(setparts(nLevels), 2, function(x) as.numeric(factor(x, levels=unique(x))))
# グルーピングの全パターン、setpartsの標準ではb,a,aなどが出るので正順に並べ替え
groupPTA <- apply(matrix(letters[1:nLevels][groupPTs], nrow=nLevels), 2, paste, collapse=" ")
・lm, glm, glm.nbなどのモデルに適用可能です。それ以外でもModel$modelを実行してモデルに使用しているデータ（y, x）が返ってくるならば使用可能、offset項にも対応。
実行例：
data(iris) # Sepal.Lengthが種によって異なるかについて多重比較を実行
Model <- glm(Sepal.Length ~ Species, data=iris)
require(partitions)
# アルファベット変換&結合
AICs <- c(AIC(update(Model,
AICs <- c(AIC(update(Model,
reformulate(c(setdiff(names(Data)[-1],Expl),1), response = names(Data)[1]))), # 帰無モデル
sapply(2:ncol(groupPTs), function(i) { # 帰無モデル以外
Data$groupLv <- factor(letters[1:nLevels][groupPTs[,i]][Data[,Expl]], levels=letters[1:nLevels])
Data$groupLv <- factor(letters[1:nLevels][groupPTs[,i]][Data[,Expl]], levels=letters[1:nLevels])
# グループ分けの変数をDataに追加
AIC(update(Model, reformulate(c(setdiff(names(Data)[-1],Expl),"groupLv"),
AIC(update(Model, reformulate(c(setdiff(names(Data)[-1],Expl),"groupLv"),
response = names(Data)[1]), data=Data)) } ) ) # 各モデル
comps <- data.frame(AIC=AICs, deltaAIC=AICs-min(AICs), Grouping=groupPTA)
print(paste(c(paste(levels(Data$Species), collapse=", "), "are labeled in this order.",
Expl, "sharing the same alphabet do not differ in the", names(Data)[1]), collapse=" "))
comps <- data.frame(AIC=AICs, deltaAIC=AICs-min(AICs), Grouping=groupPTA)
print(paste(c(paste(levels(Data$Species), collapse=", "), "are labeled in this order.",
Expl, "sharing the same alphabet do not differ in the", names(Data)[1]), collapse=" "))
return(comps[order(comps$AIC),]) }
# 関数定義ここまで（これを一度R上で実行するとカスタム関数として利用可能になる）
# 関数定義ここまで（これを一度R上で実行するとカスタム関数として利用可能になる）
・lm, glm, glm.nbなどのモデルに適用可能です。それ以外でもModel$modelを実行してモデルに使用しているデータ（y, x）が返ってくるならば使用可能、offset項にも対応。
・注意点として、グルーピングのための説明変数は "~"のすぐ後ろの変数と見なす。異なる時は文字で指定する（この例の場合は"Species"）
・Data: Model$model # モデルに使用しているデータ：y, x。これで大丈夫なはずだが、自前で指定も可能
実行例：
data(iris) # Sepal.Lengthが種によって異なるかについて多重比較を実行
Model <- glm(Sepal.Length ~ Species, data=iris)
require(partitions)
pcic(Model)
[1] "setosa, versicolor, virginica are labeled in this order. Species sharing the same alphabet do not differ in the Sepal.Length"
AIC deltaAIC Grouping
5 231.4520 0.00000 a b c
4 265.6359 34.18393 a b b
3 295.6778 64.22581 a a b
1 372.0795 140.62758 a a a
2 373.1310 141.67904 a b a
一番上がベストモデルです。setosa: a, versicolor: b, virginica: c、のように異なるアルファベットで識別されており、種によってSepal.Lengthが異なることが示されました。このアルファベットによる識別は一般的な多重比較でよく見かける表現なので、多くの人が馴染みがあるかと期待します。
結果を表形式で取り出すには、こんな感じでOKです。
write.csv(pcic(Model), "保存場所/ファイル名.csv")
最後に、post-hoc testの用途について。実験系など個々のグループの意味が明確な場合に適した手法だと思いますが、その場合は頻度論の方が適していそうです。ここまで書いておいて何ですが、個人的にはpcicが必要なケースってそんなに多く無いと思っています。生態学系でよく見かける使用例の多くは回帰分析の方が適しているケースだし（説明変数が連続変数になっている等）、2way ANOVAなどの方が適しているケースも多いように感じます。何の仮定も置かずにpost-hoc testをやってみて、その結果を見てから演繹的に考察をするような使い方をしようとしていたら、少し立ち止まって考え直してみませんか。仮説を立ててそれを検証するための解析を行うという設計さえしていれば、より適切な解析方法が見つかるはずです。
cf. lme4は今のところ非対応です。Data=Model@frameでデータを引っ張って来れそうですが、モデルのアップデートのところの改定がすぐにはできそうもなく。
2023年4月9日日曜日
Rの常用速度を新旧macで比較：Apple Silicon M1 Max vs. Intel Core i9
今更ながら、macでのRの計算速度をIntel版とApple Siliconとで比較してみました。先行するテスト結果は検索するといくつか見つかるのですが、自分にとって普段よく使いそうな実用的用途で試してみたいと思いました。（2023年04月25日追記：Intelも第12世代でCPU速度は飛躍的に向上して、以下の比較機種とは状況は異なっている模様）
比較をした計算環境は次の通りです。Rは登稿時の最新 ver. 4.2.3に揃えました。ハードウェアはintelとapple siliconとで実質1年半の差なので、ほぼCPUの差で比較できているかと思います。
（1）macbook pro 14-2021（以下、MBP_M1mと略します）
Apple Silicon M1 Max 10コアCPU、32コアGPU 64GBメモリ
（2）macbook pro 16-inch 2019（以下、MBP_Intelと略します）
Intel Core i9 2.4GHz 8コアCPU 64GBメモリ
（3）新旧macだけでもよかったのですが、Ubuntuを載せたカスタムPCも比較に加えました。
カスタムPC 2019年製（OS: Ubuntu 22.04.1 LTS）（以下、Intel_Ubuntと略します）
Intel Core i9-9980XE Extreme Edition 3.0-4.4GHz 18コアCPU 128GBメモリ
計算は次の4つのケースで実行し比較しました。概ねシングルコアの速度を試す形になっていますが、マルチコアを活用できるケースならば、単にコア数の多いマシンを使えばいいんじゃないかと思うので、制約の掛かりやすそうな事例をわざと選んでいます。
いわゆるMCMC等の計算過程で動いているマルコフ連鎖、値の更新が必要なので一括処理できない計算の例の1つかと。
x <- matrix(0, 10^6, 10)
# 10^6回分のマルコフ連鎖、10本の鎖を同時に計算
system.time( for(r in 2:(nrow(x))) x[r,] <- 0.9*x[r-1,] + rnorm(10) )
### Case 2 raster_distance: raster画像の値あり全地点からの距離計算
require(raster)
require(maptools)
data(wrld_simpl)
rgrid <- raster(ex=extent(c(-180,180,-90,90)), res=0.5)
rwrld <- raster::rasterize(wrld_simpl, rgrid)
system.time(rdst <- raster::distance(rwrld))
0.5度解像度の全世界の陸地からの距離を計算、かなり時間のかかるケース。
（ちなみにwrld_simplは気軽に呼び出しやすい粗い世界地図なのですが、{maptools}が2023年いっぱいで提供終了とのこと。高機能な{rnaturalearth}もよいですが...）
### Case 3 raster_linear_model: raster stackへの回帰モデル適用
重ねたrasterの全値について回帰モデルを適用、calcのhelpにあるコードを高解像度にしただけのもの。単純な回帰モデル（lm）に限定（逆行列を利用した方法も紹介されていますが、応用しづらいので省きます）。
help(calc)
r <- raster(nrow=100, ncol=100) # ここを10倍にした
s1 <- lapply(1:12, function(i) setValues(r, rnorm(ncell(r), i, 3)))
s2 <- lapply(1:12, function(i) setValues(r, rnorm(ncell(r), i, 3)))
s1 <- stack(s1)
s2 <- stack(s2)
# regression of values in one brick (or stack) with another
s <- stack(s1, s2)
# s1 and s2 have 12 layers; coefficients[2] is the slope
fun <- function(x) { lm(x[1:12] ~ x[13:24])$coefficients[2] }
system.time(x1 <- calc(s, fun))
### Case 4: Stan (cmdstan)
コマンドライン版のStanで単純な状態空間モデルの計算、chain毎に異なるスレッドを割り当てているので、他の例と異なり、少しだけマルチコアの恩恵があるはず。コードは長いので、後に載せておきます。
こちらが結果になります。左に行くほど高速です。
概ねMBP_M1mはMBP_intelの2倍弱ほどの速度があることが分かりました。計算時間の掛かるraster_distanceで2.33倍と最も差が大きいのは嬉しい結果です。マルチコアの恩恵があるStanでは1.5倍とさほどの差になっていないのも予想通りではあります。またOSの違いはあまり速度に関係ないようでした。
なお、マルコフ連鎖については、2012年のマシンで5.457秒という手元の記録が残っているので（Late 2012 iMac Intel Core i7 3.4GHz クアッドコア 32GB）、倍速いというのがどれほどの進化かよく分かるかと思います。
ちなみに、関連してCase 2, 3ではrasterパッケージの後継？terraパッケージでも比較をしてみたのですが、ほぼ変わりはなかったです。もしかすると連繋が進んでいたりして内部でterraを利用しているかもしれませんが、あくまで想像です。
もうM2が登場していますが、今回のようなRのシングルコア計算での比較ならば、おそらく速度は10%向上程度かなと見ています。
その他、Apple SiliconはGPUコアや機械学習用のニューラルエンジンを積んでいるので、用途次第では次元の違う速度になる期待もあります。高速なGPUを使用可能なRパッケージってまだあまり聞かないようには思いますが。
あとはSSDがかなり高速になっているそうで、物理メモリが不足した時にSSD領域を仮想メモリとして使用しても、かつてほどの速度低下はしないみたいで、もうさほどメモリに拘らなくてもよいのかもしれません。
# Case 4のcmdstanの計算コード：
stan_program <- "
data {
int T;
vector[T] Y;
}
parameters {
real mu_Y0;
vector[T] mu_Y;
real<lower=0> sigma;
real<lower=0> SIGMA;
}
model {
// state model
mu_Y[1] ~ normal(mu_Y0, sigma);
for(t in 2:T) {
mu_Y[t] ~ normal(mu_Y[t-1], sigma);
}
// observation model
for(t in 1:T) {
Y[t] ~ normal(mu_Y[t], SIGMA);
}
}"
modfile <- write_stan_file(stan_program)
mod <- cmdstan_model(modfile) # 一時ファイルへ書き出し
data(Nile)
data <- list(Y=Nile, T=length(Nile))
system.time(fit <- mod$sample(
data = data,
chains = 4,
parallel_chains = 4, # 並列化するchainの数
iter_warmup = 5000, #burn-in数
iter_sampling = 20000 #butn-in後のsampling数
)
)
2021年9月30日木曜日
問題を起こすプロセスを定期的に強制終了する
macの動作を重くしたり不安定にさせたりするプロセス（アプリケーション含む）はなるべく終了させるか排除したくなるものです。純正のものではWindowServer, Falcond, mdworkerなどが重いですが、システムに必要そうなので我慢です。しかし、非純正のプロセスならいっそ落としてしまえば動作が改善する場合があるかもしれません。
このところmacが非常に不安定になっていて、Zoom会議で落ちる、Rで負荷の掛かる処理をすると落ちる、フォルダを移動しただけでもなど、処理内容が重い方が落ちやすいけれど軽くても構わず落ちる（カーネルパニック）現象がひどい時は1〜2日に一度くらいの頻度で起きるという困った状況に陥っていました。これではまったく使い物になりません。もしやと思い、心当たりの非純正アプリケーションを朝・夜にアクティビティモニタから強制終了させてみたところ、それらのトラブルは全て解消されました。残念ながら該当のアプリを削除するわけにはいかなさそうなのですが、この定期的なリフレッシュを自動で実行する仕組みを構築してみました。以下のような手順です。Automatorとカレンダー.appを使用します。
1. まず対象プロセスを強制終了させるアプリケーションを作成する。
1-1.「アプリケーション」から「Automator.app」を立ち上げる。
1-3. 「アクション」の中から「シェルスクリプトを実行」を選び、右の枠の中へとドラッグアンドドロップする。
1-4. シェルスクリプトの入力は、「pkill プロセス名」とする。これはプロセス名のプロセス（アプリケーションである場合もある）を強制終了させるコマンドである。「ユーティリティ」の「アクティビティモニタ」にリストアップされるプロセス名で記す。プロセス名は先頭から一部の文字で構わない。また大文字と小文字を区別する必要がある。例えば、対象としたいプロセス名が「Hogehoge1」と「Hogehoge2」の場合、「pkill Hoge」で両方が対象となる。「pkill hoge」だと対象と見なされず何も起こらない。
1-5. 「ファイル」メニューから「アプリケーション」として保存する（ワークフローではなく）。保存場所はどこでも良さそう。一度、手動で立ち上げて、ちゃんと狙い通り動作（対象プロセスを強制終了）するか確かめておいた方が良いだろう。
2. 作成したアプリケーションを定期的に実行させる仕組みを作る。
6. 「カレンダー.app」の予定を作成する。アプリを実行したい時間を設定する。さらに「通知」の仕方として「ファイル（イベントの開始時刻）を開きます」を選択、開始時刻リストの下にある「カスタム」を選択する。
7. 「ファイルを開く」対象がデフォルトでは「カレンダー」のところを、「その他...」を選択し、先ほど作成したアプリケーションを選択する。
8. 「繰り返し」を毎日にすれば、決まった時刻に実行される。私の場合、朝用と夜用を作成し、1日に2回リセットさせるように設定した。カレンダー.appは落としておいても動作するようです。
難点といえば「カレンダー」に毎日この定期リセットのための予定が書き込まれてしまうことです。ただし予定のグループを専用に設けて非表示にしておけば見えなくすることができます。また、通知機能を利用しているので、定時のリセットの度に通知されるのも気になると言えば気になりますが、ちゃんと機能していることを意識はできます。
2020年12月20日日曜日
tidyverse版: 複数生物種の個体数集計データを縦に伸ばすor縮める（集計データ← →データフレーム変換）
（以前の記事をtidyverseで作り直してみました。まだ旧来のRコードの方が馴染みがあるので、tidyverseに慣れている方からすると不格好かもしれません。）
多数の種が含まれる生物群集データを解析する時、エクセルなどにデータを整理していると、各種をインデックスにするか要因のひとつとして整理するか、ケースによって変わってくるかと思います。とくに各種を応答変数とするか説明変数に使用するかあたりで必要になってくるかと。完全なデータフレーム型にしておいた方が統計解析はやりやすいけれど、膨大な行数となってエクセルの限界にあっさりと到達したりもするし…（cf. エクセルのvlookup関数）。
こんな時に、一部集計データとデータフレーム型とを自在に一発変換できたら効率がよいと思い、やり方をまとめてみました。正直、入力の手間が一番省けるのは部分的には集計されたデータだったりしますよね（下記のdata1のような形）。
データの整形にはdplyrパッケージの関数などを使用します。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
########################################################
# data1をデータフレーム型に縦に伸ばす場合
data2 <- data1 %>% pivot_longer(c(-year, -site, -depth), names_to="species", values_to="abundance")
> data2 # 出力してみます
########################################################
# もう一段階、伸ばしてみましょう。一個体あたり一行というデータセットへ。多項ロジットなど、カテゴリカルな解析に便利そうです。
data3 <- select(data2, -abundance) %>% slice(unlist(map2(1:nrow(data2), data2$abundance, rep)))
> data3
########################################################
# data3を一つ前の状態に戻してみます（観測地点あたりで集計、といったところ）
data4 <- data3 %>% count(year, site, depth, species, name="abundance")
data5 <- data4 %>% pivot_wider(names_from=species, values_from=abundance, values_fill=0)
# pivot_widerは縦長データから横長データへ変換する、見出しとなるspeciesと対応する値abundanceを展開する。生じるNAはvalues_fill=0で置き換えた。
> data5 # 出力してみます。
########################################################
data6 <- tibble(site=paste0("s", 1:10), month=rep(c(1:2), each=5),
species=c("スズメ", "","スズメ、ヒヨドリ、シジュウカラ", "ムクドリ、スズメ", "ヒヨドリ、スズメ", "スズメ、キジバト", "スズメ、ムクドリ", "スズメ、ヒヨドリ", "", "ムクドリ"))
# 例えばこんなかんじのデータセット。空欄さえあります…
# 注：このような日本語データのでの入力の場合、普通に読み込もうとすると文字化けしやすいので、readxlパッケージのread_excel関数などで読み込むのがトラブルが少なくていいと思っています。csvでは機種依存文字やmac/windows間での文字コード問題が混乱のもとなので。
# こんなデータでも次のようにすれば、通常の集計データに変換できます。data1の形態へ変換してみましょう。
sptabs <- data6 %>% # data6をベースに集計作業
data7 # 出力してみます（出来上がり）
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
########################################################
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
########################################################
data4 <- data3 %>% count(year, site, depth, species, name="abundance")
# year, site, depth, species毎にデータ頻度をカウントし、カウント結果をabundance変数として追加
> data4
# sp3 のゼロ個体データだけは消えてしまったが、当然と言えば当然
########################################################
# data4の種数部分を集計表状に横に伸ばす（最初の状態＝data1の形にする場合）
# 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
########################################################
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 "ムクドリ"
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の発表者ツールの左上に、ツールのボタンが並んでいますが、「ディスプレイの入れ替え」をクリックします。
すると、メイン画面が全画面のプレゼン、外部ディスプレイが発表者ツールに入れ替わるはずです。（2023追記：いっそここで「スライドショーの使用」を選び、両方のディスプレイにスライドショーが表示されるようにするのも良い解決方法かと思います。相手方を待たせるよりはいいでしょう。）
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日追記 → (削除ここまで)2021年09月30日現在、通常インストールの最新版で十分にうまく行くようになった模様）
# おそらくもう不要→ # 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) # 一時ファイルへ書き出し、この方が運用しやすいと思う
# あとは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 = 4, # 並列化するchainの数
refresh = 500 # コンソール上に経過を表示するiterationの更新間隔
)
# samplingの詳細設定の解説も充実## 指定可能な代表的なもの
# iter_warmup: burn-in数
# iter_sampling: butn-in後のsampling数
# thin: sampleを保存する間隔
# init: 初期値、init=function() list(theta=runif(1))のように
# 計算結果は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())
# 推定結果を保存する場合、以下のように.RDSで保存するのが推奨されている。
fit$save_object(file = "fit.RDS")
# 再読み込みの場合
fit2 <- readRDS("fit.RDS")
#（2021年12月31日補足）個々のchainに1つのCPU threadを割り当てる並列化は上記のようにparallel_chainsの指定をするだけなので容易だが、chainあたり複数のthreadを割り当てるには、モデルコード自体にも工夫が要るようだ。GPUの利用なども可能なようだが、いずれにしてもモデル構造自体の制約が大きく、またN数が非常に大きいなどのビッグデータには効果が大きいが、生態学でよく用いるような複雑なプロセスモデルでは恩恵が少なそうである。 https://mc-stan.org/docs/2_28/cmdstan-guide/parallelization.html
2020年4月3日金曜日
（2021.10月追記：もう大丈夫そう）当面要注意、mac版Rのアップデート
（2021年10月01日追記）R4.1.1へアップグレードした様子では2020年に発生していた下記の数多くの問題は既に解消されたように見えます。またApple Silicon用のRバイナリーパッケージは既にCRANからダウンロード可能になっているのも確認しました。
（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のインストール
（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）の方法が無難だと思います。
登録:
コメント (Atom)